SUBROUTINE acisabs(ear, ne, param, ifl, photar, photer) INTEGER ne, ifl REAL*4 ear(0:ne), param(8), photar(ne), photer(ne) c*** Ver 1.1 - Modified the previous ver. of acisabs.f c*** to use exponential decay instead of linear c FORTRAN version of xcont.pro, written by George Chartas c Conversion performed by Konstantin Getman 07/17/02 c c THE CONTENT OF xcont.pro c (1) Use observed decay with time in the ratio of observed c X-ray transmissions at 0.67keV and 5.895keV. c Observed decay in S3 based on Catherine Grants' analysis: c Fit parameters provided by Allyn Tennant c c R(t) = (norm*exp(-tauinf*(1.0-exp(-t/tefold)))) eq(1) c norm=0.00722+/-0.00007 c tefold=620.+/-66. c tauinf=0.582+/-0.024 c c c (2) Model decay c c Trans(E1,t) = TOBF(E1) * Tcont(E1,t) eq(2) c c where TOBF(E1) is the transmission of the OBF at energy E1 c Tcont(E1,t) is the transmission of the contamination c layer at energy E1 and time t c c The modeled decay is: c c R(t)=Trans(E1,t)/Trans(E2,t)=const*exp{-rho*d(t)*[mac(E1)-mac(E2)]} eq(3) c c c where : rho is the density of the contaminant c d(t) is the thickness of the contaminant c mac(E) is the mass absorption coefficient of the c contaminant at energy E. c mac(E) calculated from atomic scattering factor files c provided at http://www-cxro.lbl.gov/optical_constants/asf.html c c The plot of the thickness of the contaminant vs. time is c derived form eq(1) and eq(3): c c d(t)*rho = ln(R(t)/const) / [mac(E2) - mac(E1)] c c c c (3) Suggested modification of ACIS arf files to account for contamination c c Tcont(E,t) = exp(-mac(E)*rho*d(t)) c c------------------------------------------------------------------------ c THE CONTENT OF acisabs xspec model c Number of model parameters: 7 c 1 Tdays Days between launch and observation c 2 norm Normalization factor c 3 tauinf Optical depth at infinite time c 4 tefold e-folding time for build-up of contaminant c 5 nC Number of carbon atoms in hydrocarbon molecule c 6 nH Number of hydrogen atoms in hydrocarbon molecule c 7 nO Number of oxygen atoms in hydrocarbon molecule c 8 nN Number of nitrogen atoms in hydrocarbon molecule c c THIS PROGRAM REQUIRES the following data files to be present in the c $LMODDIR directory : c c.nff, h.nff, o.nff, n.nff c c Declaration integer*4 nmax, numaxC, numaxO, numaxH,numaxN, i real*4 awC, awO, awH, awN, Mac, MacE1, MacE2, thRho, decay, log real*4 Tdays, nC, nH, nO, nN parameter (nmax=1000, Ecal1=0.67, Ecal2=5.895, awC=12.01115, + awO=15.9994, awH=1.00797, awN=14.0067) real*4 tauinf, tefold, norm, muO(nmax), muH(nmax), + muC(nmax), muN(nmax), enO(nmax), enN(nmax), enH(nmax), enC(nmax), + xO, xN, xH, xC, elnorm, eninpt logical qfirst save qfirst, numaxO, enO, muO, numaxH, enH, muH save numaxN, enN, muN, numaxC, enC, muC data qfirst / .true. / c Initialization Tdays = param(1) norm = param(2) tauinf = param(3) tefold = param(4) nC = param(5) nH = param(6) nO = param(7) nN = param(8) c Read data IF ( qfirst ) THEN CALL READata ('o',awO,numaxO,enO,muO) CALL READata ('h',awH,numaxH,enH,muH) CALL READata ('n',awN,numaxN,enN,muN) CALL READata ('c',awC,numaxC,enC,muC) qfirst = .FALSE. ENDIF c Calculate the relative abundance of each element xO = nO * awO xH = nH * awH xN = nN * awN xC = nC * awC elnorm = xO + xH + xN + xC xO = xO/elnorm xH = xH/elnorm xN = xN/elnorm xC = xC/elnorm c Calculate mass absorption coefficient at Ecal1 and Ecal2 MacE1 = xO*Mac(numaxO,enO,muO,Ecal1)+xH*Mac(numaxH,enH,muH,Ecal1) ++xN*Mac(numaxN,enN,muN,Ecal1)+xC*Mac(numaxC,enC,muC,Ecal1) MacE2 = xO*Mac(numaxO,enO,muO,Ecal2)+xH*Mac(numaxH,enH,muH,Ecal2) ++xN*Mac(numaxN,enN,muN,Ecal2)+xC*Mac(numaxC,enC,muC,Ecal2) c Estimate thickness of contaminate during observation decay = norm*exp(-tauinf*(1.0 - exp(-Tdays/tefold))) thRho = log(decay/norm)/(MacE2 - MacE1) c Calculate transmission of contamination layer for input energies DO 105 i = 1, ne eninpt = (ear(i-1) + ear(i))/2. photar(i) = exp (-( + xO*Mac(numaxO,enO,muO,eninpt)+xH*Mac(numaxH,enH,muH,eninpt) + +xN*Mac(numaxN,enN,muN,eninpt)+xC*Mac(numaxC,enC,muC,eninpt + ))*thRho) 105 CONTINUE RETURN END c Additional Subroutines and Functions c----------------------------------------- c This one reads f2 and E from nff file c and calculates mu(energy) array for the specified element SUBROUTINE READata (elem,aw,numax,engrid,mu) implicit none integer numax, nnmax, ilun, lenn, ios character elem*1, contxt*300, fl_name*300 PARAMETER (nnmax=1000) real*4 aw, engrid(nnmax), mu(nnmax) INTEGER lenact EXTERNAL lenact CALL getlun(ilun) CALL getenv('LMODDIR', fl_name) lenn = lenact(fl_name) fl_name = fl_name(:lenn) // '/' // elem // '.nff' ios = 0 numax=0 CALL openwr(ilun, fl_name, 'old', ' ', ' ', 0, 0, ios) IF ( ios .NE. 0 ) THEN contxt = 'Failed to open '//fl_name(:lenact(fl_name)) CALL xwrite(contxt, 10) CALL frelun(ilun) RETURN ENDIF 15 numax=numax+1 READ(ilun,'(E12.6, 1X, E12.6)',END=25) engrid(numax),mu(numax) mu(numax) = 4.208e7 * mu(numax) / engrid(numax) / aw GOTO 15 25 CONTINUE numax=numax-1 CLOSE (ilun) CALL frelun(ilun) RETURN END c--------------------------------------- c--------------------------------------- c Function performs linear interpolation, calculates c mass absorption coefficient for specified energy enkeV FUNCTION Mac(N,engrid,mu,enkeV) integer N REAL*4 engrid(N),mu(N),enkeV, eneV, Mac eneV = 1000. * enkeV DO 100 i = 1, N-1 if ((eneV.GE.engrid(i)) .AND. (eneV.LE.engrid(i+1))) GOTO 101 100 CONTINUE 101 Mac = mu(i+1) - (engrid(i+1)-eneV)*(mu(i+1)-mu(i))/ + (engrid(i+1)-engrid(i)) RETURN END c---------------------------------------