amuse_tutorial_brussel2014: exercise1.py

File exercise1.py, 3.3 KB (added by jilkova, 4 years ago)
Line 
1# -*- coding: utf-8 -*-
2# <nbformat>3.0</nbformat>
3
4# <markdowncell>
5
6# # AMUSE units: the blackbody spectrum
7#
8# this is a small exercise to familiarize yourself with AMUSE units usage, learn to write simple AMUSE utility functions.
9
10# <markdowncell>
11
12# ## (1) first some definitions:
13#
14# 1. check that you understand the following definitions
15# 2. there is at least one bug below - if you dont find it don't worry about it just yet!
16# 3. what are the differences or similarities between an AMUSE unit, quantity and constant?
17
18# <codecell>
19
20import numpy
21
22from amuse.units import units,constants
23
24pi=numpy.pi
25e=numpy.e
26kB=constants.kB
27h=constants.h
28c=constants.c
29Ry=constants.Rydberg_constant
30sigma=constants.Stefan_hyphen_Boltzmann_constant
31
32def B_nu(nu,t):
33  return 2*h*nu**4/c**2 * 1./ (e**(h*nu/kB/t)-1)
34
35def B_lambda(l,t):
36  return 2*h*c**2/l**5*1./(e**(h*c/(l*kB*t))-1)
37
38def energy(nu):
39  return constants.h*nu
40 
41def freq(energy):
42  return energy/constants.h
43
44def freq_from_wavenumber(k):
45  return c*k
46
47def wavelength(nu):
48  return c/nu
49
50def freq_from_wavelength(l):
51  return c/l
52 
53def wiens_lambda_max(T):
54  b = 2897768.6 | units.nano(units.m)* units.K
55  return b/T
56
57def wiens_T_from_lambda_max(l):
58  b = 2897768.6 | units.nano(units.m)* units.K
59  return b/l
60   
61def total_bolometric_flux(T):
62  return sigma*T**4 
63
64# <markdowncell>
65
66# ## (2) using the above, calculate the following:
67#
68# 1. the wavelength of the maximum for a T=10000K blackbody spectrum
69# 2. the bolometric flux for a T=10000K blackbody spectrum
70# 3. the total luminosity of the sun (Teff = 5780 K), in units.LSun!
71
72# <codecell>
73
74pass
75
76# <markdowncell>
77
78# ## (3) energy flux calculation
79#
80# 1. what does the energy_flux function below do?
81# 2. what are the T, lowfreq, N arguments?
82# 3. calculate the flux for some values of T (like 5500 K, 10000K, 50000K), try to convert to W/m**2, check against total bolometric flux above,
83# 4. (spoiler alert) if you had spotted the bug: repeat 3. using the original buggy version of B_nu (with nu^4 instead of nu^3)..
84# 5. consider the numpy.trapz call: observe the use of .number and .unit, can you understand what this does?
85# 6. make a copy of the function and try to remove the .number and .unit - does the function still work? Can you think of a reason to use the form given here? (hint: what are b and b.number?)
86
87# <codecell>
88
89def energy_flux(T,lowfreq=0.|units.s**-1,N=10000):
90  nu=(numpy.arange(1.,N+1))/N*(kB*T)/h*25.+ lowfreq
91  b=pi*B_nu(nu,T)
92  return numpy.trapz(b.number,x=nu.number)| (b.unit*nu.unit)
93
94# <codecell>
95
96pass
97
98# <markdowncell>
99
100# ## (4) the photon flux
101#
102# 1. write a function that calculates the photon flux, with an optional lower frequency bound
103# 2. what would the lower bound have to be to calculate the ionizing flux?
104
105# <codecell>
106
107def photon_flux():
108    pass
109
110# <markdowncell>
111
112# ## (5) stellar evolution
113#
114# 1. instantiate a stellar evolution code (e.g. SSE or SeBa), and generate a stellar model for a 30. MSun star
115# 2. print out its attributes
116# 3. write a function that calculates the luminosity in ionizing photons Lion for a star (assuming of course it radiates like a blackbody!)
117# 4. make a plot of Lion vs time
118
119# <codecell>
120
121from amuse.community.sse.interface import SSE
122from amuse.community.seba.interface import SeBa
123from amuse.datamodel import Particle
124