;+
;
; PROJECT:  CHIANTI
;
;       CHIANTI is an atomic database package for the calculation of
;       continuum and emission line spectra from astrophysical plasmas. It is a 
;       collaborative project involving the Naval Research Laboratory
;       (Washington DC, USA), the Arcetri Observatory (Firenze, Italy), and the
;       Cambridge University (United Kingdom).
;
;
; NAME:
;	ASCII_WVL_DEM
;
; PURPOSE:
;
;	create an ascii file of predicted spectral line intensities and
;       wavelengths corresponding to a selected abundance and differential
;       emission measure (DEM).  derived from latex_wvl_dem
;
; CATEGORY:
;
;	science.
;
; CALLING SEQUENCE:
;
;       ASCII_WVL_DEM, Wmin, Wmax, Pressure
;
;
; INPUTS:
;
;	Wmin:   lower limit of the wavelength range of interest (Angstroms)	
;	Wmax:   upper limit of the wavelength range of interest (Angstroms)
;       Pressure:  electron density (cm-3) * temperature (K), assumed constant	
;
; OPTIONAL INPUTS:
;
;	Elemental abundance file:  selected with widget
;       Ionization equilibrium file:  selected with widget
;       Differential emission measure file:  selected with widget
;	
; KEYWORD PARAMETERS:
;
;	MINI:	Minimum intensity for line to be included in output
;
;	SNGL_ION:  specified a single ion to be used instead of the complete
;                  set of ions specified in !xuvtop/masterlist/masterlist.ions
;
;
;       PHOTONS:  units will be in photons rather than ergs
;
;       ALL:  all lines, including those with only theoretical wavelengths are printed
;
;       MASTERLIST:  string of a specific masterlist file or set to '' to select a 
;                particular masterlist file.  Otherwise, masterlist.ions is used
;
;       OUTFILE:  the name of the file to be written.
;
; OUTPUTS:
;
;	an ascii file:   linelist.lis  in the working directory by default
;
;
;
; COMMON BLOCKS:
;
;	common elvlc,l1a,term,conf,ss,ll,jj,ecm,eryd,ecmth,erydth,eref
;       common wgfa, wvl,gf,a_value
;       common upsilon,t_type,deu,c_ups,splups
;       common elements,abund,abund_ref,ioneq,ioneq_t,ioneq_ref,ionpot,ip_ref
;
;
; SIDE EFFECTS:
;
;	Some of the ions have 100+ energy levels and the inversion of the
;       necessary matrices can be time consuming
;
;
; EXAMPLE:
;
;             > ascii_wvl_dem,100.,200.,1.e+15,mini=10.,sngl_ion='c_2'
;
; MODIFICATION HISTORY:
; 	Written by:	Ken Dere
;	September 1996:     Version 1.0
;       September 1996:  mods to work with VMS
;       June 1998:  include keyword /photons
;       September 1999:  K. Dere, modified for Version 3.0
;       14-Jul-2000      Peter Young, now calls pop_solver
;       30-Oct-2000: V.6, Giulio Del Zanna, added the keyword outfile and changed
;                    default. Now it creates an 'linelist.lis' file in the
;                    working directory. Added minor things.
;
; VERSION     : Version 6, 30-Oct-2000
;
;
;-
pro ascii_wvl_dem,wmin,wmax,pressure,$
    minI=minI,photons=photons,all=all,$
    sngl_ion=sngl_ion,masterlist=masterlist, $
    outfile=outfile
;
;  wmin=short wavelength limit in Angstroms
;  wmax=long wavelength limit in Angstroms
;  pressure = electron pressure (Ne*T) in cm^(-3) K
;
;  set mini= minimum intensity for which lines will be output, i.e. mini=0.1
;  set sngl_ion to select only on ion, i.e. sngl_ion='fe_15' to do Fe_XV
;
;
common elvlc,l1a,term,conf,ss,ll,jj,ecm,eryd,ecmth,erydth,eref
common wgfa, wvl,gf,a_value
common upsilon,t_type,deu,c_ups,splups
common elements,abund,abund_ref,ioneq,ioneq_t,ioneq_ref
;
;
if n_params(0) lt 3 then begin
  print,'  IDL>  ascii_wvl_dem,wmin,wmax,pressure,$ '
  print,'     [outfile= , mini= , sngl_ion=, /photons, /all, masterlist=]'
  return
endif
;
if keyword_set(minI) then begin
   min_int=minI
endif else min_int=0.
;
setup_elements
;
n_ioneq_t=n_elements(ioneq_t)
dlnt=alog(10.^(ioneq_t(1)-ioneq_t(0)))
;
;
;  select differential emission measure (.dem) file
;
dir=add_subdir(!xuvtop,'dem')
dem_name=pickfile(path=dir,filter='*.dem',title='Select DEM File')
;   print,' latex_title= ',latex_title
;
read_dem,dem_name,dem_t,dem,dem_ref
;
i1=strpos(dem_name,'/dem/')
i2=strpos(dem_name,'.dem')
latex_title=strmid(dem_name,i1+5,i2-i1-5)
;
;
gdt=where((ioneq_t ge min(dem_t)) and (ioneq_t le max(dem_t)))
;
dem_int1=10.^(spline(dem_t,dem,ioneq_t(gdt)))
dem_int=fltarr(n_ioneq_t)
ngt=n_elements(gdt)
for igt=0,ngt-1 do begin
   dem_int(gdt(igt))=dem_int1(igt)
endfor
;
b=where((ioneq_t lt min(dem_t)) and (ioneq_t gt max(dem_t)))
if max(b) ge 0 then dem_int(b)=0.
;  for it=0,n_ioneq_t-1 do begin
;     print,it,ioneq_t(it),dem_int(it)
;  endfor

;
;
spd=['S','P','D','F','G','H','I','K']
jvalue=['0','1/2','1','3/2','2','5/2','3','7/2','4','9/2','5','11/2','13/2','15/2','17/2','19/2']
;
maxwvl=10000
;
mast_z=intarr(maxwvl)
mast_ion=intarr(maxwvl)
mast_lvl=intarr(maxwvl)
mast_wvl=fltarr(maxwvl)
mast_A=fltarr(maxwvl)
;
mast_term1=strarr(maxwvl)
mast_desig1=strarr(maxwvl)
mast_term2=strarr(maxwvl)
mast_desig2=strarr(maxwvl)
;
mast_abund=fltarr(maxwvl)
mast_ioneq=fltarr(maxwvl)
mast_tmax=fltarr(maxwvl)
mast_pop=fltarr(maxwvl)
mast_popazidem=fltarr(maxwvl)
;
mast_index=0
;
;
;
;  open the file that has the names of the ions
;
if ch_datatype(masterlist,2) gt 0 then begin
   mname=masterlist
endif else begin
   mname=!xuvtop+'/masterlist/masterlist.ions'
endelse
;
;
gname=''
elstage=''
;
;   main input and calculation loop  **************
;
;
;   read the name of the ions unless a single ion (sngl_ion) has been specified
;
if (not keyword_set(sngl_ion)) then begin
   read_masterlist,mname,list
;
endif else list=sngl_ion
;
nlist=n_elements(list)
;
for ilist=0,nlist-1 do begin
;
   gname=list(ilist)
   convertname,gname,iz,ion
   ion2spectroscopic,gname,snote
;
print,' current ion = ',gname  
;
abundtest=abund(iz-1)
if (abundtest gt 0.) then begin  ;********************************************
;
   elstage=elstage+snote+', '  
;
;
setup_ion,gname,wmin,wmax,wvltst,lvl1,lvl2,wvl1,gf1,a_value1
;
;
if (wvltst gt 0) then begin  ;***********************************************
;
locname=strlowcase(gname)
pos=strpos(locname,'_')
l=strlen(pos)
first=strmid(locname,0,pos)
last=strmid(locname,pos+1,l-pos-1)
;
if strpos(last,'d') ge 0 then dielectronic=1 else dielectronic=0
;
;
;
  ntrans=n_elements(lvl1)
  nlvls=max([lvl1,lvl2])

;
;  calculate level populations
;
        this_ioneq=ioneq(*,iz-1,ion-1+dielectronic)
;
        itmax=where(this_ioneq eq max(this_ioneq))
        ltmax=ioneq_t(itmax)
        tmax=10.^(ioneq_t(itmax))
        tmax=tmax(0)  ; necessary to make tmax a scalar
;
       gioneq_t=where((this_ioneq gt 0.) and (dem_int gt 0.),ngt)
;
;
        popazidem=dblarr(nlvls)
;
if ngt gt 0 then begin
;
for ig=0,ngt-1 do begin
   it=gioneq_t(ig)
   tt=10.^ioneq_t(it)
   edensity=pressure/tt
   pop_solver,tt,edensity,pop
   pop = reform(pop)
   dt=tt*dlnt
   popazidem(0)=popazidem(*)+pop(*)*this_ioneq(it)*dem_int(it)*dt/edensity
endfor
;


   for itrans=0,ntrans-1 do begin
      l1=lvl1(itrans)-1
      l2=lvl2(itrans)-1
      ww=wvl1(itrans)
;     print,ww,l1,l2
;
;     can the level be excited?
      upssiz=size(splups)
      upssiz2=upssiz(2)-1
      if(l2 le upssiz2) then begin
         maxups=max(splups(*,l2,0))
         if(maxups gt 0.) then upstst=1 else upstst=0
      endif else upstst=0
;
;
;
     if (abs(ww) gt wmin) and (abs(ww) le wmax) and (upstst gt 0) then begin
        if(abs(ww) gt 0.) then begin   ;  select only observed wvl's
           mast_z(mast_index)=iz
           mast_ion(mast_index)=ion
           mast_lvl(mast_index)=lvl2(itrans)
           mast_wvl(mast_index)=wvl1(itrans)
           mast_A(mast_index)=a_value1(itrans)
           mast_abund(mast_index)=abund(iz-1)
;
           mast_tmax(mast_index)=ltmax
           mast_ioneq(mast_index)=this_ioneq(itmax)
           mast_popazidem(mast_index)=popazidem(l2)
;
;      print,iz,ion,wvl1(itrans),popazidem(l2)
       
;
;  get lower level designation
           term0=strtrim(term(l1),2)
;
           term1=''
           blank=strpos(term0,' ')
           while blank gt 0 do begin
              term1=term1+' '+strmid(term0,0,blank)
              term0=strmid(term0,blank,100)
              term0=strtrim(term0,2)
              blank=strpos(term0,' ')
           endwhile
           ch_remchar,term1,'^'
           mast_term1(mast_index)=strlowcase(term1)
;
           jinteger=fix(2.*jj(l1))
           jstring=jvalue(jinteger)
;
           spins=strtrim(string(ss(l1),'(i2)'),2)
           mast_desig1(mast_index)= $
              spins+spd(ll(l1))+jstring
;
;
;  get upper level designation
           term0=strtrim(term(l2),2)
;
           term2=''
           blank=strpos(term0,' ')
;           print,term0
           while blank gt 0 do begin
              term2=term2+' '+strmid(term0,0,blank)
              term0=strmid(term0,blank,100)
              term0=strtrim(term0,2)
              blank=strpos(term0,' ')
           endwhile
           ch_remchar,term2,'^'
           mast_term2(mast_index)=strlowcase(term2)

;
           jinteger=fix(2.*jj(l2))
           jstring=jvalue(jinteger)
;
           spins=strtrim(string(ss(l2),'(i2)'),2)
           mast_desig2(mast_index)= $
              spins+spd(ll(l2))+jstring
;
           mast_index=mast_index+1
           if mast_index gt  maxwvl-1 then begin
              print,' need to increase size of maxwvl'
              print,' current value of maxwvl = ',maxwvl
              stop
           endif
        endif
     endif
   endfor ;  over itrans  (spectral lines)
;
endif  else begin
   print,' ionization equilibrium eq 0 in DEM temperature interval' ;  ngt gt 0
endelse
;
endif else begin
   print,' no lines in specified wavelength interval'    ; wvltest
endelse
;
endif else print,' elemental abundance = 0'   ;  abundtest  ****************
;
;
;
endfor   ;   loop over ilist  ************************
;
;   end of main input and calculation loop   *****************************
;
;
print,' number of lines =',mast_index
;
if mast_index eq 0 then begin
    print,' no lines found with the specifications listed'
    return
endif
;
if mast_index gt maxwvl then print,' must increase size of maxwvl to ',mast_index
;
;
;
mast_z=mast_z(0:mast_index-1)
mast_ion=mast_ion(0:mast_index-1)
mast_lvl=mast_lvl(0:mast_index-1)
mast_wvl=mast_wvl(0:mast_index-1)
mast_A=mast_A(0:mast_index-1)
mast_term1=mast_term1(0:mast_index-1)
mast_desig1=mast_desig1(0:mast_index-1)
mast_term2=mast_term2(0:mast_index-1)
mast_desig2=mast_desig2(0:mast_index-1)
;
mast_abund=mast_abund(0:mast_index-1)
mast_tmax=mast_tmax(0:mast_index-1)
mast_ioneq=mast_ioneq(0:mast_index-1)
mast_popazidem=mast_popazidem(0:mast_index-1)

;

srt_index=sort(abs(mast_wvl))
mast_z(0)=mast_z(srt_index)
mast_ion(0)=mast_ion(srt_index)
mast_lvl(0)=mast_lvl(srt_index)
mast_wvl(0)=mast_wvl(srt_index)
mast_A=mast_A(srt_index)
mast_term1=mast_term1(srt_index)
mast_desig1=mast_desig1(srt_index)
mast_term2=mast_term2(srt_index)
mast_desig2=mast_desig2(srt_index)
;
mast_abund=mast_abund(srt_index)
mast_tmax=mast_tmax(srt_index)
mast_ioneq=mast_ioneq(srt_index)
mast_popazidem=mast_popazidem(srt_index)
;
nlines=mast_index
;
;for i=0,nlines-1 do begin
;   print,mast_z(i),mast_ion(i),mast_lvl(i),mast_wvl(i),mast_A(i),mast_tmax(i),  $
;     format='$(3i3,f12.3,2e10.3)'
;endfor
;
lines=0
;

IF n_elements(outfile) EQ 0 THEN outfile ='linelist.lis'
openw,luo,outfile,/get_lun
;
;
printf,luo,'CHIANTI:  Predicted XUV Line Intensities'
printf,luo,' '
printf,luo,strupcase(latex_title)
printf,luo,'Pressure (NeT) = '+string(pressure,'(e10.2)')
printf,luo,'Wavelength range:  '+strtrim(string(wmin,'(i4)'),2)+' to '+strtrim(string(wmax,'(i4)'),2)
printf,luo,'Number of lines:  '+strtrim(string(mast_index,'(i4)'),2)
;
;
;
printf,luo,' Ions included:'
printf,luo,elstage
;
;
;
;
printf,luo,' Abundances:'
;
;
;
asize=size(abund_ref)
for iref=0,asize(1)-1 do begin
   if (strpos(abund_ref(iref),'%elemental') ge 0) then begin
      printf,luo,abund_ref(iref)
   endif
endfor
;
;
;
printf,luo,' Differential Emission Measure:'
;
asize=size(dem_ref)
  for iref=0,asize(1)-1 do begin
   if (strpos(dem_ref(iref),'%dem') ge 0) then begin
      printf,luo,dem_ref(iref)
   endif
  endfor
;
;
;
for i=0,nlines-1 do begin

zion2spectroscopic,mast_z(i),mast_ion(i),elstage
;
if mast_wvl(i) gt 0. then begin
   wvls=strtrim(string(abs(mast_wvl(i)),'(f12.4)'),2)
endif else wvls='*'+strtrim(string(abs(mast_wvl(i)),'(f12.3)'),2)
;
term1=strtrim(mast_term1(i),2)
desig1=strtrim(mast_desig1(i),2)
term2=strtrim(mast_term2(i),2)
desig2=strtrim(mast_desig2(i),2)
;
;   az=string(mast_abund(i),'(e9.2)')
    item1=string(10.^mast_tmax(i),'(e8.1)')
;   item2=string(mast_ioneq(i),'(e9.2)')
;   item2=string(mast_pop(i),'(e9.2)')
;
if keyword_set(photons) then begin
    hc=1.d/(4.d*3.14159d)
    format='(a10,a10,a25,a3,a25,e9.1,e10.2)'
;    format='(a10,f10.3,a25,a3,a25,e9.1,e10.2)'
endif else begin
    hc=6.626d-27*2.998d+10*1.d+8/(4.d*3.14159d*abs(mast_wvl(i)))
    format='(a10,a10,a25,a3,a25,e9.1,f10.1)'
;    format='(a10,f10.3,a25,a3,a25,e9.1,f10.1)'
endelse
;
    intensity=hc*mast_abund(i)*mast_A(i)*mast_popazidem(i)
    item2=string(intensity,'(e12.2)')
;
intest=intensity ge min_int
if keyword_set(all) then wvltest=1 else wvltest=mast_wvl(i) gt 0.
;
if intest and wvltest then begin
;     printf,luo,elstage+'  '+wvls+'  ' + term1+' '+desig1+' - ' + term2+' '+desig2 +' '+item1+' ' +item2
;
level1=term1+' '+desig1
level2=term2+' '+desig2
level2=ch_strpad(level2,25,/after)
elstage=ch_strpad(elstage,10,/after)
printf,luo,elstage,wvls,level1,' - ',level2,10.^mast_tmax(i),intensity, $
    format=format
;
endif   ; int>min_int
;
endfor   ;  iline
;
;

free_lun,luo


print,' The file '+outfile+ ' has been created !'

end
;