PRO freebound_ion, temp, wvl, int, iz, ion, ip=ip, vdata=vdata, pe=pe, $ klgfb=klgfb, noverner=noverner, kev=kev ;+ ; PROJECT: CHIANTI ; ; CHIANTI is an Atomic Database Package for Spectroscopic Diagnostics of ; Astrophysical Plasmas. It is a collaborative project involving the Naval ; Research Laboratory (USA), the University of Florence (Italy), the ; University of Cambridge and the Rutherford Appleton Laboratory (UK). ; ; NAME ; ; FREEBOUND ; ; EXPLANATION ; ; Calculates the free-bound (radiative recombination) continuum ; from a single ion. Note that the output does not contain the ion ; fraction, element abundance or differential emission measure term. ; ; INPUTS ; ; TEMP Temperature in K (can be an array). ; ; WVL Wavelength in angstroms (can be an array). If the keyword ; /keV is set, then wvl should instead be given in units of ; keV. ; ; IZ Atomic number of ion (e.g., 26 = Fe) ; ; ION Spectroscopic number of ion (e.g., 13 = XIII) ; ; OUTPUTS ; ; INT Free-bound continuum intensity. Needs to be multiplied by ; element abundance, ion fraction and DEM to obtain the final ; continuum intensity. ; ; OPTIONAL INPUTS ; ; IP The ionization potential of the ion. ; ; VDATA An array containing the Verner & Yakovlev data array. ; ; PE An array containing the PE data from READ_FBLVL ; ; KLGFB An array containing the KLGFB data from READ_FBLVL ; ; [Note: the above 3 inputs are used when calling freebound_ion from ; freebound] ; ; KEYWORDS ; ; NOVERNER If set, then the Verner & Yakovlev cross-sections will ; not be used. ; ; KEV If set, then WVL is assumed to contain energies in keV units ; rather than wavelengths in Angstroms. ; ; COMMON BLOCKS ; ; None. ; ; CALLS ; ; READ_FBLVL, ZION2FILENAME, VERNER_XS, KARZAS_XS, CONCAT_DIR, ; READ_IP, READ_KLGFB, FILE_EXIST ; ; PROGRAMMING NOTES ; ; The way I treat the exponential function in the expression for the ; emissivity may seem strange, but it saves a bit of time in the ; calculation. Basically calculating exp(E-IP+E_i/T) for each level ; was time consuming so I split it into exp(E-IP/T)*exp(E_i/T). The ; first term comes out of the for loop, while the second term is the ; exponential of a vector rather than an array. This made a ~30% ; time-saving for freebound. ; ; HISTORY ; ; Ver.1, 24-Jul-2002, Peter Young ; ; Ver.2, 26-Jul-2002, Peter Young ; Added /noverner keyword. ; ; Ver.3, 30-Jul-2002, Peter Young ; Speeded up routine by modifying treatment of exponential. ; ; Ver.4, 10-Mar-2006, Peter Young ; Added /keV keyword ; ; Ver. 5 8-Jun-2010, Ken Dere ; added check for n_params() ; ;- if n_params() lt 5 then begin print,' > freebound_ion, temp, wvl, int, iz, ion, ip=ip, vdata=vdata, pe=pe, ' print,' klgfb=klgfb, noverner=noverner, kev=kev' print,' ' return endif const=3.0992d-46 nt=n_elements(temp) nwvl=n_elements(wvl) rad=dblarr(nwvl,nt) int=rad ident_wvl=make_array(nwvl,val=1,/double) ; ; * keV units * ; IF keyword_set(kev) THEN BEGIN wvl_save=wvl efact=12.398/wvl^2 ; hc/E^2 factor to multiply continuum by efact=reverse(efact) wvl=12.398/wvl ; convert energies to angstroms wvl=reverse(wvl) ENDIF ELSE BEGIN efact=1d0 ENDELSE IF (n_elements(vdata) EQ 0) AND NOT keyword_set(noverner) THEN BEGIN vdata=dblarr(10,465) dir=concat_dir(!xuvtop,'continuum') fname=concat_dir(dir,'verner_short.txt') openr,lun,fname,/get_lun readf,lun,vdata free_lun,lun ENDIF IF n_elements(pe) EQ 0 AND n_elements(klgfb) EQ 0 THEN read_klgfb,pe,klgfb IF n_elements(ip) EQ 0 THEN BEGIN dir=concat_dir(!xuvtop,'ip') fname=concat_dir(dir,'chianti.ip') read_ip,fname,ionpot,ref ip=ionpot[iz-1,ion-1] ENDIF ewvl=1d8/wvl ; ; read .fblvl file for recombined ion ; zion2filename,iz,ion,filename ename=filename+'.fblvl' IF file_exist(ename) THEN BEGIN read_fblvl,ename,l1,conf,pqn,ll,spd,mult,ecm,ecmth,eref mult=float(mult) ind=where(ecm EQ 0. AND l1 NE 1) IF ind[0] NE -1 THEN ecm[ind]=ecmth[ind] wecm=1.e+8/(ip-ecm) necm=n_elements(ecm) ENDIF ELSE BEGIN IF keyword_set(kev) THEN wvl=wvl_save return ENDELSE ; ; if the maximum energy in WVL is not enough to ionize the ion, then don't ; bother calculating emissivity ; IF max(ewvl) LT min(ip-ecm) THEN BEGIN IF keyword_set(kev) THEN wvl=wvl_save return ENDIF ; ; read .fblvl file for recombining ion ; zion2filename,iz,ion+1,filenamer enamer=filenamer+'.fblvl' IF file_exist(enamer) THEN BEGIN read_fblvl,enamer,l1r,confr,pqnr,llr,spdr,multr,ecmr,erefr multr=float(multr) ENDIF ELSE BEGIN multr=1. ENDELSE ind=where( (wecm GE 0.) AND (max(ewvl) GE ecm) ) ni=n_elements(ind) IF ind[0] EQ -1 THEN BEGIN IF keyword_set(kev) THEN wvl=wvl_save return ENDIF FOR j=0,ni-1 DO BEGIN i=ind[j] IF (i EQ 0) AND NOT keyword_set(noverner) THEN BEGIN xs=verner_xs(iz,ion,wvl,data=vdata) ENDIF ELSE BEGIN xs=karzas_xs(wvl,pqn[i],ll[i],ip-ecm[i],pe=pe,klgfb=klgfb) ENDELSE expf=ident_wvl # exp(-1.43873 * ecm[i] / temp) rad_i=const * mult[i]/multr[0] * ((ewvl^5*xs) # (1d0/temp^(1.5))) * expf rad=rad+rad_i ENDFOR ind=where(rad NE 0.) IF ind[0] NE -1 THEN BEGIN lrad=alog(rad[ind]) f=-1.43873 * (ewvl-ip) # (1d0/temp) int[ind]=exp(lrad+f[ind])*1d40/4d0/!pi ENDIF ; ; multiply int by the hc/E^2 factor if /kev set ; IF keyword_set(kev) THEN BEGIN e_arr=efact#(dblarr(nt)+1d0) int=int*e_arr ; ; reverse the wavelength dimension to match energy ordering rather than ; wavelength ordering int=int[reverse(indgen(nwvl)),*] ; ; and set wavelengths back to energy units wvl=wvl_save ENDIF END