;__________________________________________________________________________ ; ; A fairly successful yet simple attempt at modeling a sun-like star ; ; The equations of stellar structure are standard 1st order differential ; equations, independent variable is r. We assume radiative ; transfer of energy throughout a star. ; ; The model has a specially crafted opacity vs. radius function, which ; begins near the center at the value of Thomson (electron) opacity of ; 0.4 cm^2/g. In the outer layers opacity kappa is chosen so that a total ; mass at one solar radius is very close to 1 solar mass. This opacity ; increases almost thousand-fold between the center and the surface. ; ; Another function that is chosen judiciously in order to obtain a reasonable ; model as a result of single outward directed integration of the set of ; structure equations, is the efficiency of nuclear reactions (energy per unit ; time emitted by unit mass of gas, or epsilon). It is chosen as the total ; luminosity divided by the core mass (assumed to be 0.15 Msun), and as the ; luminosity grows in course of integration, it is not allowed to ever exceed ; the Lsun value by appropriately reducing epsilon, when L(r) approaches Lsun. ; ; That's about as much is needed to succeed in stellar structure modeling ; in a very simple way. The model is still not that of our real sun, because ; opacity, nuclear reaction rate and energy transport are treated only very ; approximately (we neglect convection zone). ; Yet we obtain a very reasonable run of pressure, density and ; temperature inside our star. The central temperature is somewhat overestimated ; (17 instead of 15.7 mln K), but the central density and surface densities are ; close to real values. In particular, rho(0) = 162 g/cm^3, coinciding ; with fuller models. ; ; Students are strongly encouraged to rewrite the script in Python, Mathematica, ; or any convenient language, and modify the program, e.g., to see in ; practice how sensitive the results are to details of modeling eps(r) and ; kappa(r). ; A free implementation of IDL exists, named GDL - it can be installed & used to ; run the present script. ; ; P. Artymowicz, 24.09.2015 ;__________________________________________________________________________ ;------definition of constants, letter 'd' in the values defines double ; precision numbers, while the single precision has 'e' sig = 5.670367d-8 ; [W m^-2 K^-4] a = 7.565767d-16 ; [J m^-3 K^-4] a = 4*sig/c k_B = 1.38064852d-23 ; [m^2 kg s^-2 K^-1] mu = 1.18 ; mean mol. weight of solar gas m_H = 1.6735d-27 ; [kg] G = 6.674d-11 ; [N⋅m2/kg2] const = k_B/(mu*m_H) Rsun = 696300d3 ; [m] Msun = 1.989d30 ; [kg] Lsun = 3.846d26 ; [W] kappa_Th = 0.04d-4/1d-3 ; [m^2/kg], 0.4 cm^2/g eps_0 = Lsun/(0.15*Msun); [W/kg] rough estimate of epsilon, core mass ~15%*Msun coeff = 600 ; multiplier needed to model opacity increase near the surface ;---------starting conditions at the center L = 0. m = 0. T0 = 17e6 ; [K] rho0 = 162d3 ; [kg/m^3] 162 x water density T = T0 rho = rho0 pr = a*T^4/3 ; p_rad(0) pg = const*T*rho ; p_gas(0) p = pr + pg ; p(0) p0 = p kappa = kappa_Th ;------start at the center, and integrated from r=0 to r=Rsun N = 300 ; number of shells dr = Rsun/N ; thickness of one shell struct = fltarr(14,N+1) ; this defines an array holding all variables struct(0,0) = 0. ; r struct(1,0) = 0. ; m struct(3,0) = T ; T/T(0) struct(4,0) = p ; total central pressure p(0) struct(5,0) = pr ; p_rad(0) struct(6,0) = pg ; p_gas(0) struct(7,0) = rho ; rho(0) struct(10,0) = L ; L(0) struct(11,0) = 0. ; I(0) struct(12,0) = eps_0 ; struct(13,0) = kappa ; print,' T0=',T0/1d6, ' mln K' print,' rho0=',rho0/1d3, ' g/cm^3' ;------integrate equations of structure shell by shell and store the results for i = 1,N do begin ;-----compute last p_gas and current volume of a shell r = dr*(i-0.5) surface = 4*!pi*r*r ; surface area of a shell of radius r dvol = surface*dr ; [m^3]; volume of a shell of radius r and thickness dr ;------equation of state (ideal gas law) pg = const*T*rho ; p_gas(r) ;------mass eq dm = dvol*rho m = m + dm ; m(r) ;-----luminosity eq. eps = eps_0 *((1.-L/Lsun)^0.75>0.); make sure L does not exceed Lsun L = L + eps*dm ;------dT/dr radiative balance eq. ; opacity and radiation pressure of photon gas => ; dT/dr = -L*3*kappa*rho/(16*sig*T^3*4*!pi*r^2) ; kappa is somewhat arbitrarily modeled to increase as a function of radius ; like this: kappa = kappa_Th *(1d0 + coeff*((r/Rsun)^1.15 +60*(r/Rsun)^4.0)) dTdr = -L*3*kappa*rho/(16*sig*T^3*surface) ; standard radiative condition T = (T + dTdr*dr) > 5780. ; lower limit on T is eff. surface T ;------hydrostatic pressure eq pr = a*T^4/3 ; next p_rad(r) ; p = pr + pg dPdr = -G*m/r^2 *rho ; hydrostatic balance of total p & gravity p = p + dPdr*dr ; next p(r) pg = p - pr ; next p_gas(r) rho = pg/(const*T) ; next rho(r) ; ideal gas EOS ;-----catch wrong solutions and stop when found if (p*pr*pg LT 0.) then stop ;-----store variables for later plotting struct(0,i) = r ; r struct(1,i) = m ; m struct(3,i) = T ; T(r) struct(4,i) = p ; p(r) struct(5,i) = pr ; pr(r) struct(6,i) = pg ; pg(r) struct(7,i) = rho ; rho(r) struct(10,i) = L ; L(r) struct(11,i) = L/(4*!pi*r^2) ; I(r) struct(12,i) = eps ; struct(13,i) = kappa ; ;---------print print,i,' r T rho',float(r/Rsun),float(T/T0),float(rho/rho0),' eps',float(eps) print,i,' p pr pg',float(p/p0),float(pr/p0),float(pg/p0),' L',float(L/Lsun) print,'r T rho',float(r/Rsun),float(T/T0),float(rho/rho0),' eps',float(eps) endfor ;-----integration is done, the rest is just ;-----IDL plotting of results: this part produces plot on the screen x = struct(0,*)/Rsun lab = string(coeff) plot,x,struct(7,*)/rho0,xtitle='r/R', ytitle='normalized values',$ title='solar model with T(0)=17 mln K, rho(0)=162 g/cm3',charsize=2.3 oplot,x,struct(1,*)/Msun oplot,x,struct(3,*)/T0,linestyle=2 oplot,x,struct(4,*)/p0,linestyle=5 ; p oplot,x,struct(10,*)/Lsun,linestyle=3 print,' coeff',coeff xyouts,0.35,0.75,' T(r)',charsize=2 xyouts,0.8,0.8,' M(r)',charsize=2 xyouts,0.14,0.6,' L(r)',charsize=2 xyouts,0.3,0.11,'rho(r)',charsize=2 xyouts,0.3,0.07,'P(r)',charsize=2 end