program simpstar * computes the structures of 2 very simple stellar models in hydrostatic * equilibrium. the homogeneous model assumes constant mass density through * the star; the linear model has a linear change in density with r/R through * the star. perfect gas law and radiative transport assumed; very simple * recipe for the opacity; no central energy source; chemical abundances * constant throughout. the runs in P(r) and T(r) were prior determined with * paper and pencil; set temperature gradient equal to the radiative one at * r/R = 0.5 to estimate the star's luminosity. * ! last edited by ktk: March 2007 * implicit none integer i, nmax real xincr parameter (nmax = 1001) double precision G, mH, k, pi, c, a, Msun, Rsun, Lsun, Tsun double precision sigma, Rad, R, Mass, M double precision Xm, Ym, Zm, mu, muS, rhobarS, rhobar double precision Xsun, Ysun, Zsun double precision x(nmax), x2, x3, x4, gff * specific to model (homogeneous, linear) double precision rhoc1, Pc1, Tc1, rhor1(nmax), Mr1(nmax) double precision Pr1(nmax), Tr1(nmax), Tmean1 double precision gbft1(nmax), gbfthalf1 double precision rhoc2, Pc2, Tc2, rhor2(nmax), Mr2(nmax) double precision Pr2(nmax), Tr2(nmax), Tmean2 double precision gbft2(nmax), gbfthalf2 * these variables pertain to the radiative temperature gradient double precision Lum1, kapparho1(nmax), mfp1(nmax), kappa1(nmax) double precision Lum2, kapparho2(nmax), mfp2(nmax), kappa2(nmax) double precision rhalf1,rhohalf1,Thalf1,dTdrhalf1 double precision kapparhohalf1 double precision rhalf2,rhohalf2,Thalf2,dTdrhalf2 double precision kapparhohalf2 * some constants (cgs) CODATA 2006 G = 6.67428d-08 mH = 1.673532499d-24 k = 1.3806504d-16 c = 2.99792458d+10 sigma = 5.670400d-05 pi = 3.1415926535897932d+00 a = 4.d0*sigma/c * these have to do with the Sun (cgs) Msun = 1.9891d+33 Rsun = 6.9566d+10 ! most recent measurements Lsun = 3.841d+33 ! most recent measurements Tsun = (Lsun/(4.d0*pi*Rsun*Rsun*sigma))**0.25 ! not presently used ! the Sun's probable mean ZAMS composition (bulk composition at birth) ! Bahcall et al. (2001) Standard Solar Model Xsun = 0.7077d0 Ysun = 0.2735d0 Zsun = 1.d0 - (Xsun + Ysun) * input the key variables print* write(6,*)'enter M and R of the star relative to solar values' read(5,*) Mass, Rad print* write(6,*)'enter H and He mass fractions X and Y:' read(5,*) Xm, Ym print * Zm = 1.d0 - (Xm + Ym) ! the composition is assumed constant throughout the star * mean molecular weight of the fully ionized gas mixture, assumed constant mu = 1.d0/(2.d0*Xm + 0.75d0*Ym + 0.56d0*Zm) muS = 1.d0/(2.d0*Xsun + 0.75d0*Ysun + 0.56d0*Zsun) ! solar abundances * convert input mass of star to grams M = Mass * Msun * convert input radius of star to cm R = Rad * Rsun * the mean solar mass density gm/cm^3 rhobarS = 3.d0 * Msun / (4.d0*pi*Rsun*Rsun*Rsun) * the mean stellar mass density rhobar = rhobarS * Mass / (Rad * Rad * Rad) * * x is a variable measuring r/R, varying between 0 and 1. * zero out its array... do i = 1, nmax x(i) = 0. end do * the increment xincr = 1./float(nmax - 1) * write(6,*) xincr * now load it up x(1) = 1.d-6 do i = 2, nmax x(i) = x(i-1) + xincr end do * back off just a bit from the end point r = R x(nmax) = x(nmax) - 5.d-4 * open up the output files open(10,file='homogen.out') open(11,file='linear.out') * ******************homogeneous model****************** * * the stellar central variables of state rhoc1 = 1.d0 * rhobar Pc1 = 3.d0/(8.d0*pi) * G * M * M / (R * R * R * R) Tc1 = 0.5d0 * (G*mu*mH/k) * M / R * some variables evaluated at r = R/2 rhalf1 = 0.5d0 * R rhohalf1 = rhobar Thalf1 = 0.75d0 * Tc1 dTdrhalf1 = Tc1 / R * the loop on M(r), P(r), T(r) and kappa*rho do i = 1, nmax x2 = x(i) * x(i) x3 = x(i) * x(i) * x(i) rhor1(i) = rhoc1 ! constant density model Mr1(i) = M * (x3) Pr1(i) = Pc1 * (1. - x2) if(Pr1(i) .lt. 1.d-30) then Pr1(i) = 1.d-30 endif Tr1(i) = Tc1 * (1. - x2) if(Tr1(i) .lt. 1.d-30) then Tr1(i) = 1.d-30 endif * mean temperature from Virial Theorem, assuming ideal gas pressure Tmean1 = (1.d0/5.d0) * (G*mu*mH/k) * M / R * the opacity ! approximate gaunt, guillotine factors for b-f and f-f opacities gbft1(i) = 1.d0/(2.82d0*(rhor1(i) * (1.d0+Xm))**0.2) gbfthalf1 = 1.d0/(2.82d0*(rhohalf1 * (1.d0+Xm))**0.2) gff = 1.d0 kappa1(i) = (1.d0+Xm) * rhor1(i) c * (Tr1(i)/1.d6)**(-3.5) * ( 4.34d4*(Zm)*gbft1(i) c + 36.8d0*(1.d0-Zm)*gff ) + 0.2004d0*(1.d0+Xm) * the opacity * density or roughly the inverse mean free path kapparho1(i) = rhor1(i) * kappa1(i) * a sort of mean-free-path mfp1(i) = 1.d0/kapparho1(i) end do ! loop on x = r/R * kappa*rho evaluated at r = R/2 kapparhohalf1 = ( (1.d0+Xm) * rhohalf1 c * (Thalf1/1.d6)**(-3.5) * ( 4.34d4*(Zm)*gbfthalf1 c + 36.8d0*(1.d0-Zm)*gff ) + 0.2004d0*(1.d0+Xm) ) c * rhohalf1 * estimate the luminosity from the radiative temperature gradient * evaluated at r = R/2 (x = 0.5). Lum1 = (16.d0*pi*a*c/3.d0) * (rhalf1*rhalf1) * c (Thalf1*Thalf1*Thalf1)*dTdrhalf1 / kapparhohalf1 Lum1 = Lum1/Lsun * ************************************************ * ******************linear model****************** * * the stellar central variables of state rhoc2 = 4.d0 * rhobar Pc2 = 5.d0/(4.d0*pi) * G * M * M / (R * R * R * R) Tc2 = (5.d0/12.d0) * (G*mu*mH/k) * M / R * some variables evaluated at r = R/2 rhalf2 = 0.5d0 * R rhohalf2 = 2.d0 * rhobar Thalf2 = (31.d0/40.d0) * Tc2 dTdrhalf2 = (29.d0/20.d0) * Tc2 / R * the loop on M(r), P(r), T(r) and kappa*rho do i = 1, nmax x2 = x(i) * x(i) x3 = x(i) * x(i) * x(i) x4 = x(i) * x(i) * x(i) * x(i) rhor2(i) = rhoc2 * (1.d0 - x(i)) if(rhor2(i) .lt. 1.d-30) then rhor2(i) = 1.d-30 endif Mr2(i) = M * (4.d0*x3 - 3.d0*x4) Pr2(i) = Pc2 * ( 1.d0 - 24.d0*x2/5.d0 + c 28.d0*x3/5.d0 - 9.d0*x4/5.d0 ) if(Pr2(i) .lt. 1.d-30) then Pr2(i) = 1.d-30 endif Tr2(i) = Tc2 * ( 1.d0 + x(i) - 19.d0*x2/5.d0 + 9.d0*x3/5.d0 ) if(Tr2(i) .lt. 1.d-30) then Tr2(i) = 1.d-30 endif * mean temperature from Virial Theorem, assuming ideal gas pressure Tmean2 = (2.6d1/1.05d2) * (G*mu*mH/k) * M / R * the opacity ! approximate gaunt, guillotine factors for b-f and f-f opacities gbft2(i) = 1.d0/(2.82d0*(rhor2(i) * (1.d0+Xm))**0.2) gbfthalf2 = 1.d0/(2.82d0*(rhohalf2 * (1.d0+Xm))**0.2) gff = 1.d0 kappa2(i) = (1.d0+Xm) * rhor2(i) c * (Tr2(i)/1.d6)**(-3.5) * ( 4.34d4*(Zm)*gbft2(i) c + 36.8d0*(1.d0-Zm)*gff ) + 0.2004d0*(1.d0+Xm) * the opacity * density or roughly the inverse mean free path kapparho2(i) = rhor2(i) * kappa2(i) * a sort of mean-free-path mfp2(i) = 1.d0/kapparho2(i) end do ! loop on x = r/R * kappa*rho evaluated at r = R/2 kapparhohalf2 = ( (1.d0+Xm) * rhohalf2 c * (Thalf2/1.d6)**(-3.5) * ( 4.34d4*(Zm)*gbfthalf2 c + 36.8d0*(1.d0-Zm)*gff ) + 0.2004d0*(1.d0+Xm) ) c * rhohalf2 * estimate the luminosity from the radiative temperature gradient * evaluated at r = R/2 (x = 0.5). Lum2 = (16.d0*pi*a*c/3.d0) * (rhalf2*rhalf2) * c (Thalf2*Thalf2*Thalf2)*dTdrhalf2 / kapparhohalf2 Lum2 = Lum2/Lsun ************************************************ * punch out the data write(10,100) 100 format('# The homogeneous model (cgs units)...') write(10,101) Mass write(10,102) Rad write(10,103) Xm, Ym write(10,203) mu write(10,104) rhobar write(10,105) rhoc1 write(10,106) Pc1 write(10,107) Tc1 write(10,207) Tmean1 write(10,108) Lum1 101 format('# The mass in solar masses: ', 1p,e10.4) 102 format('# The radius in solar radii: ', 1p,e10.4) 103 format('# The H, He mass fractions are: ', 2(1p,e10.4,2x)) 203 format('# The gas mean molecular weight: ', f6.4) 104 format('# The mean density (g/cm^3): ', 1p,e10.4) 105 format('# The central density (g/cm^3): ', 1p,e10.4) 106 format('# The central pressure (dynes/cm^2): ', 1p,e10.4) 107 format('# The central temperature (K): ', 1p,e10.4) 207 format('# The mean temperature (K): ', 1p,e10.4) 108 format('# The estimated L/Lsun is: ', 1p,e10.4) write(10,109) 109 format('#') write(10,110) 110 format('# r/R density M(r)/M Pres. Temp. c kappa mfp') do i = 1, nmax write(10,111) x(i), rhor1(i), Mr1(i)/M, Pr1(i), Tr1(i), c kappa1(i), mfp1(i) 111 format(1p,e8.2,2x,6(1p,e10.4,2x)) end do close(10) * write(11,120) 120 format('# The linear model (cgs units)...') write(11,121) Mass write(11,122) Rad write(11,123) Xm, Ym write(11,223) mu write(11,124) rhobar write(11,125) rhoc2 write(11,126) Pc2 write(11,127) Tc2 write(11,227) Tmean2 write(11,128) Lum2 121 format('# The mass in solar masses: ', 1p,e10.4) 122 format('# The radius in solar radii: ', 1p,e10.4) 123 format('# The H, He mass fractions are: ', 2(1p,e10.4,2x)) 223 format('# The gas mean molecular weight: ', f6.4) 124 format('# The mean density (g/cm^3): ', 1p,e10.4) 125 format('# The central density (g/cm^3): ', 1p,e10.4) 126 format('# The central pressure (dynes/cm^2): ', 1p,e10.4) 127 format('# The central temperature (K): ', 1p,e10.4) 227 format('# The mean temperature (K): ', 1p,e10.4) 128 format('# The estimated L/Lsun is: ', 1p,e10.4) write(11,129) 129 format('#') write(11,130) 130 format('# r/R density M(r)/M Pres. Temp. c kappa mfp') do i = 1, nmax write(11,131) x(i), rhor2(i), Mr2(i)/M, Pr2(i), Tr2(i), c kappa2(i), mfp2(i) 131 format(1p,e8.2,2x,6(1p,e10.4,2x)) end do close(11) write(6,*)'data sent to homogen.out and linear.out' print * stop end