;+
;
;
; 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:
;	TEMPERATURE_RATIOS
;
; PURPOSE:
;
;	calculate and display temperature sensitivity of line intensity ratios
;
; CATEGORY:
;
;	spectral diagnostics
;
; PROCEDURE :
;
;       The intensities (Population of the upper level * A)  of the lines within
;       the selected ion are first  calculated, either at constant pressure or
;       at constant density (however specified in the input). They are plotted
;       in window 0. The intensities  relative to the brightest  reference line
;       are then plotted in window 1. A widget allows the user to select a
;       number of lines (at least one!) for the numerator of the ratio, and a
;       number of lines for the denominator. In case of multiple selections, the
;       line intensities are summed. The ratio values are plotted in window 2, and
;       optionally also saved in a text file. A postscript file can also be
;       created. The ratio values, calculated at twice and half the prescribed
;       density  are also calculated and overplotted, to show how the
;       temperature ratio also depends on the density.
;
; CALLING SEQUENCE:
;
;
;      > temperature_ratios,ion,wmin,wmax,Log10(tempmin),Log10(tempmax),temperature,ratio,description,$
;       [pressure= ,density= , psfile= , outfile= ]
;
;
; EXAMPLE:
;
;      >temperature_ratios,'o_5',1000.,1500.,4.,6.,temperature,ratio,desc,density=1.e+12,$ 
;       psfile='test.ps', outfile='test.txt'
;
;       then select  the 1371/1218 line ratio 
;
;      > print,desc
;      CHIANTI:   O V  1371.2920/1218.3900 
;
;
; INPUTS:
;
;       Ion:   the CHIANTI style name of the ion, i.e., 'c_5' for C V
;       Wmin:  minimum  wavelength  limit in Angstroms
;       Wmax:  maximum  wavelength  limit in Angstroms
;       Tempmin:  log10 of lowest temperature of interest, i.e. 4 for 10.^4 K
;       Tempmax:  log10 of highest temperature of interest
;
; OPTIONAL INPUTS:
;
;	Must specify indices of lines which are to form the ratio
;	
; OUTPUTS:
;
;	Temperature:  an array of temperatures spanning Tempmin to Tempmax
;       Ratio:  an array of the intensity ratio of the selected lines
;       Desc:  a short string description of the selected line ratio
;
;
; OPTIONAL OUTPUTS:
;
;	Ps and/or  text file with the  intensity ratio.
;
;
; KEYWORD PARAMETERS:
;
;	PRESSURE:  calculates the intensity ratios for constant pressure,
;                 i.e. calculates the intensities of the lines at each
;                 temperature T and density = pressure / T.
;		
;	DENSITY:  calculates the intensity ratios for constant density.
;
;                  If neither density or pressure are specified, a constant
;                  density of 1.e+10 cm^-3 is assumed as default.
;
;       OUTFILE:  the (optional) name of the output ascii file where a 
;                   listing of the line ratio intensities as a function of 
;                   temperature is saved.
;
;       PSFILE:  the (optional) name of the output postscript file 
;                  where a plot of the chosen temperature sensitive line
;                  ratio is saved.
;
; CALLS: 
;
;        read_ioneq, convertname, ion2spectroscopic,ion2filename, read_wgfa2,
;        read_elvlc, read_splups, pop_solver, ch_xselect_s
;
; COMMON:
;
;       elvlc,l1a,term,conf,ss,ll,jj,ecm,eryd,ecmth,erydth,eref
;       wgfa, wvl,gf,a_value
;       upsilon,t_type,deu,c_ups,splups
;
; RESTRICTIONS:
;
;               
; SIDE EFFECTS: None known yet.
;                 
;
; MODIFICATION HISTORY:
; 	Written by:	Ken Dere
;	May 1996:     Version 2.0, Ken Dere
;       April 2000:   V. 3 Ken Dere modified for V3
;       14-Jul-2000   V. 4  Peter Young, now calls pop_solver
;        2-Oct-2000   V. 5 Giulio Del Zanna, corrected an error in the creation of
;        the string list of the lines in the ratio. Also corrected a few minor errors.
;        Removed the device,window_state call, and added a few other 'cosmetic'
;        adjustements.  
;	Version 6, 21-Dec-2000, William Thompson
;		Modified for better cross-platform capability.
;
; VERSION     : 
;       Version 6, 21-Dec-2000
;
;-
pro temperature_ratios,ions,wmin,wmax,tempmin,tempmax,$
               temperature,ratio,description,$
               pressure=pressure,density=density,psfile=psfile,outfile=outfile
;
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
;
if n_params(0) lt 8 then begin
   print,' IDL> temperature_ratios,ion,wmin,wmax,Log10(tempmin),Log10(tempmax),$ '
   print,' IDL>  temperature,ratio,description,$ '
   print,' IDL>      [pressure= ,density= , psfile= , outfile= ]'
   print,' '
   print,' i.e.: '
   print,' IDL>  temperature_ratios, "c_5" ,40.,50.,5.,7.,temp,rat,desc '
   return
endif
;
default_density=1.e+10
print, ''
if (keyword_set(pressure) ne 1) and (keyword_set(density) ne 1) then begin
   print,' using constant default density = ',default_density
endif

;
charsiz=1.5
xtitle='Temperature (K)'
;
ioneq_name=!xuvtop+'/ioneq/'+!ioneq_file
  print,' using '+!ioneq_file+ ' as default ionization equilibrium file'

read_ioneq,ioneq_name,ioneq_t,ioneq,ioneq_ref
dlnt=alog(10.^(ioneq_t(1)-ioneq_t(0)))
;
convertname,ions,iz,ion
this_ioneq=ioneq(*,iz-1,ion-1)
hit=where(this_ioneq eq max(this_ioneq))
tmax=ioneq_t(hit)
tmax=10.^(tmax(0))

print,'Maximum Temperature (ionization equilibrium) = ',tmax,format='(a45,e10.3)'
;
;
dlogtemp=0.10                   ;  increment in logT for calculating ratios
;
ntemp=fix((tempmax-tempmin)/dlogtemp)+1
;
;
;
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']
;
;
ion2spectroscopic,ions,species
;
;
ion2filename,ions,fname
wname=fname+'.wgfa'
elvlcname=fname+'.elvlc'
upsname=fname+'.splups'
;
;
;  read in level information, wavelengths, gf and A values from .wgfa files
;
read_wgfa2,wname,lvl1,lvl2,wvl1,gf1,a_value1,wgfaref ;
;
ntrans=n_elements(lvl1)
nlvls=max([lvl1,lvl2])
wvl=fltarr(nlvls,nlvls)
gf=fltarr(nlvls,nlvls)
a_value=fltarr(nlvls,nlvls)
for itrans=0,ntrans-1 do begin
   wl1=lvl1(itrans)
   wl2=lvl2(itrans)
   wvl(wl1-1,wl2-1)=wvl1(itrans)
   gf(wl1-1,wl2-1)=gf1(itrans)
   a_value(wl1-1,wl2-1)=a_value(wl1-1,wl2-1)+a_value1(itrans)
endfor
;
;
read_elvlc,elvlcname,l1a,term,conf,ss,lla,jj,ecm,eryd,ecmth,erydth,eref
mult=2.*jj+1.
;
g=where(ecm eq 0.)
;
if(max(g) gt 0) then begin
   ecm(g)=ecmth(g)
endif
;
;
read_splups,upsname,t_type,gfu,deu,c_ups,splups,upsref
;
;
;
;
list=intarr(ntrans)
list_wvl=fltarr(ntrans)
list_a=fltarr(ntrans)
list_lvl2=intarr(ntrans)
list_descr1=strarr(ntrans)
list_descr2=strarr(ntrans,ntemp)
;
nlist=0
;
;
;   find out the number of excited levels
;
pop_solver,tmax,10.^tempmin,pop
pop = reform(pop)
popsiz=size(pop)
popsiz1=popsiz(1)-1
;
;   sort out list of lines
;
for itrans=0,ntrans-1 do begin
   l1=lvl1(itrans)-1            ; -1 for IDL indexing
   l2=lvl2(itrans)-1
   ww=wvl1(itrans)
;
;  can the level be excited?
   if((l2) le popsiz1) then poptst=1 else poptst=0
;
;
   if (ww gt wmin) and (ww le wmax) and (poptst gt 0) then begin
      list(nlist)=itrans
      list_wvl(nlist)=wvl1(itrans)
      list_a(nlist)=a_value(l1,l2)
      list_lvl2(nlist)=l2
;
;
;  get lower level designation
      term0=strtrim(term(l1),2)
;
      term2=''
      blank=strpos(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
      term1=strlowcase(term2)
;
      jinteger=fix(2.*jj(l1))
      jstring=jvalue(jinteger)
;
      spins=strtrim(string(ss(l1),'(i2)'),2)
      list_descr1(nlist)=term1+' '+spins+spd(lla(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
      term2=strlowcase(term2)
;
      jinteger=fix(2.*jj(l2))
      jstring=jvalue(jinteger)
;
      spins=strtrim(string(ss(l2),'(i2)'),2)
      list_descr2(nlist)=term2+' '+spins+spd(lla(l2))+jstring
;
;
      nlist=nlist+1
;
   endif                        ; wvl within wmin, wmax
;
endfor                          ;  loop over itrans
;
;
;
;   list=list(0:nlist-1)
list_wvl=list_wvl(0:nlist-1)
list_a=list_a(0:nlist-1)
list_lvl2=list_lvl2(0:nlist-1)
list_descr1=list_descr1(0:nlist-1)
list_descr2=list_descr2(0:nlist-1)
;
;
;
srt_wvl=sort(list_wvl)
;
;   list=list(srt_wvl)
list_wvl=list_wvl(srt_wvl)
list_a=list_a(srt_wvl)
list_lvl2=list_lvl2(srt_wvl)
list_descr1=list_descr1(srt_wvl)
list_descr2=list_descr2(srt_wvl)
;
list_int=fltarr(nlist,ntemp)
temperature=fltarr(ntemp)
;
;  for sorted list, calculated level populations and intensities
;
for itemp=0,ntemp-1 do begin    ; loop over densities
;
   temperature(itemp)=10.^(tempmin+itemp*dlogtemp)
;
   if keyword_set(pressure) then begin
      edensity=pressure/temperature(itemp)
   endif else if keyword_set(density) then begin
      edensity=density
   endif else edensity=default_density
;
   pop_solver,temperature[itemp],edensity,pop
   pop = reform(pop)
;
   for ilist=0,nlist-1 do begin
      list_int(ilist,itemp)=pop(list_lvl2(ilist))*list_a(ilist)
   endfor
;
endfor                          ;   loop over densities
;
ytitle='Relative Intensity (Population of upper level)*A'
;
intmax=max(list_int)
intmin=intmax/1000.
;
;device,window_state=ws
;if(ws(0) ne 1) then  window,0,ysize=425,xsize=525 else wset,0

window,0,ysize=425,xsize=525 

;
plot_oo,fltarr(2),xr=[temperature(0),temperature(ntemp-1)],yr=[intmin,intmax], $
  xtitle=xtitle,ytitle=ytitle

;
;
for ilist=0,nlist-1 do begin
   oplot,temperature,list_int(ilist,*)
endfor
;
rand1=(ntemp-1)*randomu(seed,nlist) ; to try to help with problem of labeling lines
rand2=(ntemp-1)*randomu(seed,nlist)
;
;
for ilist=0,nlist-1 do begin
   itemp=rand1(ilist)
   xyouts,temperature(itemp),list_int(ilist,itemp),string(ilist,'(i3)'),size=charsiz,align=0.5,noclip=0
   xyouts,temperature(1),list_int(ilist,1),string(ilist,'(i3)'),charsiz=charsiz ,align=0.5,noclip=0
   xyouts,temperature(ntemp-2),list_int(ilist,ntemp-2),string(ilist,'(i3)'),charsiz=charsiz,align=0.5,noclip=0
endfor
;
;
;  get a reference line
;
imax=where(list_int eq max(list_int))
ref=imax-fix(imax/nlist)*nlist
oplot,temperature,list_int(ref,*),thick=2
;
refstr=strtrim(string(ref,'(i3)'),2)
wvlstr=strtrim(string(list_wvl(ref),'(f10.4)'),2)
ytitle='Intensity relative to line '+refstr+' at '+wvlstr
;
;device,window_state=ws
;if(ws(2) ne 1) then  window,2,ysize=425,xsize=525 else wset,2

window,1,ysize=425,xsize=525 

;
;  plot ratio to reference line
;
plot_oo,fltarr(2),xr=[temperature(0),temperature(ntemp-1)],yr=[1.e-3,1.e+1], $
  xtitle=xtitle,ytitle=ytitle

for ilist=0,nlist-1 do begin
   oplot,temperature,list_int(ilist,*)/list_int(ref,*)
endfor
;
;
for ilist=0,nlist-1 do begin
   itemp=rand2(ilist)
   xyouts,temperature(itemp),list_int(ilist,itemp)/list_int(ref,itemp),string(ilist,'(i3)'),charsiz=charsiz,align=0.5,noclip=0
   xyouts,temperature(1),list_int(ilist,1)/list_int(ref,1),string(ilist,'(i3)'),charsiz=charsiz,align=0.5,noclip=0
   xyouts,temperature(ntemp-2),list_int(ilist,ntemp-2)/list_int(ref,ntemp-2),string(ilist,'(i3)'),charsiz=charsiz,align=0.5,noclip=0
endfor
;

print, ''
print,' id # lambda (A) intensity         transition'


print,' ------------------------------------------------------'
;
dash=' - '
;
;
options=strarr(nlist)
nstatus=intarr(nlist)
dstatus=intarr(nlist)
;
for i=0,nlist-1 do begin
   print,i,list_wvl(i),list_int(i,ntemp-1),list_descr1(i),dash,ch_strpad(list_descr2(i),15,/after), $
     format='(i4,f10.4,e10.2,a15,a3,a15)'
   options(i)=string(i,list_wvl(i),list_int(i,ntemp-1),list_descr1(i), $
                     dash,ch_strpad(list_descr2(i),15,/after), $
                     format='(i4,f10.4,e10.2,a15,a3,a15)')
endfor
;
;
ch_xselect_s,options,nstatus,abort,title='Select lines for numerator'
;
ch_xselect_s,options,dstatus,abort,title='Select lines for denominator'
;
numerator=fltarr(ntemp)
for i=0,nlist-1 do begin
   if nstatus(i) eq 1 then begin
      numerator=numerator+list_int(i,*)/list_wvl(i)
   endif
endfor
;
denominator=fltarr(ntemp)
for i=0,nlist-1 do begin
   if dstatus(i) eq 1 then begin
      denominator=denominator+list_int(i,*)/list_wvl(i)
   endif
endfor
;
ratio=numerator/denominator

;
;
description='CHIANTI:   '+species+'  '
;
i=where(nstatus eq 1,nnum)

IF  nnum EQ  1 THEN  BEGIN 
   description=description+strtrim(string(list_wvl(i(0)),'(f10.4)'),2)+'/'
ENDIF  ELSE   BEGIN
   description=description+' ('+strtrim(string(list_wvl(i(0)),'(f10.4)'),2)
   for j=1,nnum-1 DO  description=description+'+'+strtrim(string(list_wvl(i(j)),'(f10.4)'),2)
   description=description+')/'
END 
;
;
i=where(dstatus eq 1,dnum)
IF  dnum EQ  1 THEN  BEGIN
   description=description+strtrim(string(list_wvl(i(0)),'(f10.4)'),2)
ENDIF  ELSE   BEGIN
   description=description+'('+strtrim(string(list_wvl(i(0)),'(f10.4)'),2)
   for j=1,dnum-1 DO  description=description+'+'+strtrim(string(list_wvl(i(j)),'(f10.4)'),2)
   description=description+')'
END 

;
print,'-------------------------------------------------------'

;read,line1,line2, $
;   prompt='type indices (0:'+strtrim(string(nlist-1,'(i3)'),2)+') for level1 (numerator), level2 ;(denominator):  '
;
;ratio=list_int(line1,*)*list_wvl(line2)/(list_int(line2,*)*list_wvl(line1))
;ratio=reform(ratio)
;
;description=species+'  ;'+string(list_wvl(line1),'(f10.3)')+'/'+strtrim(string(list_wvl(line2),'(f10.3)'),2)
;

;device,window_state=ws
;if(ws(1) ne 1) then  window,1,ysize=425,xsize=525 else wset,1

window,2,ysize=425,xsize=525

ytitle='Intensity (Ergs) Ratio'
plot_oi,temperature,ratio,xr=[temperature(0),temperature(ntemp-1)],yr=[min(ratio),max(ratio)],title=description, $
  xtitle=xtitle,ytitle=ytitle , xstyle=1, ystyle=1 

;
if keyword_set(pressure) then begin
   xyouts,temperature(1),max(ratio)/2.,'Constant Pressure = '+string(pressure,'(e10.2)')+' (cm!e-3!n !eo!nK)',$
     charsiz=charsiz
endif else if keyword_set(density) then begin
   xyouts,temperature(1),max(ratio)/2.,'Constant Density = '+string(edensity,'(e10.2)')+' (cm!e-3!n)',$
     charsiz=charsiz
endif else begin
   xyouts,temperature(1),max(ratio)/2.,'Constant Density = '+string(edensity,'(e10.2)')+' (cm!e-3!n)',$
     charsiz=charsiz
endelse
;


;
if keyword_set(outfile) then begin
   openw,luo,outfile,/get_lun
;
   printf,luo, outfile
   printf,luo,'  intensity ratios (ergs) calculated by CHIANTI'
   printf,luo,description
   if keyword_set(pressure) then begin
      printf,luo,'Constant Pressure = '+string(pressure,'(e10.2)')
   endif else if keyword_set(density) then begin
      printf,luo,'Constant Density = '+string(edensity,'(e10.2)')
   endif else begin
      printf,luo,'Constant Density = '+string(edensity,'(e10.2)')
   endelse
   printf,luo, ' '
;
   for i=0,ntemp-1 do begin
      printf,luo, temperature(i),ratio(i)
   endfor
;
   close,luo
;
ENDIF 
;
;
;  also include effects of factor of 2 changes in pressure
;
for itemp=0,ntemp-1 do begin    ; loop over densities
;
   temperature(itemp)=10.^(tempmin+itemp*dlogtemp)
;
   if keyword_set(pressure) then begin
      edensity=pressure/temperature(itemp)
   endif else if keyword_set(density) then begin
      edensity=density
   endif else edensity=default_density
;
   edensity=2.*edensity

;

   pop_solver,temperature[itemp],edensity,pop
   pop = reform(pop)
;
   for ilist=0,nlist-1 do begin
      list_int(ilist,itemp)=pop(list_lvl2(ilist))*list_a(ilist)
   endfor
;
endfor                          ;   loop over temperatures
;
numerator=fltarr(ntemp)
for i=0,nlist-1 do begin
   if nstatus(i) eq 1 then begin
      numerator=numerator+list_int(i,*)/list_wvl(i)
   endif
endfor
;
denominator=fltarr(ntemp)
for i=0,nlist-1 do begin
   if dstatus(i) eq 1 then begin
      denominator=denominator+list_int(i,*)/list_wvl(i)
   endif
endfor
;
ratiohi=numerator/denominator
;
;ratiohi=list_int(line1,*)*list_wvl(line2)/(list_int(line2,*)*list_wvl(line1))
;ratiohi=reform(ratiohi)
;
;oplot,temperature, ratiohi,linestyle=2

;
for itemp=0,ntemp-1 do begin    ; loop over densities
;
   temperature(itemp)=10.^(tempmin+itemp*dlogtemp)
;
;
   if keyword_set(pressure) then begin
      edensity=pressure/temperature(itemp)
   endif else if keyword_set(density) then begin
      edensity=density
   endif else edensity=default_density
;
   edensity=0.5*(edensity)
;
;
   pop_solver,temperature[itemp],edensity,pop
   pop = reform(pop)
;
   for ilist=0,nlist-1 do begin
      list_int(ilist,itemp)=pop(list_lvl2(ilist))*list_a(ilist)
   endfor
;
endfor                          ;   loop over temperatures
;
numerator=fltarr(ntemp)
for i=0,nlist-1 do begin
   if nstatus(i) eq 1 then begin
      numerator=numerator+list_int(i,*)/list_wvl(i)
   endif
endfor
;
denominator=fltarr(ntemp)
for i=0,nlist-1 do begin
   if dstatus(i) eq 1 then begin
      denominator=denominator+list_int(i,*)/list_wvl(i)
   endif
endfor
;
ratiolo=numerator/denominator
;
;ratiolo=list_int(line1,*)*list_wvl(line2)/(list_int(line2,*)*list_wvl(line1))
;ratiolo=reform(ratiolo)
;
oplot,temperature, ratiolo,linestyle=2
oplot,temperature,ratiohi,linestyle=2 ;linestyle=1
;

print, ' '
print,  ' The dashed lines overplotted are the ratios calculated at twice and half the density'
print, ' '

if keyword_set(psfile) then begin
   dname = !d.name
   set_plot,'ps'
   device,filename=psfile
   landscape
   !p.font=0
   device,/times,/isolatin1,font_size=16
   thick=3
   plot_oi,temperature,ratio,xr=[temperature(0),temperature(ntemp-1)],title=description, $
     xthick=thick,ythick=thick,thick=thick, $
     xtitle=xtitle,ytitle=ytitle , xstyle=1, ystyle=1 

   if keyword_set(pressure) then begin
      xyouts,temperature(1),max(ratio)/2.,'Constant Pressure = '+string(pressure,'(e10.2)')+' (cm!e-3!n !eo!nK)',$
        charsiz=charsiz
   endif else if keyword_set(density) then BEGIN
      xyouts,temperature(1),max(ratio)/2.,'Constant Density = '+string(density,'(e10.2)')+' (cm!e-3!n)',$
        charsiz=charsiz
   endif else begin
      xyouts,temperature(1),max(ratio)/2.,'Constant Density = '+string(default_density,'(e10.2)')+' (cm!e-3!n)',$
        charsiz=charsiz
   endelse
;
   !p.font=-1
   device,/close
   set_plot,dname	;Return to previous device
endif
;

end