PROGRAM BALMER C ***** C THIS COMPUTES THE FIRST FOUR EMERGENT PROFILES OF THE BALMER SERIES C OF HYDROGEN C SET UP FOR ATLAS9. READS FILE VCSBALMER ON 25 C C THE VCS THEORY IS USED C C IN AN APPROXIMATE WAY WE HAVE ACCOUNTED FOR RADIATIVE AND RESONANCE C IN THAT LORENTZ PROFILES ARE ADDED IN OUTSIDE THE CORE C ***** IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (kw=99) COMMON /ABROSS/ABROSS(kw),TAUROS(kw) COMMON /ABTOT/ABTOT(kw),ALPHA(kw) COMMON /CONT/ABTOTC(kw),ALPHAC(kw),TAUNUC(kw),SNUC(kw),HNUC(kw), 1 JNUC(kw),JMINSC(kw),RESIDC(kw) REAL*8 JNUC,JMINSC COMMON /CONV/DLTDLP(kw),HEATCP(kw),DLRDLT(kw),VELSND(kw), 1 GRDADB(kw),HSCALE(kw),FLXCNV(kw),VCONV(kw),MIXLTH, 2 OVERWT,FLXCNV0(kw),FLXCNV1(kw),IFCONV REAL*8 MIXLTH COMMON /DEPART/BHYD(kw,6),BMIN(kw),NLTEON COMMON /ELEM/ABUND(99),ATMASS(99),ELEM(99) COMMON /EDENS/EDENS(kw),IFEDNS COMMON /FILENAME/FILENAME CHARACTER*60 FILENAME COMMON /FLUX/ FLUX,FLXERR(kw),FLXDRV(kw),FLXRAD(kw) COMMON /FREQ/FREQ,FREQLG,EHVKT(kw),STIM(kw),BNU(kw) COMMON /FRESET/FRESET(1563),RCOSET(1563),NULO,NUHI,NUMNU COMMON /HEIGHT/HEIGHT(kw) COMMON /IF/IFCORR,IFPRES,IFSURF,IFSCAT,TAUSCAT,IFMOL COMMON /IFOP/IFOP(20) COMMON /IONS/XNFPH(kw,2),XNFPHE(kw,3),XNFH(kw,2),XNFHE(kw,3) COMMON /ITER/ ITER,IFPRNT(15),IFPNCH(15),NUMITS COMMON /JUNK/TITLE(74),FREQID(6),WLTE,XSCALE COMMON /MUS/ANGLE(20),SURFI(20),NMU COMMON /OPS/AHYD(kw),AH2P(kw),AHMIN(kw),SIGH(kw),AHE1(kw), 1 AHE2(kw),AHEMIN(kw),SIGHE(kw),ACOOL(kw),ALUKE(kw),AHOT(kw), 2 SIGEL(kw),SIGH2(kw),AHLINE(kw),ALINES(kw),SIGLIN(kw), 3 AXLINE(kw),SIGXL(kw),AXCONT(kw),SIGX(kw),SHYD(kw), 4 SHMIN(kw),SSSSSS(kw),SXLINE(kw),SXCONT(kw) DIMENSION SHLINE(kw) COMMON /OPTOT/ACONT(kw),SCONT(kw),ALINE(kw),SLINE(kw),SIGMAC(kw), 1 SIGMAL(kw) COMMON /PTOTAL/PTOTAL(kw) COMMON /PUT/PUT,IPUT COMMON /PZERO/PZERO,PCON,PRADK0,PTURB0,KNU(kw),PRADK(kw),RADEN(kw) REAL*8 KNU COMMON /RAD/ ACCRAD(kw),PRAD(kw) COMMON /RHOX/RHOX(kw),NRHOX COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw) COMMON /TAUSHJ/TAUNU(kw),SNU(kw),HNU(kw),JNU(kw),JMINS(kw) REAL*8 JNU,JMINS COMMON /STEPLG/STEPLG,TAU1LG,KRHOX COMMON /TEFF/TEFF,GRAV,GLOG COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),ITEMP COMMON /TURBPR/VTURB(kw),PTURB(kw),TRBFDG,TRBCON,TRBPOW,TRBSND, 1 IFTURB COMMON /WAVEY/WBEGIN,DELTAW,IFWAVE COMMON /XABUND/XABUND(99),WTMOLE COMMON /XNF/XNFC(kw,6),XNFN(kw,6),XNFO(kw,6),XNFNE(kw,6), 1 XNFMG(kw,6),XNFSI(kw,6),XNFS(kw,6),XNFFE(kw,5) COMMON /XNFP/XNFPC(kw,4),XNFPN(kw,5),XNFPO(kw,6),XNFPNE(kw,6), 1 XNFPAL(kw,1),XNFPMG(kw,2),XNFPSI(kw,2),XNFPCA(kw,2), 2 XNFPFE(kw,1),XNFPCH(kw),XNFPOH(kw) DIMENSION PART(kw,6) COMMON /IFPOP/IFPOP COMMON /IFEQUA/IFEQUA(101),KCOMPS(450),LOCJ(161),CODE(160), 1 EQUIL(6,160),IDEQUA(25),NEQUA,NUMMOL,NLOC COMMON /XNSAVE/XNSAVE(kw,25) CHARACTER*60 BFILENAME DIMENSION EQWID(5),FMINUS(5),FPLUS(5) DIMENSION DELWT(95),DEL(95),PR1(95),PR2(95),PR(95,kw) DIMENSION SPECTRUM(99),RESID(99,5) DIMENSION FPM(5),FN(5),RADAMP(5),ALAM(5) DIMENSION XN2(kw) DATA NDEL/95/ DATA DEL/ 1-100.,-90.,-80.,-70.,-60.,-55.,-50.,-45.,-40.,-36.,-32.,-30.,-28., 2 -26.,-24.,-22.,-20.,-18.,-16.,-14.,-12.,-10., -9., -8., -7., -6., 3 -5.0,-4.5,-4.0,-3.5,-3.0,-2.5,-2.0,-1.8,-1.6,-1.4,-1.2,-1.0, -.9, 4 -.8, -.7, -.6, -.5, -.4, -.3, -.2, -.1, 0.0, .1, .2, .3, .4, 5 .5, .6, .7, .8, .9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 6 3.5, 4.0, 4.5, 5.0, 6., 7., 8., 9., 10., 12., 14., 16., 18., 7 20., 22., 24., 26., 28., 30., 32., 36., 40., 45., 50., 55., 60., 8 70., 80., 90.,100./ C INTEGRATION WEIGHTS ARE IN UNITS OF 1/3 DATA DELWT/ 1 10., 40., 20., 40., 15., 20., 10., 20., 9., 16., 7., 5., 8., 2 4., 8., 4., 8., 4., 8., 4., 8., 3.5, 2.5, 4., 2., 4., 3 1.5, 2., 1., 2., 1., 2., .8, .5, .8, .4, .8, .3, .4, 4 .2, .4, .2, .4, .2, .4, .2, .4, .2, .4, .2, .4, .2, 5 .4, .2, .4, .2, .4, .3, .8, .4, .8, .5, .8, 2., 1., 6 2., 1., 2., 1.5, 4., 2., 4., 2.5, 3.5, 8., 4., 8., 4., 7 8., 4., 8., 4., 8., 5., 7., 16., 9., 20., 10., 20., 15., 8 40., 20., 40., 10./ DATA FPM/.39206,.11932,.037656,.022093,.011394/ DATA FN/.6408,.11932,.04467,.022093,.01271/ DATA RADAMP/0.000626,.000313,.000239,.000211,.000236/ ITEMP=0 ITER=1 1 CALL READIN(2) IFOP(11)=0 IFOP(14)=0 IFOP(15)=0 IFPRES=0 5 ITEMP=ITEMP+ITER CALL POPS(0.D0,1,XNE) CALL POPS(1.01D0,11,XNFPH) CALL POPS(2.02D0,11,XNFPHE) CALL POPS(6.03D0,11,XNFPC) CALL POPS(7.04D0,11,XNFPN) CALL POPS(8.05D0,11,XNFPO) CALL POPS(10.05D0,11,XNFPNE) CALL POPS(12.01D0,11,XNFPMG) CALL POPS(13.00D0,11,XNFPAL) CALL POPS(14.01D0,11,XNFPSI) CALL POPS(20.01D0,11,XNFPCA) CALL POPS(26.00D0,11,XNFPFE) IF(IFMOL.EQ.0)THEN CALL POPS(1.01D0,12,XNFH) CALL POPS(2.02D0,12,XNFHE) CALL POPS(6.05D0,12,XNFC) CALL POPS(7.05D0,12,XNFN) CALL POPS(8.05D0,12,XNFO) CALL POPS(10.05D0,12,XNFNE) CALL POPS(12.05D0,12,XNFMG) CALL POPS(14.05D0,12,XNFSI) CALL POPS(16.05D0,12,XNFS) CALL POPS(26.04D0,12,XNFFE) ENDIF IF(IFMOL.EQ.1)THEN CALL POPS(106.00D0,11,XNFPCH) CALL POPS(108.00D0,11,XNFPOH) C CALL W(6HXNFPCH,XNFPCH,NRHOX) C CALL W(6HXNFPOH,XNFPOH,NRHOX) C THE POPS WILL NOT RETURN NUMBER DENSITIES WHEN MOLECULES ARE ON C SO WE COMPUTE NUMBER DENSITIES/PART FUNCTIONS AND PART FUNCTIONS CALL POPS(6.05D0,11,XNFC) CALL POPS(7.05D0,11,XNFN) CALL POPS(8.05D0,11,XNFO) CALL POPS(10.05D0,11,XNFNE) CALL POPS(12.05D0,11,XNFMG) CALL POPS(14.05D0,11,XNFSI) CALL POPS(16.05D0,11,XNFS) CALL POPS(26.04D0,11,XNFFE) DO 444 J=1,NRHOX CALL PFSAHA(J,1,1,3,PART) XNFH(J,1)=XNFPH(J,1)*PART(J,1) XNFH(J,2)=XNFPH(J,2) CALL PFSAHA(J,2,2,13,PART) XNFHE(J,1)=XNFPHE(J,1)*PART(J,1) XNFHE(J,2)=XNFPHE(J,2)*PART(J,2) XNFHE(J,3)=XNFPHE(J,3) CALL PFSAHA(J,6,6,13,PART) XNFC(J,1)=XNFC(J,1)*PART(J,1) XNFC(J,2)=XNFC(J,2)*PART(J,2) XNFC(J,3)=XNFC(J,3)*PART(J,3) XNFC(J,4)=XNFC(J,4)*PART(J,4) XNFC(J,5)=XNFC(J,5)*PART(J,5) XNFC(J,6)=XNFC(J,6)*PART(J,6) CALL PFSAHA(J,7,6,13,PART) XNFN(J,1)=XNFN(J,1)*PART(J,1) XNFN(J,2)=XNFN(J,2)*PART(J,2) XNFN(J,3)=XNFN(J,3)*PART(J,3) XNFN(J,4)=XNFN(J,4)*PART(J,4) XNFN(J,5)=XNFN(J,5)*PART(J,5) XNFN(J,6)=XNFN(J,6)*PART(J,6) CALL PFSAHA(J,8,6,13,PART) XNFO(J,1)=XNFO(J,1)*PART(J,1) XNFO(J,2)=XNFO(J,2)*PART(J,2) XNFO(J,3)=XNFO(J,3)*PART(J,3) XNFO(J,4)=XNFO(J,4)*PART(J,4) XNFO(J,5)=XNFO(J,5)*PART(J,5) XNFO(J,6)=XNFO(J,6)*PART(J,6) CALL PFSAHA(J,10,6,13,PART) XNFNE(J,1)=XNFNE(J,1)*PART(J,1) XNFNE(J,2)=XNFNE(J,2)*PART(J,2) XNFNE(J,3)=XNFNE(J,3)*PART(J,3) XNFNE(J,4)=XNFNE(J,4)*PART(J,4) XNFNE(J,5)=XNFNE(J,5)*PART(J,5) XNFNE(J,6)=XNFNE(J,6)*PART(J,6) CALL PFSAHA(J,12,6,13,PART) XNFMG(J,1)=XNFMG(J,1)*PART(J,1) XNFMG(J,2)=XNFMG(J,2)*PART(J,2) XNFMG(J,3)=XNFMG(J,3)*PART(J,3) XNFMG(J,4)=XNFMG(J,4)*PART(J,4) XNFMG(J,5)=XNFMG(J,5)*PART(J,5) XNFMG(J,6)=XNFMG(J,6)*PART(J,6) CALL PFSAHA(J,14,6,13,PART) XNFSI(J,1)=XNFSI(J,1)*PART(J,1) XNFSI(J,2)=XNFSI(J,2)*PART(J,2) XNFSI(J,3)=XNFSI(J,3)*PART(J,3) XNFSI(J,4)=XNFSI(J,4)*PART(J,4) XNFSI(J,5)=XNFSI(J,5)*PART(J,5) XNFSI(J,6)=XNFSI(J,6)*PART(J,6) CALL PFSAHA(J,16,6,13,PART) XNFS(J,1)=XNFS(J,1)*PART(J,1) XNFS(J,2)=XNFS(J,2)*PART(J,2) XNFS(J,3)=XNFS(J,3)*PART(J,3) XNFS(J,4)=XNFS(J,4)*PART(J,4) XNFS(J,5)=XNFS(J,5)*PART(J,5) XNFS(J,6)=XNFS(J,6)*PART(J,6) CALL PFSAHA(J,26,5,13,PART) XNFFE(J,1)=XNFFE(J,1)*PART(J,1) XNFFE(J,2)=XNFFE(J,2)*PART(J,2) XNFFE(J,3)=XNFFE(J,3)*PART(J,3) XNFFE(J,4)=XNFFE(J,4)*PART(J,4) XNFFE(J,5)=XNFFE(J,5)*PART(J,5) 444 CONTINUE ENDIF C PRINT 7777,XNE C PRINT 7777,XNFH C PRINT 7777,XNFPH C PRINT 7777,XNFPHE C PRINT 7777,XNFHE 7777 FORMAT(1P10E12.4) C IFSURF=1 DO 30 J=1,NRHOX 30 XN2(J)=EXP(-11.8353E4/T(J))*8.*XNFPH(J,1)/RHO(J) C PRINT 7777,XN2 DO 230 LINE=1,4 M=LINE+2 TM=FLOAT(M*M) GNM=1./4.-1./TM VNM=GNM*3.28805E15 XLAM=2.99792458E18/VNM ALAM(LINE)=XLAM DO 90 J=1,NRHOX C RESON=XLAM**2*1.71E-27*XNATOM(J)*XNF(J) RESON=XLAM**2*1.71E-27*XNFPH(J,1)*2. A=RADAMP(LINE)+RESON DO 60 I=1,NDEL 60 PR2(I)=A/3.1416/(DEL(I)**2+A*A) CALL VCS(PR1,XNE(J),T(J),DEL,NDEL,2,M) E=EXP(-4.7993E-11*VNM/T(J)) STIML=BHYD(J,2)-E*BHYD(J,M) SHLINE(J)=.524946*GNM**3/(BHYD(J,2)/BHYD(J,M)/E-1.) DO 90 I=1,NDEL PR(I,J)=PR1(I)*FN(LINE) IF(ABS(DEL(I)).GT..2) PR(I,J)=PR(I,J)+PR2(I)*FN(LINE) PR(I,J)=PR(I,J)*STIML*XN2(J)/VNM**2*7.957E16 90 CONTINUE DO 99 I=1,NDEL WAVE=XLAM+DEL(I) FREQ=2.99792458E18/WAVE FREQLG=LOG(FREQ) FREQ15=FREQ/1.E15 DO 95 J=1,NRHOX EHVKT(J)=EXP(-FREQ*HKT(J)) STIM(J)=1.-EHVKT(J) 95 BNU(J)=1.47439E-02*FREQ15**3*EHVKT(J)/STIM(J) CALL KAPP(0,NSTEPS,STEPWT) DO 96 J=1,NRHOX ALINE(J)=PR(I,J) 96 SLINE(J)=SHLINE(J) CALL JOSH(IFSCAT,IFSURF) 99 SPECTRUM(I)=HNU(1) FMINUS(LINE)=4.*SPECTRUM(1) FPLUS(LINE)=4.*SPECTRUM(NDEL) EQWID(LINE)=0. DO 190 I=1,NDEL CONT=(DEL(I)-DEL(1))/(DEL(NDEL)-DEL(1))* 1(SPECTRUM(NDEL)-SPECTRUM(1))+SPECTRUM(1) RESID(I,LINE)=SPECTRUM(I)/CONT EQWID(LINE)=EQWID(LINE)+DELWT(I)*(1.-RESID(I,LINE)) 190 CONTINUE C INTEGRATION WEIGHTS ARE IN UNITS OF 1/3 EQWID(LINE)=EQWID(LINE)/3. 230 CONTINUE PRINT 245,TEFF,GLOG,TITLE,FILENAME 245 FORMAT ('1TEFF',F9.1,2X,'LOG G',F9.5,6X/1X,74A1/1X,A60// 1 14X,'DEL',10X,'HALPHA',11X,'HBETA',12X,'HGAMMA',11X,'HDELTA') BFILENAME='B'//FILENAME(2:)//'.DAT' OPEN(UNIT=7,NAME=BFILENAME,STATUS='NEW') WRITE(7,250)TEFF,GLOG,FILENAME,TITLE 250 FORMAT(4HTEFF,F9.1,6H LOG G,F9.5,3X,A60/74A1) WRITE(7,302)(EQWID(LINE),LINE=1,4) 302 FORMAT('EQ WIDTH(A) HALPHA,HBETA,HGAMMA,HDELTA ',4F8.3) WRITE(7,303)(FMINUS(LINE),LINE=1,4) 303 FORMAT('FNU(-100A) ',1P4E11.3) WRITE(7,304)(FPLUS(LINE),LINE=1,4) 304 FORMAT('FNU(+100A) ',1P4E11.3) DO 290 I=48,95 PRINT 270,DEL(96-I),DEL(I),(RESID(I,LINE),RESID(96-I,LINE), 1 LINE=1,4) WRITE(7,280)DEL(96-I),DEL(I),(RESID(I,LINE),RESID(96-I,LINE), 1 LINE=1,4) 270 FORMAT(8X,2F6.1,4(F9.3,F7.3,1X)) 280 FORMAT('DEL-,DEL+(A),R-,R+ ',2F6.1,8F6.3) 290 CONTINUE PRINT 300, (EQWID(LINE),LINE=1,4) 300 FORMAT ('0EQ WIDTH(A)',4X,4(10X,F7.3)) PRINT 301,(FMINUS(LINE),LINE=1,4) 301 FORMAT('0FNU(-100A)',12X,1P4(E12.3,5X)) PRINT 311,(FPLUS(LINE),LINE=1,4) 311 FORMAT(' FNU(+100A)',12X,1P4(E12.3,5X)) CLOSE(UNIT=7) GO TO 1 END SUBROUTINE VCS (PR,XNE,T,DEL,II,N,M) C CALCULATES VIDAL COOPER AND SMITH PROFILES FOR FIRST FOUR BALMER C LINES. THEY ARE RETURNED IN ANGSTROM UNITS. C ASSUMES RAW PROFILES ARE IN ALPHA UNITS. IMPLICIT REAL*8 (A-H,O-Z) DIMENSION PR(40),DEL(40) C DIMENSION ALPHA(40),PRALPH(40) DIMENSION PRALPH(40) DIMENSION SVCS(6,17,40,4),ALPHA0(4) DATA SVCS(1,1,1,1)/0./,ALPHA0/-3.,-3.,-3.,-3./ EXP10(X)=EXP(2.30258509299405E0*X) IF(SVCS(1,1,1,1).NE.0.)GO TO 3 C READ IN VCS ARRAYS READ(25,10) DO 20 LINE=1,4 READ(25,10) READ(25,10)(((SVCS(I,J,K,LINE),J=1,17),I=1,6),K=1,40) 10 FORMAT (10F8.4) 20 CONTINUE 3 LINE=M-N C TEMPERATURE AND ELECTRON DENSITY INTERPOLATION AT=LOG10(T) BTEMP=(AT-3.3979411)/.3010300+1. ITEMP=BTEMP C ITEMP=(LOG10(T)-LOG10(2500.))/LOG10(2.)+1. ITEMP=MAX(MIN(ITEMP,5),1) WTTEMP=BTEMP-ITEMP ANE=LOG10(XNE) C TO AVOID BAD EXTRAPOLATION AT LOW DENSITY ANE=MAX(10.D0,ANE) BNE=(ANE-10.)/.5+1. INE=BNE INE=MAX(MIN(INE,16),1) WTXNE=BNE-INE DO 21 I=1,40 21 PRALPH(I)=(1.-WTXNE)*(1.-WTTEMP)*SVCS(ITEMP ,INE ,I,LINE)+ 1 (1.-WTXNE)*WTTEMP*SVCS(ITEMP+1,INE ,I,LINE)+ 2 WTXNE*(1.-WTTEMP)*SVCS(ITEMP ,INE+1,I,LINE)+ 3 WTXNE*WTTEMP*SVCS(ITEMP+1,INE+1,I,LINE) C NOW ALPHA INTERPOLATION FO=1.25E-9*XNE**.66666667 DO 50 I=1,II IF (DEL(I).NE.0.) GO TO 25 PR(I)=EXP10(PRALPH(1))/FO GO TO 50 25 AALP=LOG10(ABS(DEL(I))/FO) BALP=(AALP-ALPHA0(LINE))/.2+1. IALP=BALP IF (IALP.GT.1) GO TO 30 PR(I)=EXP10(PRALPH(1))/FO GO TO 50 30 IF (IALP.GT.39) GO TO 40 WT=BALP-IALP PR(I)=EXP10((1.-WT)*PRALPH(IALP)+WT*PRALPH(IALP+1))/FO GO TO 50 40 PR(I)=EXP10(PRALPH(40)+2.5*(ALPHA0(LINE)+7.8-AALP))/FO 50 CONTINUE RETURN END SUBROUTINE ATLAS9 RETURN END SUBROUTINE LINOP RETURN END