;+
; Project     : SOHO - CDS     
;                   
; Name        : CALC_DMM_TR
;               
; Purpose     : Calculates CHIANTI temperature sensitive line ratios.
;               
; Explanation : 
;               
; Use         : Called from DMM_TE
;    
; Inputs      :  ciz  - element number
;                cion - ion number
;                wmin, wmax - wavelength range
;                temmin, temmax - log temperature range
;                mess_win - message widget ID for DMM_TE use.
;               
; Opt. Inputs :  None
;               
; Outputs     :  None
;               
; Opt. Outputs:  None
;               
; Keywords    :  None
;
; Calls       :  To read_xxx files for CHIANTI database
;
; Common      :
;               
; Restrictions:  None
;               
; Side effects:  None
;               
; Category    :  Spectral
;               
; Prev. Hist. :  Based on temperature_ratios by K Dere.
;
; Written     :  C D Pike, RAL, 22-Jan-96,
; Modified    :  H E Mason, 03-Oct-96 (density to temperature)
;               
; Modified    :  Send warning message to widget.  CDP, 27-Jan-96
;                Update list of elements and cut call to 
;                       read_elvl.                CDP, 18-Jun-99
;                V.3. Fix hiccupp in above.            CDP, 13-Jul-99
;
;                V.4. Rewritten completely, adding possibility to pass on the
;                density at which the intensities are calculated, and making
;                this routine compatible ith CHIANTI v.3.
;                Giulio Del Zanna (DAMTP), 10 Oct-2000
;
;
; Version     :  Version 4, 10 Oct-2000
;-                             

pro calc_dmm_tr,ciz,cion,wmin,wmax,temmin,temmax,mess_win, $
    density=density

;

;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_te.pro :

common dmm_lines_te, list_wvl, list_int, list_descr1, list_descr2, species,$
                  temperature, ntem, ratio, description, nlist, savetext, $
                  dens, intens

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

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



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


IF  density EQ  0 THEN BEGIN 
density=10.^8
ENDIF 

if n_elements(mess_win) gt 0 then begin
   widget_control, mess_win, $
      set_val='calculating intensities at constant Ne = '+$
      string(density,format='(e10.3)')
END 
;
print,'calculating intensities at constant Ne = '+string(density,format='(e10.3)')

;
dlogtem=0.1       ;  increment in log T for calculating ratios
;
ntem=fix((temmax-temmin)/dlogtem) +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,ntem)
;
nlist=0

;
;   find out the number of excited levels
;
;    populate,10.^temmin,density,pop

;CALL to pop_solver:

pop_solver,10.^temmin,density,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,ntem)
;
;  for sorted list, calculated level populations and intensities
;


temperature=10.^(temmin+findgen(ntem) *dlogtem)

pop_solver,temperature,density,pop
pop = reform(pop)

;help, temperature,density,pop


FOR ilist=0,nlist-1 DO $
   list_int[ilist,*]=pop[*,list_lvl2[ilist]]*list_a[ilist]  ;/temperature


;temperature=fltarr(ntem)
;for item=0,ntem-1 do begin   ; loop over temperatures
;   temperature(item)=10.^(temmin+item*dlogtem)
;   populate,temperature(item),density,pop
;   for ilist=0,nlist-1 do begin
;      list_int(ilist,item)=pop(list_lvl2(ilist))*list_a(ilist)/temperature(item)
;   endfor
;endfor;   loop over temperatures


;widget_control, mess_win, /append,$
;                set_val=string(' Ne = ',density,format='(a8,e10.3)')


;now exclude all lines weaker than  intens*maximum(intensity)
;
dd = list_int(*,ntem-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,ntem-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
print, ''

end