# -*- coding: utf-8 -*- # 3.0 # # # AMUSE units: the blackbody spectrum # # this is a small exercise to familiarize yourself with AMUSE units usage, learn to write simple AMUSE utility functions. # # ## (1) first some definitions: # # 1. check that you understand the following definitions # 2. there is at least one bug below - if you dont find it don't worry about it just yet! # 3. what are the differences or similarities between an AMUSE unit, quantity and constant? # import numpy from amuse.units import units,constants pi=numpy.pi e=numpy.e kB=constants.kB h=constants.h c=constants.c Ry=constants.Rydberg_constant sigma=constants.Stefan_hyphen_Boltzmann_constant def B_nu(nu,t): return 2*h*nu**4/c**2 * 1./ (e**(h*nu/kB/t)-1) def B_lambda(l,t): return 2*h*c**2/l**5*1./(e**(h*c/(l*kB*t))-1) def energy(nu): return constants.h*nu def freq(energy): return energy/constants.h def freq_from_wavenumber(k): return c*k def wavelength(nu): return c/nu def freq_from_wavelength(l): return c/l def wiens_lambda_max(T): b = 2897768.6 | units.nano(units.m)* units.K return b/T def wiens_T_from_lambda_max(l): b = 2897768.6 | units.nano(units.m)* units.K return b/l def total_bolometric_flux(T): return sigma*T**4 # # ## (2) using the above, calculate the following: # # 1. the wavelength of the maximum for a T=10000K blackbody spectrum # 2. the bolometric flux for a T=10000K blackbody spectrum # 3. the total luminosity of the sun (Teff = 5780 K), in units.LSun! # pass # # ## (3) energy flux calculation # # 1. what does the energy_flux function below do? # 2. what are the T, lowfreq, N arguments? # 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, # 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).. # 5. consider the numpy.trapz call: observe the use of .number and .unit, can you understand what this does? # 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?) # def energy_flux(T,lowfreq=0.|units.s**-1,N=10000): nu=(numpy.arange(1.,N+1))/N*(kB*T)/h*25.+ lowfreq b=pi*B_nu(nu,T) return numpy.trapz(b.number,x=nu.number)| (b.unit*nu.unit) # pass # # ## (4) the photon flux # # 1. write a function that calculates the photon flux, with an optional lower frequency bound # 2. what would the lower bound have to be to calculate the ionizing flux? # def photon_flux(): pass # # ## (5) stellar evolution # # 1. instantiate a stellar evolution code (e.g. SSE or SeBa), and generate a stellar model for a 30. MSun star # 2. print out its attributes # 3. write a function that calculates the luminosity in ionizing photons Lion for a star (assuming of course it radiates like a blackbody!) # 4. make a plot of Lion vs time # from amuse.community.sse.interface import SSE from amuse.community.seba.interface import SeBa from amuse.datamodel import Particle