PROGRAM HW5A c homework 5A implicit none integer m,i,j,k,ip,jmax,kmax double precision Pg,Prad,PN,Pe,kB,arad,T,sumj,sumk double precision A,alpha,x,f,ntot,n,nN,ne,nk,njk,nejk,nejkmax double precision mu,muN,mue,rho,rhoN,rhoe,beta,ma,me,RH dimension alpha(2),x(2),A(2),f(3,3),nk(2),njk(3,2) c declarations for problem 2 double precision gi,gip,chii,chiip,ripi,rip1 c declarations for problem 3 double precision U,IE,phij,fj26,P,S,CT,Y dimension U(3),IE(2),fj26(3),P(3),Y(3) c useful constants arad = 7.5676e-15 ! radiation constant kB = 1.3807e-16 ! Boltzman constant [erg/K] ma = 1.6605e-24 ! atomic mass unit [g] me = 9.1095e-28 ! electron mass [g] RH = 13.598 ! Rydberg constant for H A(1) = 1.00794 ! H atomic mass A(2) = 4.00260 ! He atomic mass c --------------------------------------------------------------------------- c PROBLEM 1 c --------------------------------------------------------------------------- c gas properties data PN = 100.0 ! nuclear pressure [dynes/cm^2] T = 10000. ! temperature [K] x(1) = 0.70 ! H mass fraction x(2) = 0.30 ! He mass fraction f(1,1) = 0.15 ! HI ionization fraction f(2,1) = 0.85 ! HII ionization fraction f(1,2) = 0.01 ! HeI ionization fraction f(2,2) = 0.40 ! HeII ionization fraction f(3,2) = 0.59 ! HeIII ionization fraction c (1a) nuclear number density nN = PN/(kB*T) c (1b) abundance fractions sumk = 0.0 DO k=1,2 sumk = sumk + x(k)/A(k) ! sum the denominator ENDDO alpha(1) = x(1)/A(1)/sumk alpha(2) = x(2)/A(2)/sumk c (1c,1d) n_k and n_jk DO k=1,2 nk(k) = alpha(k)*nN DO j=1,k+1 njk(j,k) = f(j,k)*nk(k) ENDDO ENDDO c (1e) electron density, ne ne = 0.0 nejkmax = 0.0 DO k=1,2 DO j=1,k+1 nejk = float(j-1)*njk(j,k) ! compute contribution from ion j,k ne = ne + nejk ! sum over ions nejkmax = max(nejk,nejkmax) ! store the maximum contribution IF (nejk.eq.nejkmax) then ! determine the ion j,k with max nejk jmax = j kmax = k ENDIF ENDDO ENDDO c (1f) total particle density n = nN + ne c (1g) gas pressure (particles) Pg = (nN+ne)*kB*T c (1h) radiation pressure, beta Prad = (arad/3.0)*(T**4) beta = Pg/(Pg+Prad) c (1i) mean molecular weights muN = 0.0 mue = 0.0 DO k=1,2 muN = muN + x(k)/a(k) sumj = 0.0 DO j=1,k+1 sumj = sumj + float(j-1)*f(j,k) ENDDO mue = mue + x(k)/a(k) * sumj ENDDO mu = muN + mue ! sum the inverses mu = 1.0/mu muN = 1.0/muN mue = 1.0/mue c (1j) mass densities rhoN = PN*muN*ma/(kB*T) rhoe = ne*me ! 1st way C rhoe = rhoN/(mue*ma) ! 2nd way rho = rhoN + rhoe sumj = rhoe/rhoN c communicate P1 results WRITE(6,*) ' ' WRITE(6,*) 'PROBLEM 1' WRITE(6,*) ' ' WRITE(6,100) '(a) n_N =',nN WRITE(6,100) '(b) alpha_1 =',alpha(1) WRITE(6,100) ' alpha_2 =',alpha(2) WRITE(6,100) '(c) n_1 =',nk(1) WRITE(6,100) ' n_2 =',nk(2) WRITE(6,100) '(d) n_11 =',njk(1,1) WRITE(6,100) ' n_21 =',njk(2,1) WRITE(6,100) ' n_12 =',njk(1,2) WRITE(6,100) ' n_22 =',njk(2,2) WRITE(6,100) ' n_32 =',njk(3,2) WRITE(6,100) '(e) n_e =',ne WRITE(6,101) ' j,k max =',jmax,kmax WRITE(6,100) ' %max =',100.*nejkmax/ne WRITE(6,100) '(f) n =',n WRITE(6,100) '(g) P_g =',Pg WRITE(6,100) '(h) P_rad =',Prad WRITE(6,100) ' beta =',beta WRITE(6,100) '(i) mu_N =',muN WRITE(6,100) ' mu_e =',mue WRITE(6,100) ' mu =',mu WRITE(6,100) '(j) rho_N =',rhoN WRITE(6,100) ' rho_e =',rhoe WRITE(6,100) ' rho =',rho WRITE(6,100) ' re/rN =',sumj 100 FORMAT(1x,a13,1pe12.3) 101 FORMAT(1x,a13,i4','i1) c --------------------------------------------------------------------------- c PROBLEM 2 c --------------------------------------------------------------------------- WRITE(6,*) ' ' WRITE(6,*) 'PROBLEM 2' WRITE(6,*) ' ' c (2a,2b) Boltzman ratios of different excited states WRITE(6,*) '(a,b) Boltzmann Excitation Ratios T=10,000 K' WRITE(6,201) 201 FORMAT(1x,'i/1',3x,'n_i/n_1',5x,'i+1/i',3x,'n_i+1/n_i') T = 10000. kB = 8.6173e-5 ! Boltzman constant eV/K c loop over excitation states i and i+1 DO i=1,4 ip = i + 1 gi = 2.0*float(i**2) ! stat weight i gip = 2.0*float(ip**2) ! stat weight i+1 chii = RH*(1.0-(1.0/float(i*i))) ! excite energy i chiip = RH*(1.0-(1.0/float(ip*ip))) ! excite energy i+1 ripi = (gip/gi)*exp(-(chiip-chii)/(kB*T)) ! density ratio i+1/i rip1 = (gip/2.0)*exp(-chiip/(kB*T)) ! density ratio i+1/1 (i+1=2-5) WRITE(6,202) ip,1,rip1,ip,i,ripi ! communicate ENDDO 202 FORMAT(1x,i1,'/',i1,1pe12.3,5x,i1,'/',i1,1pe12.3) c (2c) partition functions WRITE(6,*) ' ' WRITE(6,*) '(c) Partition Function T=10,000 K' U(1) = 0.0 c loop over excitation states i DO i=1,5 gi = 2.0*float(i**2) ! stat weight i chii = RH*(1.0-(1.0/float(i*i))) ! excite energy i U(1) = U(1) + gi*exp(-chii/(kB*T)) ! sum U WRITE(6,203) i,gi*exp(-chii/(kB*T)) ! write U for level i ENDDO WRITE(6,204) U(1) 203 FORMAT(1x,' i = ',i1,1x,' U_i11 = ',1pe11.3) 204 FORMAT(1x,' total i=1-5 U_11 = ',1pe11.3) c (2d) Boltzman ratios of excited states to total neutral hydrogens WRITE(6,*) ' ' WRITE(6,*) '(d) Boltzmann Excitation/Total Ratios T=10,000 K' c loop over excitation states i for hydrogenic atom DO i=1,5 gi = 2.0*float(i**2) ! stat weight i chii = RH*(1.0-(1.0/float(i*i))) ! excite energy i rip1 = (gi/U(1))*exp(-chii/(kB*T)) ! ratio WRITE(6,205) i,rip1 ! write U for level i ENDDO 205 FORMAT(1x,' i = ',i1,1x,' n_i11/n11 = ',1pe12.4) c --------------------------------------------------------------------------- c PROBLEM 3 c --------------------------------------------------------------------------- WRITE(6,*) ' ' WRITE(6,*) 'PROBLEM 3' WRITE(6,*) ' ' WRITE(6,*) 'Ionization fractions for FeII' WRITE(6,300) 300 FORMAT(1x,t12,'ne',t23,'T',t37,'f_1,26',t48,'f_2,26',t59,'f_3,26') c properties of Fe ionization IE(1) = 7.9024 ! Fe ionization E j=1 [eV] IE(2) = 16.1877 ! Fe ionization E j=2 [eV] c (3a) G2 V star ne = 1.0e14 T = 5800. U(1) = 10.0**0.736 ! Fe partition function j=1 U(2) = 10.0**0.938 ! Fe partition function j=2 U(3) = 1.0 ! Fe partition function j=3 P(1) = 1.0 ! set up P_1=1 for recursive Saha S = P(1) ! first term of the sum CT = 4.83e15*(T**(1.5)) ! C_phi * T^3/2 DO j=1,2 ! execute recursive Saha phij = CT * (U(j+1)/U(j)) * exp(-IE(j)/(kB*T)) Y(j) = phij/ne ! Saha P(j+1) = P(j)*Y(j) ! recursive Saha S = S + P(j+1) ! sum denominator ENDDO c compute ionization fractions fj26(1) = 1.0/S ! j=1 ionization fraction DO j=2,3 fj26(j) = fj26(j-1)*Y(j-1) ! recursion for j>1 ENDDO WRITE(6,301) 'G2 V ',ne,T,(fj26(j),j=1,3) 301 FORMAT(1x,a5,2x,1p2e11.3,3x,1p3e11.3) c store for % diff calc; use problem 1 arrays for storage f(1,1) = fj26(1) f(2,1) = fj26(2) f(3,1) = fj26(3) c (3b) G2 Ia star (same T, so same U) ne = 1.0e12 P(1) = 1.0 ! set up P_1=1 for recursive Saha S = P(1) ! reset first term of the sum DO j=1,2 ! execute recursive Saha phij = CT * (U(j+1)/U(j)) * exp(-IE(j)/(kB*T)) Y(j) = phij/ne P(j+1) = P(j)*Y(j) S = S + P(j+1) ENDDO fj26(1) = 1.0/S ! j=1 ionization fraction DO j=2,3 fj26(j) = fj26(j-1)*Y(j-1) ! recursion for j>1 ENDDO WRITE(6,301) 'G2 Ia',ne,T,(fj26(j),j=1,3) c store for % diff calc; use problem 1 arrays for storage f(1,2) = fj26(1) f(2,2) = fj26(2) f(3,2) = fj26(3) c (3c) A0 V star (reset ne, T, and U) ne = 1.0e14 T = 10000. U(1) = 10.0**0.770 U(2) = 10.0**0.962 U(3) = 1.0 CT = 4.83e15*(T**(1.5)) P(1) = 1.0 ! set up P_1=1 for recursive Saha S = P(1) ! reset first term of the sum DO j=1,2 ! execute recursive Saha phij = CT * (U(j+1)/U(j)) * exp(-IE(j)/(kB*T)) Y(j) = phij/ne P(j+1) = P(j)*Y(j) S = S + P(j+1) ENDDO fj26(1) = 1.0/S ! j=1 ionization fraction DO j=2,3 fj26(j) = fj26(j-1)*Y(j-1) ! recursion for j>1 ENDDO WRITE(6,301) 'A0 V ',ne,T,(fj26(j),j=1,3) c store for % diff calc f(1,3) = fj26(1) f(2,3) = fj26(2) f(3,3) = fj26(3) c (3d.i) compute percent differences WRITE(6,*) ' ' WRITE(6,*) 'Ratio from G2 V star' WRITE(6,310) 310 FORMAT(1x,t12,'ne',t23,'T',t37,'r_1,26',t48,'r_2,26',t59,'r_3,26') WRITE(6,301) 'G2 V ',1.e14,5800.,(f(j,1)/f(j,1),j=1,3) WRITE(6,301) 'G2 Ia',1.e12,5800.,(f(j,2)/f(j,1),j=1,3) WRITE(6,301) 'A0 V ',1.e14,10000.,(f(j,3)/f(j,1),j=1,3) STOP END