function rd_cont,root,pres=pres,logP=logP,n_e=n_e,wrange=wrange,dbdir=dbdir,$ abund=abund,metal=metal,wvl=wvl,logT=logT,fcstr=fcstr,verbose=verbose,$ chdwvl=chdwvl,twoph=twoph,noff=noff,nofb=nofb,cocofil=cocofil,$ xlogT=xlogT, _extra=e ;+ ;function rd_cont ; routine to read in the output of WRT_CONT_* and return continuum ; emissivities [1e-23 erg cm^3/s/Ang] for specified criteria ; ;syntax ; emis=rd_cont(root,pres=pres,logP=logP,n_e=n_e,wrange=wrange,$ ; dbdir=dbdir,abund=abund,metal=metal,wvl=wvl,logT=logT,fcstr=fcstr,$ ; chdwvl=chdwvl,/twoph,/noff,/nofb,cocofil=cocofil,/xlogT) ; ;parameters ; root [INPUT; required] prefix of files containing the emissivities ; * if set to 'chianti', will carry out a real-time computation ; of the continua using CHIANTI routines FREEFREE, FREEBOUND, ; and TWO_PHOTON ; * if set to 'apec', will read in the emissivities from the ; APEC continuum database table, which includes both regular ; continuum and the so-called pseudo-continuum ;keywords ; pres [INPUT; default: 1e15] pressure [cm^-3 K] ; logP [INPUT; default: 15] log(pressure [cm^-3 K]) ; * LOGP takes precedence over PRES ; * if no match is found, uses nearest available dataset ; n_e [INPUT] electron density [cm^-3]. if set, ; * overrides LOGP ; * appends a "D" to DBDIR ; wrange [INPUT; default: all] wavelength range ; * if 2-element array, wrange==[wmin,wmax] ; * if 1-element array or scalar, wrange==[wr(0),wr(0)] ; dbdir [INPUT; default: "cont"] directory to look for files in ; * $ : look for environment variable. if missing, ; following special case is hardcoded -- ; -- $CONT : see variable CONT_DB ; * ignored when ROOT='chianti' or ROOT='apec' ; abund [INPUT; default: anders & grevesse] abundances ; metal [INPUT; default: 1] metallicity ; * assumes ABUND to be of metallicity 1 ; * if set, scales abundances from ABUND appropriately ; wvl [OUTPUT] bin beginning wavelengths and final ending ; wavelength [Angstroms] ; logT [OUTPUT] log10(Temperatures [K]) ; fcstr [OUTPUT] if specified in call on input, will contain all ; the output variables ; {CONT_INT,LOGT,WVL,REH,TKEV,EKEV,MIDWVL} ; in that order in a structure. ; chdwvl [INPUT] bin size to use for ROOT='chianti' ; * default is 0.005 Angstrom ; twoph [INPUT] set this to _include_ two-photon continuum ; noff [INPUT] set this to _exclude_ free-free continuum ; nofb [INPUT] set this to _exclude_ free-bound continuum ; * TWOPH, NOFF, NOFB are ignored unless ROOT='chianti' ; cocofil [INPUT] full path to the location of the APEC ; continuum emissivities file ; * default is to use '$APED', which points to ; /data/atomdb/apec_coco.fits ; * used only if CEROOT points to APEC ; xlogT [INPUT] if set, returns the emissivity in whatever ; native logT grid as is in the databases ; * the default is to rebin it with REBINX() and put ; it on the grid logT=findgen(81)*0.05+4 ; verbose [INPUT] controls chatter ; _extra [INPUT ONLY] pass defined keywords to RD_IONEQ(), ; FREE_FREE_CH, FREE_BOUND_CH, and TWO_PHOTON_CH. ; ;restrictions ; * WRT_CONT_* must have been run first ; * cannot handle logT arrays that change from file to file ; * requires subroutines ; INICON ; GETABUND ; SYMB2ZION [LAT2ARAB] ; SETSYSVAL [LEGALVAR] ; RD_IONEQ ; FREE_FREE_CH [ITOH_CH,SUTHERLAND_CH] ; FREE_BOUND_CH ; TWO_PHOTON_CH [RD_CHIANTI] ; RD_COCO ; GETPOADEF ; from CHIANTI: READ_IONEQ, READ_GFFGU, BILINEAR, READ_IP, ; READ_KLGFB, CONCAT_DIR, GET_IEQ, FREEBOUND_ION, ZION2FILENAME, ; FILE_EXIST, DATATYPE, FILE_STAT, DATA_CHK, SINCE_VERSION, ; READ_FBLVL, IS_DIR, EXIST, GET_DELIM, OS_FAMILY, TAG_EXIST, ; IDL_RELEASE, VERNER_XS, KARZAS_XS ; ; ;history ; vinay kashyap (Dec97; modified from RD_LINE) ; allowed dbdir to be "$CONT" or something that may be read in ; via an environmental variable (VK; Dec98) ; added missing file close statement (VK; 1999Feb) ; bug was preventing reading nearest available database (VK; Mar99) ; bug in last bug correction, wasn't reading element files (VK; Apr99) ; allowed FCSTR output if set in command line (VK; AugMM) ; changes suggested by Antonio Maggio: *_DB defaults (VK; NovMM) ; added byte-order check; added calls to SWAP_ENDIAN() as suggested ; by Antonio Maggio; changed CONT_DB to CONT_CDB environment ; variable; added call to SETSYSVAL; now deciphers default ; emissivity directory location from where this file is located ; (VK; FebMMI) ; added keyword VERBOSE (VK; Aug01) ; corrected bug where first bin was being chopped off (VK; Jan04) ; check to see whether input abundances are already ratios (VK; Apr04) ; added option to use chianti routines to compute continuum, added ; TWO_PHOT keyword to toggle inclusion of two photon continuum when ; calculating via CHIANTI to save time (LL; Feb05) ; added NOFF and NOFB keywords to toggle inclusion of CHIANTI free-free ; and free-bound continuum, allowed system variable !XUVTOP to be set ; if not present (LL; Jun05) ; changed keyword TWO_PHOT to TWOPH, added keyword CHDWVL, and ; cleaned up interface and internals (VK; Jun05) ; bug fix: $CONT/cont was missing the "/" on some systems (VK; May06) ; added option to read in APEC_COCO files and added keyword XLOGT ; (VK; Apr09) ; fixed complication with abund being input as ratio (VK; Nov13) ; added call to GETPOADEF (VK; Aug15) ;- ; usage if n_elements(root) ne 1 then begin print,'Usage: fc=rd_cont(root,pres=pres,logP=logP,n_e=n_e,wrange=wrange,$' print,' dbdir=dbdir,abund=abund,metal=metal,wvl=wvl,logT=logT,$' print,' chdwvl=chdwvl,fcstr=fcstr,/twoph,/noff,/nofb,cocofil=cocofil,$ print,' /xlogT)' print,' read in (or compute) continuum spectra' return,dblarr(1,1)-1 endif ; initialize elem=[ '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'] ;elements from 1-30 ; verbosity vv=0 ivar=0 & defsysv,'!VERBOSE',exists=ivar ;if !VERBOSE exists if ivar ne 0 then setsysval,'VERBOSE',vv,/getval if n_elements(verbose) gt 0 then vv=long(verbose[0]) ; databases ;zTOPDIR=filepath('',root_dir='data',subdirectory=['fubar','SCAR']) ;zTOPDIR='/data/fubar/SCAR/' zTOPDIR=getpoadef() ivar=0 & defsysv,'!TOPDIR',exists=ivar ;if !TOPDIR exists if ivar ne 0 then setsysval,'TOPDIR',zTOPDIR,/getval else begin tmp=routine_info('RD_CONT',/source,/functions) if n_tags(tmp) ne 0 then begin scardir=tmp.PATH jvar=strpos(scardir,'pro') if jvar ge 0 then begin if strpos(scardir,'rd_cont.pro',jvar) ge 0 then $ zTOPDIR=strmid(scardir,0,jvar) ;we know where RD_CONT.PRO is endif endif if vv ge 10 then message,'Default emissivity directory is: '+zTOPDIR,/info endelse CONT_DB=getenv('CONT_CDB') if CONT_DB eq '' then $ CONT_DB=filepath('cont',root_dir=filepath('emissivity',root_dir=zTOPDIR)) ;CONT_DB=zTOPDIR+filepath('',root_dir='emissivity')+'cont' ; check keywords if keyword_set(pres) then pr=pres else pr=1e15 ;pressure if keyword_set(logP) then pr=10.^(logP) ;log10(pressure) logP=alog10(pr) case n_elements(wrange) of ;wavelength range 1: wr=[wrange(0),wrange(0)] ;pick only one 2: wr=wrange ;range else: wr=[0.,10000.] ;all possible endcase if not keyword_set(dbdir) then dbdir='cont' ;directory with files ; remove trailing '/' and/or trailing 'D' c=dbdir(0) if strmid(c,strlen(c)-1,1) eq '/' then c=strmid(c,0,strlen(c)-1) if strmid(c,strlen(c)-1,1) eq 'D' then c=strmid(c,0,strlen(c)-1) if strmid(c,0,1) eq '$' then begin if c eq '$CONT' then c=CONT_DB else begin cdb=getenv(strmid(c,1,strlen(c)-1)) ;environment variable if cdb eq '' then begin ;($WHATEVER was not set c0=strupcase(c) & c1=CONT_DB ;default if strpos(c0,'CONT',1) gt 0 then c1=CONT_DB message,c+': not defined; using '+c1,/info c=c1 endif ;CDB='') endelse endif dbasedir=c ; append 'D' if necessary if keyword_set(n_e) then dbasedir=c+'D' ;if keyword_set(n_e) then dbasedir=dbdir+'D' else dbasedir=dbdir if vv gt 0 then print,'reading from '+dbasedir defabu=getabund('anders & grevesse') if n_elements(abund) lt 30 then abnd=defabu else abnd=abund ;abundances if n_elements(metal) ne 1 then met=1. else met=metal(0) ;metallicity abnd(2:*)=abnd(2:*)*met ;include metallicity ; check C, N, O, and Fe abundances to figure out whether ; abundances already appear to be ratios if abnd(6-1) gt 0.05 or $ abnd(7-1) gt 0.01 or $ abnd(8-1) gt 0.01 or $ abnd(26-1) gt 0.005 then begin message,'input abundances appear to be ratios',/informational factor=abnd abnd=factor*defabu ;bug fix 2014-nov-13 endif else factor=abnd/defabu case root of ;{special cases 'apec': begin if not keyword_set(cocofil) then cocofil='$APED' jnk=rd_coco(cocofil,wrange=wrange,abund=abnd,$ cemis=emis,wvl=wvl,wmid=wmid,logT=logT,verbose=vv, _extra=e) nw=n_elements(wvl) if not keyword_set(xlogT) then begin oldlogT=logT & logT=findgen(81)*0.05+4. emis=(rebinx(emis,oldlogT,logT)>0) endif end 'chianti': begin ;(compute with CHIANTI xivar = 0 & defsysv,'!xuvtop',exists=xivar ;if !XUVTOP exists if xivar eq 0 then begin defsysv, '!chidir', exists=chivar if chivar eq 1 then defsysv,'!xuvtop',!chidir else $ defsysv,'!xuvtop',filepath('',root_dir=zTOPDIR,subdirectory=['CHIANTI','dbase']) endif dwvl = 0.005 ; set the default wavelength binning if keyword_set(chdwvl) then dwvl=chdwvl[0] if dwvl le 0 then dwvl=0.005 nw = (wr(1)-wr(0))/dwvl ; get number of wavelength bins wvl = findgen(nw+1)*dwvl + wr(0) ; create wavelength grid wmid=0.5*(wvl(1:*)+wvl) ivar = 0 & defsysv,'!LOGT',exists=ivar if ivar ne 0 then logt = !logt else logt=findgen(81)*0.05+4.0 nt=n_elements(logt) ; define outputs (note that actual output will be the ; transpose of these -- i.e., fltarr(nt,nw)) ffint=fltarr(nw,nt) fbint=fltarr(nw,nt) tpint=fltarr(nw,nt) ioneq = rd_ioneq(findgen(31)+1,logt,_extra=e) ; accnt for ion and z arrays exchanged btwn poa and chi in ioneq chioneq = transpose(ioneq,[0,2,1]) if keyword_set(noff) and keyword_set(nofb) and not keyword_set(twoph) then begin message, $ 'No continuum requested, please check keywords NOFF,NOFB, and TWOPH',/info return, ffint+fbint+tpint endif if not keyword_set(noff) then $ freefree_ch,10.^logt,wmid,ffint,ioneq=chioneq,ieq_logt=logt,$ abund=[abnd,0],/no_setup,_extra=e if not keyword_set(nofb) then $ freebound_ch,10.^logt,wmid,fbint,ioneq=chioneq,ieq_logt=logt,$ abund=[abnd,0],/no_setup,_extra=e if keyword_set(twoph) then begin if keyword_set(n_e) then begin two_photon_ch,10.^logt,wmid,tpint,ioneq=chioneq,ieq_logt=logt,$ abund=[abnd,0],edensity=n_e,/no_setup,_extra=e endif else begin message,'WARNING: N_E not set, TWO_PHOTON not included.',/info endelse endif emis = transpose(tpint+fbint+ffint) inicon, atom=atm atmn = findgen(30)+1 ; unit conversions 1e-40 ergs/s/sr/ang -> 1e-23 ergs/s/ang emis = 4.*!pi*emis/1d17 ;4.*!pi*emis/1d17 end ;CHIANTI) else: begin ;(read from database ; output arrays wvl=[-1.] & logT=wvl ; find the right suffixes for pressure/density fnum=logp(0) & if keyword_set(n_e) then fnum=alog10(n_e(0)) if fnum ge 0 and fnum lt 10 then fsfx='0'+string(fnum,'(f3.1)') if fnum lt 0 or fnum ge 10 then fsfx=string(fnum,'(f4.1)') ; find out if there are any data files at all nfil=n_elements(fils) & fils=findfile(dbasedir+'/'+root+'*_'+fsfx,count=nfil) if nfil lt 4 then begin ;(no files present? ; find the nearest suffix filz=findfile(dbasedir+'/'+root+'_*.?',count=kfil) if kfil eq 0 then begin message,' no spectra available for '+dbasedir+'/'+root+'_*',/info return,dblarr(1,1)-1. endif sufx=strarr(kfil) for i=0,kfil-1 do sufx(i)=strmid(filz(i),strlen(dbasedir+'/'+root+'_'),4) sufx=float(sufx) & tmp=min(abs(sufx-fnum),jj) message,'dataset with suffix '+fsfx+' not available. using '+$ string(sufx(jj),'(f4.1)')+' instead',/info fnum=sufx(jj) ;reset to nearest available if fnum ge 0 and fnum lt 10 then fsfx='0'+string(fnum,'(f3.1)') if fnum lt 0 or fnum ge 10 then fsfx=string(fnum,'(f4.1)') ; ; ok, and files with this new suffix do exist.. fils=findfile(dbasedir+'/'+root+'*_'+fsfx,count=nfil) if nfil lt 4 then message,'BUG!' endif ;NFIL<4) ; need to figure out whether to swap byte order or not if n_elements(set_byte_swap) eq 0 then begin openr,uw,dbasedir+'/'+root+'_wvl',/get_lun nw=0L & readu,uw,nw if nw ge 16777216 or nw lt 0 then set_byte_swap=1 else set_byte_swap=0 if float(strmid(!version.RELEASE,0,3)) lt 5.3 then $ set_byte_swap=(-1)*set_byte_swap close,uw & free_lun,uw ; ;Antonio Maggio points out that keyword /SWAP_ENDIAN is not valid ;for earlier OPENRs, so the *function* SWAP_ENDIAN() must be used. ;however, SWAP_ENDIAN() is very slow for large double precision ;arrays, so we will do that only for older IDL versions, hence the ;multiplication by -1 above. ; endif ; first read the wavelength array... wfil=dbasedir+'/'+root+'_wvl' if set_byte_swap gt 0 then openr,uw,wfil,/get_lun,/swap_endian else $ openr,uw,wfil,/get_lun nw=0L & readu,uw,nw if set_byte_swap lt 0 then nw=swap_endian(nw) wvl=dblarr(nw) & readu,uw,wvl if set_byte_swap lt 0 then wvl=swap_endian(wvl) close,uw & free_lun,uw ; ... and the temperatures array ... tfil=dbasedir+'/'+root+'_tmp' if set_byte_swap gt 0 then openr,ut,tfil,/get_lun,/swap_endian else $ openr,ut,tfil,/get_lun nt=0L & readu,ut,nt if set_byte_swap lt 0 then nt=swap_endian(nt) logT=fltarr(nt) & readu,ut,logT if set_byte_swap lt 0 then logT=swap_endian(logT) close,ut & free_lun,ut ; ... and the "parameter history" array ... sfil=dbasedir+'/'+root+'_src' if set_byte_swap gt 0 then openr,us,sfil,/get_lun,/swap_endian else $ openr,us,sfil,/get_lun c1='' & src=[''] while not eof(us) do begin readf,us,c1 if set_byte_swap lt 0 then c1=swap_endian(c1) & src=[src,c1] endwhile & if n_elements(src) gt 1 then src=src(1:*) close,us & free_lun,us ; ... and the base emissivities. efil=dbasedir+'/'+root+'_'+fsfx if set_byte_swap gt 0 then openr,ue,efil,/get_lun,/swap_endian else $ openr,ue,efil,/get_lun mt=nt & mw=nw & emis=dblarr(nt,nw-1) readu,ue,mt & readu,ue,mw & readu,ue,emis if set_byte_swap lt 0 then mt=swap_endian(mt) if set_byte_swap lt 0 then mw=swap_endian(mw) if set_byte_swap lt 0 then emis=swap_endian(emis) close,ue & free_lun,ue ; and finally look through the element list for i=0,n_elements(elem)-1 do begin ;{the rootZ_##.# files atm=elem(i) & efil=dbasedir+'/'+root+atm+'_'+fsfx oo=where(fils eq efil,moo) if moo gt 0 then begin ;(emissivities exist if set_byte_swap gt 0 then openr,ue,efil,/get_lun,/swap_endian else $ openr,ue,efil,/get_lun emZ=dblarr(nt,nw-1) readu,ue,mt,mw,emZ if set_byte_swap lt 0 then mt=swap_endian(mt) if set_byte_swap lt 0 then mw=swap_endian(mw) if set_byte_swap lt 0 then emZ=swap_endian(emZ) close,ue & free_lun,ue emis=emis+emZ*factor(i) ;add 'em up! oi=where(wvl gt 20 and wvl lt 50) ; or should it be emis=emis+emZ*abnd(i) ?? ;print,efil,total(emZ),factor(i) endif ;MOO>0) endfor ;I=0,N_ELEMENTS(ELEM)-1} end ;read from database) endcase ;ROOT} ; select in wavelength range ;ww=wvl(0:nw-1-1) ;use bin-beginning values ww=wvl(0:nw-1) ;use bin-beginning values oo=where(ww ge wr(0) and ww le wr(1),moo) if moo gt 0 then begin if oo[0] ne 0L then oo=[oo[0]-1L,oo] wvl=wvl([oo,oo(moo-1)+1]) & emis=emis(*,oo) endif else begin wvl=[-1.] & emis=dblarr(nt,1)-1. message,' No bins in this wavelength range',/info return,emis endelse ; get TkeV and EkeV TkeV=10.^(logT)*1.380662e-16/1.6021892e-9 ;[K]*[erg/K]/[erg/keV] EkeV=12.3985/wvl ;[Ang]->[keV] ; get mid-bin values wmid=0.5*(wvl(1:*)+wvl) if arg_present(fcstr) then begin fcstr=create_struct('CONT_INT',emis,'LOGT',logt,'WVL',wvl,$ 'TkeV',TkeV,'EkeV',EkeV,'midWVL',wmid) endif return,emis end