C C CODE FOR CALCULATION OF LTE POPULATION OF DIFFERENT IONIZATION C STAGES OR BOTH LOWER AND UPPER LEVELS OF A SPECTRAL LINE C THROUGHOUT THE STELLAR ATMOSPHERES c c authors: J. Budaj, M.M. Dworetsky c contact: budaj@ta3.sk c published: Budaj J., Dworetsky M.M., 2002, MNRAS, accepted c special requirements: INCLUDE 'pfdwor.inc' c non original routins: 'pfdwor.inc' from: c Budaj J., Dworetsky M.M., Smalley B, 2002, c Comm. Univ. London Obs. No. 82, c http://www.ulo.ucl.ac.uk/ulo\_comms/82/index.html c C C PROGRAM ION C C INPUT: -(3) MODEL OF THE STELLAR ATMOPHERE C (NEW KURUCZ OR HUBENY) C -(7) ATOMIC DATA: FIRST 8 IONIZ.POTENTIALS, C DATA FOR PARTITION FUNCTION CALCULATION C OF FIRST 9 IONIZATION STAGES OF THE ELEMENT OF INTEREST C (AFTER IRWIN,1981,APJS 45,621) C -(8) ATOMIC LEVEL DATA FROM KURUCZ LINELIST FORMAT C OUTPUT: -(2)ION.OUT C C DIMENZATION PARAMETER=NDIM (NUMBER OF THE DEPT POINTS) PARAMETER (NDIM=90,NION1=9) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*30 INMOD,INATOM,INLEV,INN DIMENSION ZIP(NION1-1),STA(6,NION1),STAVS(NDIM,NION1) P ,ZI(NDIM,3),SA(NDIM,2),F1T(NDIM,3) P ,DM(NDIM),HLB(NDIM),TEMP(NDIM),PRESS(NDIM),EKONC(NDIM) P ,AKONC(NDIM),DENS(NDIM) P ,RION(NION1),R1(NION1),RR(NDIM,NION1) C SAHA EQUATION (RATIO OF NUMBER DENSITIES OF TWO ADJACENT C IONS) C ZP -IN [EV],BOLT.CONST.-IN [ERG/K] SAHA(TE,EK,ZP,S1,S2)=1.2020D9*S2/(S1*(5039.77D0/TE)**2.5D0 P *DEXP(2.302585D0*5039.77D0/TE*ZP)*EK*TE*1.38062D-16) C----------------------------------------------------------------------- WRITE(6,'(A,$)')'NEW KURUCZ [0] OR HUBENY MODEL [1]? => ' READ(5,*)IMODEL WRITE(6,'(A,$)')'INPUT FILE WITH THE ATMOSPHERE MODEL => ' READ(5,5)INMOD 5 FORMAT(A30) WRITE(6,'(A,$)')'INPUT FILE WITH ATOMIC (IONIZ.) DATA => ' READ(5,5)INATOM WRITE(6,'(A,$)')'IONIZ. FRACTIONS [0] OR LEVEL POPUL. [1]? => ' READ(5,*)IONLEV IF(IONLEV.EQ.1)THEN WRITE(6,'(A,$)')'FILE WITH LEVEL DATA -KURUCZ LINE FORMAT => ' READ(5,5)INLEV ENDIF WRITE(6,'(A)')' OUTPUT: ION.OUT ' OPEN(7,FILE=INATOM,STATUS='OLD') OPEN(2,FILE='ion.out') C COMMENTS READ(7,'(36X,A30)')INN READ(7,*) READ(7,*)IZMAZ,IZMAZ,IZMAZ,NION IF(NION.GT.NION1)THEN WRITE(6,*)' ERROR, TOO MANY IONIZATION DEGREES CONSIDERED' STOP ENDIF READ (7,'(////////)') READ(7,*)IAT,ZMAZ READ(7,*) C 8-IONIZ.POTENCIALS READ (7,*)(ZIP(I),I=1,8) READ(7,*) READ(7,*)IPART READ(7,*) DO 10 J=1,NION1 C IRWIN COEF. FOR PARTITION F. FOR NION1 STATES READ (7,*)(STA(I,J),I=1,6) 10 CONTINUE CLOSE(7) IF(IONLEV.EQ.1)THEN OPEN(8,FILE=INLEV,STATUS='OLD') READ(8,*)ZMAZ,COD,ZMAZ,ELO,QLO,EUP,QUP GLO=2.D0*QLO+1.D0 GUP=2.D0*QUP+1.D0 STIO=(COD-DINT(COD))*100.D0+1.D0 IION=INT(STIO+0.5D0) CLOSE(8) ENDIF C IF (IMODEL.EQ.1) THEN C READ HUBENY'S MODEL OF THE ATMOSPHERE C DM - MASS DEPTH COORDINATE C TEMP - TEMPERATURE C EKONC - ELECTRON NUMBER DENSITY C AKONC - (ALL) ATOM NUMBER DENSITY C CALL HMODEL(INMOD,NDIM,NBOD P ,DM,TEMP,EKONC,AKONC,DENS) TEFF=0.D0 GLOG=0.D0 ENDIF IF (IMODEL.EQ.0) THEN CALL AMODEL(INMOD,NDIM,NBOD,TEFF,GLOG, P DM,HLB,TEMP,PRESS,EKONC,AKONC,DENS) C MODEL OF THE ATMOSPHERE (NEW KURUCZ) C DM - MASS DEPTH COORDINATE C HLB - GEOMETRICAL DEPTH C TEMP - TEMPERATURE C PRESS - GASS PRESSURE C EKONC - ELECTRON NUMBER DENSITY C AKONC - (ALL) ATOM NUMBER DENSITY C DENS - DENSITY ENDIF WRITE(2,16)TEFF,GLOG,INN 16 FORMAT(2F10.2,5X,A30) C----------------------------------------------------------------------- C PARTITION FUNCTIONS AT DIFFERENT DEPTHS C STAVS(I,J) : I-DEPTH POINT, J- ION IF(IPART.EQ.1)THEN CALL PFDWOR(NDIM,NBOD,NION1,NION,IAT,TEMP,EKONC,STAVS) ELSE CALL PFIRW(NDIM,NBOD,NION1,NION,STA,TEMP,STAVS) ENDIF C----------------------------------------------------------------- C CALCULATION OF THE THREE MOST IMPORTANT IONIZATION STAGES C AT I-TH DEPTH (NUMERICALLY INDENDENT OF OTHER DEPTH POINTS) C ZI(I,1) -MOST POPULATED ION C ZI(I,2) -SECOND MOST POPULATED ION C ZI(I,3) -THIRD MOST POPULATED ION C SA(I,1) -N(MIDDLE)/N(LOWEST IONIZED ION) C SA(I,2) -N(HIGHEST ION)/N(MIDDLE IONIZED ION) C----------------------------------------------------------------- DO 190 I=1,NBOD DO 180 J=1,8 IF (ZIP(J).EQ.0.0) GO TO 180 PSAH=SAHA(TEMP(I),EKONC(I),ZIP(J),STAVS(I,J),STAVS(I,J+1)) IF (PSAH-1.)150,170,180 150 IF (J-1)1200,170,160 160 ZI(I,1)=FLOAT(J) PSAH1=SAHA(TEMP(I),EKONC(I),ZIP(J-1),STAVS(I,J-1) P ,STAVS(I,J)) IF (PSAH.GE.1./PSAH1) THEN ZI(I,2)=FLOAT(J+1) IF(J.GT.7)GOTO 200 PSAH2=SAHA(TEMP(I),EKONC(I),ZIP(J+1),STAVS(I,J+1) P ,STAVS(I,J+2)) IF(PSAH*PSAH2.GE.1./PSAH1)THEN ZI(I,3)=FLOAT(J+2) SA(I,1)=PSAH SA(I,2)=PSAH2 ELSE ZI(I,3)=FLOAT(J-1) SA(I,1)=PSAH1 SA(I,2)=PSAH ENDIF ELSE ZI(I,2)=FLOAT(J-1) IF(J.EQ.2)THEN ZI(I,3)=3. SA(I,1)=PSAH1 SA(I,2)=PSAH GOTO 190 ENDIF PSAH2=SAHA(TEMP(I),EKONC(I),ZIP(J-2),STAVS(I,J-2) P ,STAVS(I,J-1)) IF(1./PSAH.LE.PSAH1*PSAH2)THEN ZI(I,3)=FLOAT(J+1) SA(I,1)=PSAH1 SA(I,2)=PSAH ELSE ZI(I,3)=FLOAT(J-2) SA(I,1)=PSAH2 SA(I,2)=PSAH1 ENDIF ENDIF GO TO 190 170 ZI(I,1)=FLOAT(J) ZI(I,2)=FLOAT(J+1) IF (J.EQ.1)THEN ZI(I,3)=FLOAT(3) SA(I,1)=PSAH SA(I,2)=SAHA(TEMP(I),EKONC(I),ZIP(2),STAVS(I,2) P ,STAVS(I,3)) ELSE IF(J.GT.7)GOTO 200 PSAH1=SAHA(TEMP(I),EKONC(I),ZIP(J-1),STAVS(I,J-1) P ,STAVS(I,J)) PSAH2=SAHA(TEMP(I),EKONC(I),ZIP(J+1),STAVS(I,J+1) P ,STAVS(I,J+2)) IF(PSAH2.GE.1./PSAH1)THEN ZI(I,3)=FLOAT(J+2) SA(I,1)=PSAH SA(I,2)=PSAH2 ELSE ZI(I,3)=FLOAT(J-1) SA(I,1)=PSAH1 SA(I,2)=PSAH ENDIF ENDIF GO TO 190 180 CONTINUE IF (PSAH.GT.1.)GOTO 200 190 CONTINUE GOTO 210 200 WRITE(6,*)'YOU SHOULD CONSIDER HIGHER IONIZATION STAGES' STOP 210 CONTINUE C----------------------------------------------------------------------- C POPULATIONS OF THE THREE MOST IMPORTANT IONS/ TOTAL ELEMENT C POPULATION C F1T(I,1) C F1T(I,2) C F1T(I,3) C----------------------------------------------------------------------- DO 420 I=1,NBOD IF (ZI(I,1).GT.ZI(I,2).AND.ZI(I,1).GT.ZI(I,3))THEN F1T(I,1)=1./(1.+1./SA(I,2)+1./SA(I,1)/SA(I,2)) F1T(I,2)=1./(1.+SA(I,2)+1./SA(I,1)) ELSE IF (ZI(I,1).LT.ZI(I,2).AND.ZI(I,1).LT.ZI(I,3))THEN F1T(I,1)=1./(1.+SA(I,1)+SA(I,1)*SA(I,2)) F1T(I,2)=1./(1.+1./SA(I,1)+SA(I,2)) ELSE IF (ZI(I,2).GT.ZI(I,1))THEN F1T(I,1)=1./(1.+SA(I,2)+1./SA(I,1)) F1T(I,2)=1./(1.+1./SA(I,2)+1./SA(I,1)/SA(I,2)) ELSE F1T(I,1)=1./(1.+1./SA(I,1)+SA(I,2)) F1T(I,2)=1./(1.+SA(I,1)+SA(I,1)*SA(I,2)) ENDIF ENDIF ENDIF 420 CONTINUE DO 430 I=1,NBOD F1T(I,3)=1.D0-F1T(I,1)-F1T(I,2) 430 CONTINUE C----------------------------------------------------------------------- C WRITE(*,*)' DM * TEMP * Z1*Z2*Z3 * N1/N * N2/N * ' C P ,' N3/N' WRITE(2,*)' DM * TEMP * Z1*Z2*Z3 * N1/N * N2/N * ' P ,' N3/N' DO 260 I=1,NBOD C WRITE(*,270)DM(I),TEMP(I),ZI(I,1),ZI(I,2),ZI(I,3),F1T(I,1) C P ,F1T(I,2),F1T(I,3) WRITE(2,270)DM(I),TEMP(I),ZI(I,1),ZI(I,2),ZI(I,3),F1T(I,1) P ,F1T(I,2),F1T(I,3) 260 CONTINUE 270 FORMAT (E10.4,1X,F8.1,1X,3F3.0,3(E12.4)) C----------------------------------------------------------------------- C RR(I,J)= POPULATION OF THE J-TH ION AT I-TH DEPTH/TOTAL C ELEMENT POPULATION (IONIZATION FRACTION) C R1(J)= N(J)/N(1) C RION(J)= N(J)/N(J-1) C----------------------------------------------------------------------- CALL ZASTUP(DM,TEMP,EKONC,ZIP,STAVS,RR,RION,R1,NDIM,NBOD P ,NION1,NION) IF(IONLEV.EQ.1)THEN CALL POPLEV(DM,TEMP,IION,ELO,GLO,EUP,GUP,STAVS,RR P ,NDIM,NBOD,NION1) ENDIF GO TO 1220 1200 WRITE (6,1210) 1210 FORMAT (' ERROR IN THE CODE!!!') 1220 CONTINUE CLOSE(2) STOP END C----------------------------------------------------------------------- SUBROUTINE POPLEV(DM,TEMP,IION,ELO,GLO,EUP,GUP,STAVS,RR P ,NDIM,NBOD,NION1) C RR(I,J)= POPULATION OF THE J-TH ION AT I-TH DEPTH/TOTAL C ELEMENT POPULATION (IONIZATION FRACTION) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION STAVS(NDIM,NION1) P ,DM(NDIM),TEMP(NDIM) P ,RR(NDIM,NION1) WRITE(2,'(A,A)')' DM TEMP N_LOW/N N_UP/N ', P ' N_LOW/N_ION N_UP/N_ION' DO 10 ID1=1,NBOD POPLO=BOLT(TEMP(ID1),STAVS(ID1,IION),GLO,ELO,1)*RR(ID1,IION) POPUP=BOLT(TEMP(ID1),STAVS(ID1,IION),GUP,EUP,1)*RR(ID1,IION) WRITE(2,20)DM(ID1),TEMP(ID1),POPLO,POPUP P ,POPLO/RR(ID1,IION),POPUP/RR(ID1,IION) 20 FORMAT(6(D11.4,1X)) 10 CONTINUE RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION BOLT(TEMP,G1,G2,POT,IRIAD) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C TEMP -TEMPERATURE [KELVIN] C C G1 -STATISTICAL WEIGHT OF THE LOWER LEVEL-m C C G2 -STATISTICAL WEIGHT OF THE UPPER LEVEL-n C C POT -ENERGY DIFFERENCE (UPPER-LOWER LEVEL) [EV], POT>0 C C OUTPUT -POPULATION RATIO N(n)/N(m) C C---------------------OR IF: -------------------------------C C TEMP -TEMPERATURE [KELVIN] C C G2 -STATISTICAL WEIGHT OF THE LEVEL-n C C G1 -PARTITION FUNCTION OF THE ION C C POT -EXCITATION POTENCIAL OF THE n-TH LEVEL [EV] C C (WITH RESPECT TO GROUND STATE), POT>0 C C OUTPUT -RATIO N(n)/N, (n-TH LEVEL POPULATION TO C C TO THE ION POPULATION C C-----------------------OR IF------------------------------------------C C IRIAD=1 POT IS IN [1/CM] C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(IRIAD.EQ.1)THEN BOLT=G2/(G1*DEXP(2.302585D0*5039.77D0/TEMP*POT/8065.9D0)) RETURN ENDIF BOLT=G2/(G1*DEXP(2.302585D0*5039.77D0/TEMP*POT)) RETURN END C----------------------------------------------------------------------- SUBROUTINE ZASTUP(DM,TEMP,EKONC,ZIP,STAVS,RR,RION,R1,NDIM,NBOD P ,NION1,NION) C RR(I,J)= POPULATION OF THE J-TH ION AT I-TH DEPTH/TOTAL C ELEMENT POPULATION (IONIZATION FRACTION) C R1(J)= N(J)/N(1) C RION(J)= N(J)/N(J-1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION ZIP(NION1-1),STAVS(NDIM,NION1) P ,DM(NDIM),TEMP(NDIM),EKONC(NDIM) P ,RION(NION1),R1(NION1),RR(NDIM,NION1) DO 10 J=1,NION1 DO 5 I=1,NBOD RR(I,J)=1.D-31 5 CONTINUE 10 CONTINUE DO 40 I=1,NBOD RION(1)=1 R1(1)=1 SUM1=1.D0 DO 20 J=2,NION RION(J)=SAH(TEMP(I),EKONC(I),ZIP(J-1),STAVS(I,J-1) P ,STAVS(I,J)) R1(J)=R1(J-1)*RION(J) SUM1=SUM1+R1(J) 20 CONTINUE DO 30 J=1,NION RR(I,J)=R1(J)/SUM1 30 CONTINUE 40 CONTINUE WRITE(2,'(A)')' MORE PRECISE IONIZATION FRACTIONS' WRITE(2,'(A,40I11)')' DM',(I,I=1,NION) DO 50 I=1,NBOD WRITE(2,130)DM(I),(RR(I,J),J=1,NION) 130 FORMAT(40(D10.3,1X)) 50 CONTINUE RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION SAH(TE,EK,ZP,S1,S2) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C SAHA EQUATION (RATIO OF THE TWO NEIGHBOUR ION NUMBER C DENSITIES) C ZP -IV [EV],BOLT.CONST.-IN [ERG/K] SAH=1.2020D9*S2/(S1*(5039.77D0/TE)**2.5D0 P *DEXP(2.302585D0*5039.77D0/TE*ZP)*EK*TE*1.38062D-16) RETURN END C----------------------------------------------------------------------- SUBROUTINE AMODEL(INN,NM1,NM,TEFF,GLOG P ,VDM,VHLB,VTEP,VPRESS,VEKONC,VAKONC,VDENS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION VHLB(NM1),VAKONC(NM1),VEKONC(NM1),VPRESS(NM1) p ,VTEP(NM1),VDM(NM1),VDENS(NM1) CHARACTER*30 INN C C MODEL OF THE ATMOSPHERE (NEW KURUCZ) C C DM - MASS DEPTH COORDINATE C HLB - GEOMETRICAL DEPTH C TEP - TEMPERATURE C PRESS - GASS PRESSURE C EKONC - ELECTRON NUMBER DENSITY C AKONC - (ALL) ATOM NUMBER DENSITY C DENS - DENSITY C OPEN(3,FILE=INN,STATUS='OLD') READ(3,30) TEFF,GLOG READ(3,35) NM 30 FORMAT(4X,F8.0,9X,F8.5) 35 FORMAT(/////////////////////10X,I3) BOL=1.38062D-16 C WMM=SUM_e (M_e*ABUNDANCE_e) IS MEAN THE ATOM MASS C WHERE ABUNDANCE_e IS ELEMENT NUMBER DENSITY/AKONC C FOR 90% H + 10 % HE WMM=2.2034D-24 VHLB(1)=0.D0 DO 40 I=1,NM READ(3,*)VDM(I),VTEP(I),VPRESS(I),VEKONC(I) VAKONC(I)=VPRESS(I)/VTEP(I)/BOL VAKONC(I)=VAKONC(I)-VEKONC(I) VDENS(I)=VAKONC(I)*WMM IF(I.GT.1) THEN DELTA=2.D0*(VDM(I)-VDM(I-1))/(VDENS(I)+VDENS(I-1)) VHLB(I)=DELTA+VHLB(I-1) ENDIF 40 CONTINUE CLOSE(3) RETURN END C----------------------------------------------------------------------- SUBROUTINE HMODEL(INN,NM1,NDPTH P ,DEPTH,TEMP,ELEC,AKONC,DENS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION DEPTH(NM1),TEMP(NM1),ELEC(NM1), P AKONC(NM1),DENS(NM1),X(500) CHARACTER*30 INN C C MODEL OF THE ATMOSPHERE (HUBENY) C C DEPTH - MASS DEPTH COORDINATE C TEMP - TEMPERATURE C ELEC - ELECTRON NUMBER DENSITY [CGS] C AKONC - (ALL) ATOM NUMBER DENSITY [CGS] C DENS - DENSITY [CGS] C C WMM=SUM_e (M_e*ABUNDANCE_e) IS THE MEAN ATOM MASS C WHERE ABUNDANCE_e IS N_e/AKONC C FOR 90% H + 10 % HE : WMM=2.2034D-24 OPEN(3,FILE=INN,STATUS='OLD') READ(3,*) NDPTH,NUMPAR READ(3,*) (DEPTH(I),I=1,NDPTH) DO 30 ID=1,NDPTH READ(3,*) (X(I),I=1,NUMPAR) TEMP(ID)=X(1) ELEC(ID)=X(2) DENS(ID)=X(3) AKONC(ID)=DENS(ID)/WMM C write(9,'(4e11.4)')depth(id),temp(id),elec(id),dens(id) 30 CONTINUE CLOSE(3) RETURN END C----------------------------------------------------------------------- SUBROUTINE PFIRW(NDIM,NBOD,NION1,NION,STA,TEMP,STAVS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER(NION1=9) DIMENSION STA(6,NION1),STAVS(NDIM,NION1),TEMP(NDIM) DO 80 I=1,NBOD DO 70 J=1,NION IF (ABS(STA(1,J)-0.D0).LT.1.D-30) GO TO 60 SALA=DLOG(TEMP(I)) ARG=STA(1,J)+STA(2,J)*SALA+STA(3,J)*SALA*SALA+ P STA(4,J)*SALA**3.D0+STA(5,J)*SALA**4.D0+STA(6,J)*SALA**5.D0 STAVS(I,J)=DEXP(ARG) GO TO 70 60 STAVS(I,J)=1.D0 70 CONTINUE 80 CONTINUE WRITE (2,'(A,A)')' TEM STAVS1 STAVS2 STAVS3 STAVS4' P,' STAVS5 STAVS6 STAVS7 STAVS8 STAVS9' DO 140 I=1,NBOD WRITE(2,130)TEMP(I),(STAVS(I,J),J=1,NION) 130 FORMAT(F7.0,20E10.3) 140 CONTINUE RETURN END C----------------------------------------------------------------------- subroutine potlow(TEMP,ane,eplow) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c TEMP - Temperature[K] c ane - electron number density [cm^{-3}] c lowering of the ionization potential for neutrals, z=1 in [eV] c for ions(z): eplow(z)=z*eplow(1) BOL=1.38062D-16 PI=3.14159D0 e=4.80325d-10 pom=bol*TEMP/8.d0/pi/ane c debye lenght in [CGS] debye=sqrt(pom)/e eplow=1.44d-7/debye RETURN END C----------------------------------------------------------------------- SUBROUTINE PFDWOR(NDIM,NBOD,NION1,NION,IAT,TEMP,EKONC,STAVS) c c partition functions after: c M.M.Dworetsky & B. Smalley, private communication c IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION STAVS(NDIM,NION1),EKONC(NDIM),TEMP(NDIM),U(5) DO 80 I=1,NBOD CALL POTLOW(TEMP(I),EKONC(I),EPLOW) THETA=5039.77D0/TEMP(I) C PARTFN IS IN 'pfdwor.inc' CALL PARTFN(IAT,THETA,U,EPLOW) DO 70 J=1,NION IF(J.LE.5)THEN STAVS(I,J)=U(J) ELSE STAVS(I,J)=1.D0 ENDIF 70 CONTINUE 80 CONTINUE WRITE (2,'(A,A)')' TEM STAVS1 STAVS2 STAVS3 STAVS4' P,' STAVS5 STAVS6 STAVS7 STAVS8 STAVS9' DO 140 I=1,NBOD WRITE(2,130)TEMP(I),(STAVS(I,J),J=1,NION) 130 FORMAT(F7.0,20E10.3) 140 CONTINUE RETURN END C----------------------------------------------------------------------- INCLUDE 'pfdwor.inc'