;+
; Project     : SOHO - CDS     
;                   
; Name        : CALC_DMM_DR
;               
; Purpose     : Calculates CHIANTI density sensitive line ratios.
;               
; Explanation : 
;               
; Use         : Called from DMM_NE
;    
; Inputs      :  ciz  - element number
;                cion - ion number
;                denmin, denmax - log density range
;                mess_win - message widget ID for DMM_NE use.
;               
; Opt. Inputs :  None
;               
; Outputs     :  None
;               
; Opt. Outputs:  None
;               
; Keywords    :  Temperature - if set use this temperature rather than 
;                              temp at max of ioneq.
;
; Calls       :  To read_xxx files for CHIANTI database
;
; Common      :
;               
; Restrictions:  None
;               
; Side effects:  None
;               
; Category    :  Spectral
;               
; Prev. Hist. :  Based on density_ratios by K Dere.
;
; Written     :  C D Pike, RAL, 22-Jan-96
;               
; Modified    :  Send warning message to widget.  CDP, 27-Jan-96
;                Update generally.  CDP, 7-Jun-97
;                Corrected typo in reading elvlc file.  CDP, 14-Jul-97
;                Added minimum line ratio factor.  CDP, 17-Jul-97
;                Added multiple wavelength ranges. CDP, 18-Jul-97
;                Fix typo for XXII stage.          CDP, 05-Mar-98
;                Update list of elements.          CDP, 18-Jun-99
;                v. 9. Correct hiccup in above.          CDP, 13-Jul-99
;
;                Update for compatibility with CHIANTI v.3. Replaced populate
;                with pop_solver. Left out 'dielectronic' lines. ouput
;                temperature. plus various small additions.
;                Now ratios are calculated with  0.1 increment in log Ne
;                Giulio Del Zanna  10-Oct-2000 
;
; Version     :  Version 10, 10-Oct-2000 
;-                             

pro calc_dmm_dr,ciz,cion,denmin,denmax,mess_win,$
        temperature=temperature

;
;  denmin= log 10 of minimum density to be considered
;  denmax= log 10 of maximum density to be considered
;

;these COMMONS are used by pop_solver:

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 with the caller routine, chianti_ne.pro :

common dmm_lines, list_wvl, list_int, list_descr1, list_descr2, species,$
  density, nden, ratio, description, nlist, savetext,$
  temper, intens

;commmon with the read_* routines:
common dmm_refs, ioneq_ref, wgfaref, copy_eref, upsref

;common with the caller routine, chianti_ne.pro :

common wavmm, wminw1, wminw2, wminw3, wminw4, $
  wmaxw1, wmaxw2, wmaxw3, wmaxw4, $
  minw1, minw2, minw3, minw4, $
  maxw1, maxw2, maxw3, maxw4

defsysv,'!xuvtop', EXISTS = EXISTS 
IF NOT EXISTS THEN use_chianti
defsysv,'!xuvtop', EXISTS = EXISTS 
IF NOT EXISTS THEN $
message, 'system variable !xuvtop must be set using e.g. use_chianti,directory '
xuvtop = !xuvtop


element=['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']

;
ionstage=['I','II','III','IV','V','VI','VII','VIII','IX','X','XI','XII',$
          'XIII','XIV','XV','XVI','XVII','XVIII','XIX','XX','XXI','XXII',$
          'XXIII','XXIV','XXV','XXVI','XXVII','XXVIII']


iz = where(strlowcase(element) eq strlowcase(ciz))
if iz(0) eq -1 then begin
   bell
   if n_elements(mess_win) gt 0 then begin
      widget_control, mess_win, /append, set_v='Unrecognised element'
   endif else begin 
      print,'Unrecognised element' 
   endelse
   return
endif else iz = iz(0)+1

ion = where(ionstage eq strupcase(cion))
if ion(0) eq -1 then begin
   bell
   if n_elements(mess_win) gt 0 then begin
      widget_control, mess_win, /append, set_v='Unrecognised ionization stage'   
   endif else begin 
      print,'Unrecognised ionization stage'   
   endelse
   return
endif else ion = ion(0)+1

;
;xuvtop = getenv('CDS_SS_DERE')
;defsysv,'!xuvtop',getenv('CDS_SS_DERE')


;
;  Use T at max ioneq or user supplied
;
if temperature gt 0 then begin
   tmax = temperature
endif else begin


   ioneqdir=add_subdir(!xuvtop,'ioneq')
   if !version.os eq 'vms' then begin
      ioneq_name=ioneqdir+!ioneq_file
   endif else ioneq_name = ioneqdir+'/'+!ioneq_file

;ioneq_name=xuvtop+'/ioneq/'+'arnaud_rothenflug_lmf.ioneq'

   read_ioneq,ioneq_name,ioneq_t,ioneq,ioneq_ref

;dlnt=alog(10.^(ioneq_t(1)-ioneq_t(0)))
;

;   this_ioneq=ioneq(*,iz-1,ion-1+dielectronic)

   this_ioneq=ioneq(*,iz-1,ion-1)
   hit=where(this_ioneq eq max(this_ioneq))
   tmax=ioneq_t(hit)
;
   if n_elements(mess_win) gt 0 then begin
      widget_control, mess_win, $
        set_val='  deriving maximum T from the ionization eq. file: '+$
        !ioneq_file
   END 
   print,  ' deriving maximum T from the ionization eq. file: '+$
     !ioneq_file

   tmax=10.^(tmax(0))

;define temperature to pass on to caller
   temperature = tmax

endelse

;tmax=10.^(tmax(0))

if n_elements(mess_win) gt 0 then begin
   widget_control, mess_win, $
     set_val= string('Using Tmax= ',tmax,format='(a12,e10.3)'), /append
ENDIF
print,'Using Tmax=',tmax,format='(a12,e10.3)'

;
;
;dlogden=0.25                    ;  increment in log Ne for calculating ratios

dlogden=0.1                    ;  increment in log Ne for calculating ratios

;
nden=fix((denmax-denmin)/dlogden) +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']
;
;

species=element(iz-1)+' '+ionstage(ion-1)
; if dielectronic=1 then 

;ion2spectroscopic,ions,species
;
;
zion2filename,iz,ion,fname      ;diel=
wname=fname+'.wgfa'
;elvlname=fname+'.elvl'
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
copy_eref = eref

mult=2.*jj+1.
;
g=where(ecm eq 0.)
;
if(max(g) gt 0) then ecm(g) = ecmth(g)
;
;
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,nden)
;
nlist=0
;
;   find out the number of excited levels
;
;populate,tmax,10.^denmin,pop

;CALL to pop_solver:

pop_solver,tmax,10.^denmin,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
;
;

;also CHECK FOR LINES WITH WAVELENGTH=0 
 
   if ( ((ww ge minw1) and (ww le maxw1)) or $
        ((ww ge minw2) and (ww le maxw2)) or $
        ((ww ge minw3) and (ww le maxw3)) or $
        ((ww ge minw4) and (ww le maxw4)) ) and $ 
     (ww NE 0) 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

;
;  If no lines found then exit
;
if nlist eq 0 then begin
   bell
   if n_elements(mess_win) gt 0 then begin
      widget_control, mess_win, set_val='No lines found in that region'
   endif else begin
      print,'No lines found in that region'
   endelse
   return
endif
;
;
;
;   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,nden)

;
;  for sorted list, calculate level populations and intensities
;

density = 10.^(denmin+findgen(nden)*dlogden)
pop_solver,tmax,density,pop
pop = reform(pop)
FOR ilist=0,nlist-1 DO $
   list_int[ilist,*]=pop[*,list_lvl2[ilist]]*list_a[ilist]/density


;density=fltarr(nden)
;for iden=0,nden-1 do begin   ; loop over densities
;   density(iden)=10.^(denmin+iden*dlogden)
;   populate,tmax,density(iden),pop
;   for ilist=0,nlist-1 do begin
;      list_int(ilist,iden)=pop(list_lvl2(ilist))*list_a(ilist)/density(iden)
;   endfor
;endfor;   loop over densities

;widget_control, mess_win, /append, $
;  set_v=string(' Tmax=',tmax,format='(a8,e10.3)')


;now exclude all lines weaker than  intens*maximum(intensity)

dd = list_int(*,nden-1)
ddmax = max(dd)
n = where(dd ge intens*ddmax)
list_wvl = list_wvl(n)
list_int = list_int(n,*)
list_descr1 = list_descr1(n)
list_descr2 = list_descr2(n)
nlist = n_elements(list_wvl)


;  print results

dash=' - '
;
widget_control, mess_win, /append,$
  set_val='----------------------------------------'
print,'------------------------------------------'

widget_control, mess_win, set_val= '#    Wavelength  Intensity    Transition ',/append
widget_control, mess_win, set_val= '        (A)   ',/append
widget_control, mess_win, set_val= ' ',/append

print, '#    Wavelength  Intensity    Transition '
print, '        (A)   '
print, ''

savetext = strarr(nlist)

for i=0,nlist-1 do begin
   text = string(i,list_wvl(i),list_int(i,nden-1),list_descr1(i),dash,$
                 strpad(list_descr2(i),15,/after), $
                 format='(i4,f10.4,e10.2,a15,a3,a15)')
   savetext(i) = text
   widget_control, mess_win, set_val=text,/append
   print,text
endfor
end