SUBROUTINE neispec(ntemp, temp, nion, ionfrac, ionel, ionstage, & ne, ear, z, linesh, nabun, abund, abunel, & qcalc, callr, photar) IMPLICIT NONE INTEGER ne, nion, nabun, ntemp REAL temp(ntemp), ionfrac(nion,ntemp), ear(0:ne), z, linesh(2) REAL abund(nabun), photar(ne) INTEGER ionel(nion), ionstage(nion), abunel(nabun) CHARACTER*(*) callr LOGICAL qcalc c Routine to return the spectrum from a multi-temperature NEI plasma. c This is a wrap-up routine to call either nneispec (the default) or c oneispec if NEIVERS is set to 1.x. c c Arguments c temp r i: temperature (keV) c nion i i: number of different ions c ionfrac r i: fraction in ion c ionel i i: element for this ion c ionstage i i: level for this ion (ionel+1 = fully stripped) c ne i i: number of response energies c ear r i: response energies c z r i: redshift c linesh r i: lineshape parameters c nabun i i: number of elements with abundances c abund r i: abundances c abunel i i: element for this abundance c qcalc l i: if true force a recalculation of all arrays c callr c i: name of calling routine c photar r r: output spectrum INTEGER TOTEL PARAMETER(TOTEL=30) REAL andgrv(TOTEL) INTEGER lenn, iold, newnion, newnabun, j CHARACTER pname*128, pvalue*128, ovalue*128, tmpstr*255 CHARACTER*2 celt(TOTEL) LOGICAL qdone INTEGER lenact REAL fgabnd CHARACTER fgmstr*128 EXTERNAL lenact, fgmstr, fgabnd SAVE newnion, newnabun, ovalue DATA celt/'H ', 'He', 'Li', 'Be', 'B ', 'C ', 'N ', 'O ', 'F ', & 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P ', 'S ', 'Cl', 'Ar', & 'K ', 'Ca', 'Sc', 'Ti', 'V ', 'Cr', 'Mn', 'Fe', 'Co', & 'Ni', 'Cu', 'Zn'/ DATA andgrv/1.00e+00, 9.77e-02, 1.45e-11, 1.41e-11, 3.98e-10, & 3.63e-04, 1.12e-04, 8.51e-04, 3.63e-08, 1.23e-04, & 2.14e-06, 3.80e-05, 2.95e-06, 3.55e-05, 2.82e-07, & 1.62e-05, 3.16e-07, 3.63e-06, 1.32e-07, 2.29e-06, & 1.26e-09, 9.77e-08, 1.00e-08, 4.68e-07, 2.45e-07, & 4.68e-05, 8.32e-08, 1.78e-06, 1.62e-08, 3.98e-08/ DATA ovalue / '1.1' / c Check for the NEIVERS string parameter. pname = 'NEIVERS' pvalue = fgmstr(pname) IF ( pvalue .NE. ' ' .AND. pvalue .NE. ovalue ) THEN ovalue = pvalue tmpstr = 'Changing to NEI models version '// & pvalue(:lenact(pvalue)) CALL xwrite(tmpstr,10) ENDIF tmpstr = 'Running NEI models version '//ovalue(:lenact(ovalue)) CALL xwrite(tmpstr, 20) c Convert from whatever the current abundance ratios are to the Anders & c Grevesse ratios which are used internally DO j = 1, nabun abund(j) = abund(j) * fgabnd(celt(abunel(j))) & / andgrv(abunel(j)) ENDDO IF ( ovalue(1:1) .EQ. '2' ) THEN CALL nneispec(ntemp, temp, nion, ionfrac, ionel, ionstage, & ne, ear, z, linesh, nabun, abund, abunel, & qcalc, callr, photar) ELSEIF ( ovalue(1:1) .EQ. '1' ) THEN c The old version which had no Ar so we have to strip this out of the ion* c and ab* arrays (they may already have been removed from the ion* arrays c if they have not been recalculated since the last invocation) qdone = .TRUE. DO iold = 1, nion IF ( ionel(iold) .EQ. 18 ) qdone = .FALSE. ENDDO IF ( .NOT.qdone ) THEN newnion = 0 DO iold = 1, nion IF ( ionel(iold) .NE. 18 ) THEN newnion = newnion + 1 DO j = 1, ntemp ionfrac(newnion, j) = ionfrac(iold, j) ENDDO ionel(newnion) = ionel(iold) ionstage(newnion) = ionstage(iold) ENDIF ENDDO ENDIF newnabun = 0 DO iold = 1, nabun IF ( abunel(iold) .NE. 18 ) THEN newnabun = newnabun + 1 abund(newnabun) = abund(iold) abunel(newnabun) = abunel(iold) ENDIF ENDDO CALL oneispec(ntemp, temp, nion, ionfrac, ionel, ionstage, & ne, ear, z, linesh, newnabun, abund, abunel, & qcalc, callr, photar) ENDIF RETURN END c ************************************************************************ SUBROUTINE nneispec(ntemp, temp, nion, ionfrac, ionel, ionstage, & ne, ear, z, linesh, nabun, abund, abunel, & qcalc, callr, photar) INTEGER ne, nion, nabun, ntemp REAL temp(ntemp), ionfrac(nion,ntemp), ear(0:ne), z, linesh(2) REAL abund(nabun), photar(ne) INTEGER ionel(nion), ionstage(nion), abunel(nabun) CHARACTER*(*) callr LOGICAL qcalc c Routine to return the spectrum from a multi-temperature NEI plasma. c This uses the APEC tables. This routine loads the APEC data arrays c then calls xneispec to actually do the work. c Arguments c temp r i: temperature (keV) c nion i i: number of different ions c ionfrac r i: fraction in ion c ionel i i: element for this ion c ionstage i i: level for this ion (ionel+1 = fully stripped) c ne i i: number of response energies c ear r i: response energies c z r i: redshift c linesh r i: lineshape parameters c nabun i i: number of elements with abundances c abund r i: abundances c abunel i i: element for this abundance c qcalc l i: if true force a recalculation of all arrays c callr c i: name of calling routine c photar r r: output spectrum INCLUDE '../../inc/xspec.inc' C Pointers to local arrays C itval Tabulated temperatures C icont Tabulated continua C icenerg Energies for tabulated continua C iconEl Atomic numbers for tabulated continua c iconIon Ion stages for tabulated continua C iconpt Offsets in the cont array for each ion and kT C iconsiz Sizes in the cont array for each ion and kT C ipseudo Tabulated pseudo-continua C ipenerg Energies for tabulated pseudo-continua C ipsupt Offsets in the pseudo array for each ion and kT C ipsusiz Sizes in the pseudo array for each ion and kT C ilineE Tabulated line energies C ilineEm Tabulated line emissivities C ilineEl Tabulated line element (by atomic number) c ilineIon Tabulated line ion C ilinpt Offsets in the line arrays for each kT C ilinsiz Sizes in the line arrays for each kT c nmtval The number of tabulated temperatures c nmions The number of tabulated ions INTEGER itval, icont, iconEl, iconIon, icenerg, iconpt INTEGER iconsiz, ipseudo, ipenerg, ipsupt, ipsusiz, ilineE INTEGER ilineEm, ilineEl, ilineIon, ilinpt, ilinsiz REAL velocity INTEGER nmtval, nmions INTEGER status, lenp CHARACTER contxt*72, wrtstr*255 CHARACTER*128 pname, pvalue, apecroot LOGICAL qfirst, qtherm SAVE qfirst, nmtval, nmions, apecroot SAVE itval, icont, iconEl, iconIon, icenerg, iconpt, iconsiz SAVE ipseudo, ipenerg, ipsupt, ipsusiz, ilineE, ilineEm, ilineEl SAVE ilineIon, ilinpt, ilinsiz INTEGER lenact CHARACTER*128 fgmstr EXTERNAL fgmstr, lenact DATA qfirst, apecroot / .TRUE., ' ' / status = 0 C If the user has switched NEI APEC files then we reset the qfirst logical pname = 'NEIAPECROOT' pvalue = fgmstr(pname) IF ( pvalue .NE. apecroot ) THEN qfirst = .TRUE. apecroot = pvalue ENDIF C If the first time through then load the model file. If this file has C already been loaded then LDNEIAPED just returns the pointers to the C appropriate memory. IF ( qfirst ) THEN CALL ldneiaped(itval, nmtval, nmions, icont, icenerg, iconEl, & iconIon, iconpt, iconsiz, ipseudo, ipenerg, & ipsupt, ipsusiz, ilineE, ilineEm, ilineEl, & ilineIon, ilinpt, ilinsiz, status) contxt = 'Failure in LDNEIAPED' IF ( status .NE. 0 ) GOTO 999 qfirst = .FALSE. ENDIF C Check for values of APECTHERMAL and APECVELOCITY xset parameters c Check the APECTHERMAL variable - if the variable has not been set then c assume .false. pname = 'APECTHERMAL' pvalue = fgmstr(pname) lenp = lenact(pvalue) CALL upc(pvalue) qtherm = .FALSE. IF ( lenp .GT. 0 ) THEN IF ( (pvalue(:lenp) .EQ. 'YES' .OR. & pvalue(:lenp) .EQ. 'ON') ) THEN qtherm = .TRUE. ELSEIF ((pvalue(:lenp) .EQ. 'NO' .OR. & pvalue(:lenp) .EQ. 'OFF') ) THEN qtherm = .FALSE. ENDIF ENDIF c and the APECVELOCITY variable - if it doesn't exist then use 0. pname = 'APECVELOCITY' pvalue = fgmstr(pname) lenp = lenact(pvalue) velocity = 0.0 IF ( lenp .GT. 0 ) THEN READ(pvalue(:lenp), *, iostat=ierr) velocity IF ( ierr .NE. 0 ) THEN WRITE(wrtstr,*) & 'Failed to read value from APECVELOCITY - assuming 0.' CALL xwrite(wrtstr, 10) velocity = 0.0 ENDIF ENDIF WRITE(wrtstr, *) 'Thermal broadening : ', qtherm CALL xwrite(wrtstr, 25) WRITE(wrtstr, *) 'Velocity broadening : ', velocity CALL xwrite(wrtstr, 25) c Call the routine to do the work CALL xneispec(ntemp, temp, nion, ionfrac, ionel, ionstage, & ne, ear, z, linesh, nabun, abund, abunel, & qcalc, callr, photar, nmtval, MEMR(itval), nmions, & MEMR(icont), MEMR(icenerg), MEMI(iconEl), & MEMI(iconIon), MEMI(iconpt), MEMI(iconsiz), & MEMR(ipseudo), MEMR(ipenerg), MEMI(ipsupt), & MEMI(ipsusiz), MEMR(ilineE), MEMR(ilineEm), & MEMI(ilineEl), MEMI(ilineIon), MEMI(ilinpt), & MEMI(ilinsiz), qtherm, velocity, status) contxt = 'Failure in XNEISPEC' 999 CONTINUE IF ( status .NE. 0 ) THEN WRITE(wrtstr,'(''Error '', i6, '' in NNEISPEC'')') status CALL xwrite(wrtstr, 10) CALL xwrite(contxt, 10) ENDIF RETURN END c ************************************************************************ SUBROUTINE xneispec(ntemp, temp, nion, ionfrac, ionel, ionstage, & ne, ear, z, linesh, nabun, abund, abunel, & qcalc, callr, photar, nmtval, tval, nmions, & cont, cenerg, conEl, conIon, conpt, & consiz, pseudo, penerg, psupt, psusiz, lineE, & lineEm, lineEl, lineIon, linpt, linsiz, & qtherm, velocity, status) IMPLICIT NONE INTEGER ne, nion, nabun, ntemp, nmtval, nmions, status REAL temp(ntemp), ionfrac(nion,ntemp), ear(0:ne), z, linesh(2) REAL abund(nabun), photar(ne) REAL tval(nmtval), velocity REAL cont(*), cenerg(*), pseudo(*), penerg(*) REAL lineE(*), lineEm(*) INTEGER lineEl(*), lineIon(*) INTEGER conpt(nmions, nmtval), consiz(nmions, nmtval) INTEGER psupt(nmions, nmtval), psusiz(nmions, nmtval) INTEGER conEl(nmions, nmtval), conIon(nmions, nmtval) INTEGER linpt(nmtval), linsiz(nmtval) INTEGER ionel(nion), ionstage(nion), abunel(nabun) CHARACTER*(*) callr LOGICAL qcalc, qtherm c Routine to return the spectrum from a multi-temperature NEI plasma. c Arguments c ntemp i i: number of temperatures c temp r i: temperature (keV) c nion i i: number of different ions c ionfrac r i: fraction in ion c ionel i i: element for this ion c ionstage i i: level for this ion (ionel+1 = fully stripped) c ne i i: number of response energies c ear r i: response energies c z r i: redshift c linesh r i: lineshape parameters c nabun i i: number of elements with abundances c abund r i: abundances c abunel i i: element for this abundance c qcalc l i: if true force a recalculation of all arrays c callr c i: name of calling routine c photar r r: output spectrum c nmtval I i: number of tabulated temperatures C tval R i: tabulated temperatures C nmions I i: number of tabulated ions C cont R i: Tabulated continua C cenerg R i: Energies for tabulated continua C conEl I i: Atomic numbers for tabulated continua c conIon I i: Ion stages for tabulated continua C conpt I i: Offsets in the cont array for each ion and kT C consiz I i: Sizes in the cont array for each ion and kT C pseudo R i: Tabulated pseudo-continua C penerg R i: Energies for tabulated pseudo-continua C psupt I i: Offsets in the pseudo array for each ion and kT C psusiz I i: Sizes in the pseudo array for each ion and kT C lineE R i: Tabulated line energies C lineEm R i: Tabulated line emissivities C lineEl I i: Tabulated line elements c lineIon I i: Tabulated line ion C linpt I i: Offsets in the line arrays for each kT C linsiz I i: Sizes in the line arrays for each kT C qtherm L i: If true then apply thermal broadening to lines C velocity R i: Gaussian velocity broadening to apply C status I r: 0==OK c This code is adapted from SUMAPE but loops over ions rather than elements. REAL c_cm_s, keV_erg, amu_g, km_cm c Speed of light in cm s^-1 PARAMETER (c_cm_s=2.9979246e10) c kT in erg corresponding to temperature of 1 keV PARAMETER (keV_erg=1.60217646e-9) c Unified atomic mass constant in g PARAMETER (amu_g=1.660538e-24) c km in cm PARAMETER (km_cm=1e5) INTEGER NUMZ PARAMETER (NUMZ=28) REAL abunZ(NUMZ), atomM(NUMZ) REAL tdelta, tcoeff(2), energy, edelta, zf, coeff REAL width, lineflux INTEGER atomZ(13) INTEGER itemp(2), lim(2) INTEGER ien, it, ielt, itabt, il, ihigh, ilow, ion INTEGER ionout, i, itpt CHARACTER wrtstr*72 SAVE atomZ, atomM DATA atomZ /2, 6, 7, 8, 10, 12, 13, 14, 16, 18, 20, 26, 28/ c Atomic masses in AMU DATA atomM /1.00794, 4.002602, 6.941, 9.012182, 10.811, & 12.0107, 14.0067, 15.9994, 18.9984032, 20.1797, & 22.989770, 24.3050, 26.981538, 28.0855, 30.973761, & 32.065, 35.453, 39.948, 39.0983, 40.078, & 44.955910, 47.867, 50.9415, 51.9961, 54.938049, & 55.845, 58.933200, 58.6934/ status = 0 c Initialize output array DO ien = 1, ne photar(ien) = 0.0 ENDDO c set up the abundances in terms of atomic number so we can match them c to the lineEl array DO ielt = 1, NUMZ abunZ(ielt) = 0.0 ENDDO DO ielt = 1, nabun abunZ(abunel(ielt)) = abund(ielt) ENDDO c Loop over required temperatures DO it = 1, ntemp c Find tabulated temperatures to be used in the interpolation and c calculate interpolation coefficients IF ( temp(it) .LT. tval(1) ) THEN WRITE(wrtstr, '(a,g12.5)') & ' Warning: temperature < minimum in model file: ', tval(1) CALL xwrite(wrtstr, 10) itemp(1) = 1 itemp(2) = 1 ELSEIF ( temp(it) .GT. tval(nmtval) ) THEN WRITE (wrtstr, '(a,g12.5)') & ' Warning: temperature > maximum in model file: ', & tval(nmtval) CALL xwrite(wrtstr, 10) itemp(1) = nmtval itemp(2) = nmtval ELSE itemp(2) = 1 DO WHILE ( temp(it) .GE. tval(itemp(2)) ) itemp(2) = itemp(2) + 1 ENDDO itemp(1) = itemp(2) - 1 IF ( temp(it) .EQ. tval(itemp(1)) ) itemp(2) = itemp(1) ENDIF tdelta = tval(itemp(2)) - tval(itemp(1)) IF ( tdelta .GT. 0. ) THEN tcoeff(1) = (tval(itemp(2)) - temp(it))/tdelta tcoeff(2) = (temp(it) - tval(itemp(1)))/tdelta ELSE tcoeff(1) = 1.0 tcoeff(2) = 0.0 ENDIF WRITE(wrtstr, '(a, f8.4, a, f8.4, a, f8.4)' ) & 'Interpolating between ', tval(itemp(1)), & ' and ', tval(itemp(2)), ' for ', temp(it) CALL xwrite(wrtstr, 15) WRITE(wrtstr, '(a, f8.4, a, f8.4)' ) & 'with coefficients ', tcoeff(1), ' and ', tcoeff(2) CALL xwrite(wrtstr, 15) c Loop over the two temperatures to use in interpolation DO itabt = 1, 2 itpt = itemp(itabt) c Loop over tabulated ions. There are at most nmions - if there are less c then the arrays are set to -1. DO ion = 1, nmions IF ( conEl(ion,itpt) .EQ. -1 ) GOTO 100 coeff = tcoeff(itabt) * abunZ(conEl(ion,itpt)) c If we are not using this element (ie abunZ is zero then jump to end c of loop IF ( coeff .EQ. 0. ) GOTO 100 c Find which entry in the ionfrac array this ion matches ionout = 0 DO i = 1, nion IF ( (ionel(i) .EQ. conEl(ion,itpt)) .AND. & (ionstage(i) .EQ. conIon(ion,itpt)) ) ionout = i ENDDO IF ( ionout .EQ. 0 ) GOTO 100 coeff = coeff * ionfrac(ionout,it) c Interpolate the continuum onto the output energy array. Note that the c input is flux/keV values at the energies tabulated so need to multiply c by the bin width to convert to the units that XSPEC requires. Also note c that this algorithm assumes that the response energies are monotonically c increasing. lim(1) = conpt(ion, itpt) lim(2) = lim(1) + consiz(ion, itpt) - 1 ihigh = lim(1) DO ien = 1, ne energy = ( (ear(ien-1) + ear(ien))/2. ) * (1 + z) IF ( energy .GE. cenerg(lim(1)) .AND. & energy .LE. cenerg(lim(2)) ) THEN DO WHILE ( energy .GT. cenerg(ihigh) .AND. & ihigh .LE. lim(2) ) ihigh = ihigh + 1 ENDDO ilow = ihigh - 1 photar(ien) = photar(ien) + coeff * & ( (cenerg(ihigh)-energy)*cont(ilow) & +(energy-cenerg(ilow))*cont(ihigh) ) & * ( ear(ien) - ear(ien-1) ) / (1 + z) & / (cenerg(ihigh)-cenerg(ilow)) ENDIF ENDDO c Repeat for the pseudo-continuum lim(1) = psupt(ion, itpt) lim(2) = lim(1) + psusiz(ion, itpt) - 1 ihigh = lim(1) DO ien = 1, ne energy = ( (ear(ien-1) + ear(ien))/2. ) * (1 + z) IF ( energy .GE. penerg(lim(1)) .AND. & energy .LE. penerg(lim(2)) ) THEN DO WHILE ( energy .GT. penerg(ihigh) .AND. & ihigh .LE. lim(2) ) ihigh = ihigh + 1 ENDDO ilow = ihigh - 1 photar(ien) = photar(ien) + coeff * & ( (penerg(ihigh)-energy)*pseudo(ilow) & +(energy-penerg(ilow))*pseudo(ihigh) ) & * ( ear(ien) - ear(ien-1) ) / (1 + z) & / (penerg(ihigh)-penerg(ilow)) ENDIF ENDDO c End loop over elements. 100 CONTINUE ENDDO c Now do the lines. Loop over the lines for this temperature DO il = linpt(itpt), linpt(itpt)+linsiz(itpt)-1 IF ( abunZ(lineEl(il)) .NE. 0. ) THEN energy = lineE(il) / (1 + z) c Find the width of the line using the equation c W = E/c sqrt [ 2 kT / (A mp ) + velocity^2 ] if thermal broadening c is on else c W = (E/c) velocity IF ( qtherm ) THEN width = energy / c_cm_s * & sqrt( 2. * temp(it) * keV_erg / & atomM( lineEl(il) ) / amu_g + & velocity*velocity*km_cm*km_cm) ELSE width = (energy / c_cm_s) * velocity * km_cm ENDIF c Find which entry in the ionfrac array this line is from ionout = 0 DO i = 1, nion IF ( (ionel(i) .EQ. lineEl(il)) .AND. & (ionstage(i) .EQ. lineIon(il)) ) ionout = i ENDDO c Add in multiplying by abundance and ion fraction IF ( ionout .GT. 0 ) THEN lineflux = tcoeff(itabt) * abunZ(lineEl(il)) & * lineEm(il) * ionfrac(ionout,it) CALL gaussline(ear, ne, energy, width, & lineflux, 3.0, .TRUE., photar) ENDIF ENDIF c End loop over lines. ENDDO c End loop over two temperatures used in the interpolation ENDDO c End loop over temperatures. ENDDO c Fix the normalization by multiplying by 1e14 DO ien = 1, ne photar(ien) = photar(ien) * 1.e14 ENDDO c Correct for time dilation IF ( z .NE. 0 ) THEN zf = 1./(1.+z) DO ien = 1, ne photar(ien) = photar(ien) * zf ENDDO ENDIF RETURN END