function eis_chianti_dens_ratio, input, denom, numer, drange=drange, temp=temp

;+
; NAME:
;      EIS_CHIANTI_DENS_RATIO
;
; PROJECT:
;      Hinode/EIS
;
; PURPOSE:
;      Computes the theoretical variation of a density sensitive line
;      ratio and outputs it to an IDL structure. The output is
;      intended to be sent to the routine eis_density.pro.
;
;      The CHIANTI database is used to generate the theoretical ratios
;      so the user must have the CHIANTI package in their Solarsoft
;      distribution.
;
;      eis_chianti_dens_ratio can be called in two ways. If the ratio
;      to be studied is a standard EIS diagnostic then the routine may
;      already contain the parameters necessary to calculate the
;      ratio. Simply call the routine without any inputs and a list of
;      the standard diagnostics will appear. Select the ratio you want
;      and that is it.
;
;      For non-standard diagnostics you will need to know the CHIANTI
;      indices of the wavelengths of the diagnostic. This is more
;      complicated and an example is shown below.
;
; INPUTS:
;      INPUT   Can either be the index of a standard EIS diagnostic,
;              or the CHIANTI format name of an ion, e.g., 'fe_13' for
;              Fe XIII. If INPUT is not specified then a menu will be
;              printed to the screen asking the user to select one of
;              the standard EIS diagnostics.
;
;      DENOM   If an ion name has been specified for INPUT, then DENOM
;              needs to be set to the CHIANTI index/indices of the
;              denominator line of the ratio. See the example below
;              for how to determine the index.
;
;      NUMER   If an ion name has been specified for INPUT, then NUMER
;              needs to be set to the CHIANTI index/indices of the
;              numerator line of the ratio. See the example below
;              for how to determine the index.
;
; OPTIONAL INPUTS:
;      DRANGE  By default the routine computes ratio values for
;              densities between 10^8 and 10^12. By specifying DRANGE
;              this range can be modified. E.g., DRANGE=[7,10] sets
;              the range to 10^7 to 10^10.
;
;      TEMP    The routine automatically computes the ratio for the
;              temperature of maximum ionization of the ion (as
;              determined from the default CHIANTI ion balance
;              file). TEMP can be used to specify an alternative
;              value. Note the LOG value must be specified.
;
; OUTPUTS:
;      A structure with the following tags (using Si X 258/261 as an
;      example): 
;
;      DENS            FLOAT     Array[41]
;      RATIO           DOUBLE    Array[41]
;      NUM_WVL         STRING    '258.374'
;      DEN_WVL         STRING    '261.057'
;      ION             STRING    'si_10'
;      EM_SUM          DOUBLE    Array[41]
;      TEMP            DOUBLE           6.1500001
;      NUM_IND         INT           2378
;      DEN_IND         INT           2403
;      CHIANTI_VER     STRING    '7.1'
;      TIME_STAMP      STRING    'Mon Dec  3 11:44:22 2012'
;      QUICK           INT              0
;
; EXAMPLES:
;      First do:
;
;      IDL> ratio=eis_chianti_dens_ratio()
;
;      which will print a list of the standard diagnostics. If your
;      ratio is listed, then just type the index. To avoid the
;      interactive input one can do, e.g.,
;
;      IDL> ratio=eis_chianti_dens_ratio(3)
;
;      which will load the Mg VII 280.74/278.40 diagnostic.
;
;      If your ratio is not one of the standard CHIANTI ones, then the
;      following procedure needs to be applied. Consider the Fe XIII
;      ratio 203.80+203.83/202.04.
;
;      Do:
;
;      IDL> em=emiss_calc(26,13)
;      IDL> em1=emiss_select(em,wrange=[201,203],sel=sel)
;
;      a widget appears showing all Fe XIII lines between 201 and 203
;      angstroms. Choose the transition corresponding to the 202.04
;      line. The output 'sel' contains the CHIANTI reference index for
;      the transition. For CHIANTI 7.1 the index is 39127, but this
;      may change in future versions.
;
;      Do the same thing for the 203.80+203.83 self-blend by
;      specifying wrange=[203,205]. This time you will have to choose
;      the two transitions. The output 'sel' should contain the
;      numbers 39479 and 39485 (for CHIANTI 7.1).
;
;      With these index values we now call eis_chianti_dens_ratio:
;
;      IDL> ratio=eis_chianti_dens_ratio('fe_13',39127,[39479,39485])
;
;      where it can be seen that the indices for the denominator and
;      numerator lines are specified in the 2nd and 3rd inputs,
;      respectively. The blend is identified by specifying a 2 element
;      array. 
;
; CALLS:
;      EMISS_CALC, GET_TMAX, CONVERTNAME
;
; MODIFICATION HISTORY:
;      Version 1, 4-May-2011, Peter Young
;        added Fe X 184/257 ratio (no. 8)
;      Version 2, 3-Feb-2012, Peter Young
;        updated code to reflect new atomic models in CHIANTI 7, but
;        also retained compatibility with CHIANTI 6.
;      Version 3, 3-Dec-2012, Peter Young
;        updated for CHIANTI 7.1; added chianti_ver and time_stamp
;        tags to output structure; updated header.
;      Version 4, 18-Mar-2013, Peter Young
;        added Fe XI 182.17/188.22.
;      Version 5, 25-Mar-2013, Peter Young
;        the output structure now contains the 'quick' tag; the output
;        ratios are now calculated at 0.2 dex intervals (instead of
;        0.1) in order to speed up the routine for the coronal iron
;        ions (there should be minimal loss of accuracy in doing
;        this). 
;      Version 6, 11-Jun-2013, Peter Young
;        added Ar XIV ratio.
;      Version 7, 9-Apr-2014, Peter Young
;        added Ca XV 181/182 ratio.
;      Version 8, 15-Jul-2014, Peter Young
;        updated header; no change to code.
;      Version 9, 30-Jun-2015, Peter Young
;        updated ratios for CHIANTI 8.
;      Version 10, 24-May-2016, Peter Young
;        it turns out that the emiss_calc index does not uniquely
;        specify a transition when two or more transitions have the
;        same wavelength; I've implemented a fix for CHIANTI 8;
;        thanks to Peter Gomory for identifying this problem.
;      Version 11, 6-Jun-2019, Peter Young
;        I've updated the transition information for CHIANTI 9.
;      Version 12, 03-Mar-2021, Peter Young
;        Updated for CHIANTI 10.
;      Version 13, 27-Jun-2023, Peter Young
;        Updated for CHIANTI 10.1; fixed bug for ratio 9 (Fe X 175/174).
;-

;
; First check the version of CHIANTI that user has. 
;
vfile=concat_dir(!xuvtop,'VERSION')
chck=file_search(vfile)
IF chck[0] NE '' THEN BEGIN
  openr,lin,vfile,/get_lun
  a=0.0
  readf,lin,a
  free_lun,lin
  CASE 1 OF
    a GE 6 AND a LT 7: vnum='6'
    a GE 7 AND a LT 7.1: vnum='7'
    a GE 7.1 AND a LT 8: vnum='7.1'
    a GE 8.0 AND a LT 9: vnum='8'
    a GE 9.0 AND a LT 10: vnum='9'
    a GE 10.0 AND a LT 10.1: vnum='10'
    a GE 10.1: vnum='10.1'
    ELSE: BEGIN
      print,'%EIS_CHIANTI_DENS_RATIO: This routine is not compatible with CHIANTI version '+trim(a)+'.'
      print,'                         Please contact Peter Young to update it for this version.'
      return,-1
    END 
  ENDCASE 
ENDIF ELSE BEGIN
  print,'%EIS_CHIANTI_DENS_RATIO: it seems that CHIANTI is not part of your Solarsoft distribution.'
  print,'                         Please check your setup. Returning...'
  return,-1
ENDELSE 


IF n_elements(drange) EQ 0 THEN drange=[8.0,12.0]

ans=0

IF n_elements(input) EQ 0 THEN BEGIN
  print,''
  print,'Please choose a density sensitive ratio from the list of standard diagnostics below:'
  print,'   1.  Fe XIII  (203.80+203.83)/202.04'
  print,'   2.  Fe XII   (186.85+186.89)/195.12'
  print,'   3.  Mg VII   280.74/278.40'
  print,'   4.  Si X     258.37/261.04'
  print,'   5.  Fe XIV   264.79/274.20'
  print,'   6.  Fe XIII  196.54/202.04'
  print,'   7.  Fe XII   196.64/195.12'
  print,'   8.  Fe X     184.54/(257.259+257.263)'
  print,'   9.  Fe X     175.27/174.53'
  print,'  10.  Fe XI    182.17/188.22'
  print,'  11.  Ar XIV   187.97/194.40'
  print,'  12.  Ca XV    181.90/182.87'
  print,''
  read,'Input number (type 0 to exit): ',ans
  IF datatype(ans) NE 'INT' THEN BEGIN
    print,'**Invalid input. Returning...'
    return,-1
  ENDIF
ENDIF

IF datatype(input) EQ 'INT' THEN BEGIN
  ans=input
ENDIF 

quick=0


;
; My old way of doing the transition identification was to specify
; NUMER and DENOM to be indices within the emiss_calc structure that
; identify the transitions I need. However, when two transitions have
; the same wavelength it's possible that a transition's index will
; vary depending on the computer on which the code runs
; (!). I've thus implemented NUM_I and NUM_J for the problem
; cases of Fe XII and Fe XIII, which specify the level indices (i -
; lower level; j - upper level). These are then used to find the NUMER
; values later in the code, after emiss_calc has been run. 
;
IF datatype(input) EQ 'STR' THEN BEGIN
  IF n_params() LT 3 THEN BEGIN
    print,'** Please check the input parameters **'
    print,'When specifying an ion name, the calling sequence is, e.g.,'
    print,"  IDL> eis_chianti_dens_ratio, 'fe_13', denom, numer"
    print,'where denom and numer are CHIANTI indices identifying the denominator and numerator'
    print,'of the line ratio.'
    return,-1
  ENDIF
  ionname=input
ENDIF ELSE BEGIN
  IF ans GE 1 AND ans LE 12 THEN BEGIN
    CASE ans OF
      1: BEGIN
        ionname='fe_13'
        CASE 1 OF
          vnum EQ '10.1': BEGIN
            den_i=1
            den_j=20
            num_i=[3,3]
            num_j=[24,25]
;            quick=1
          END 
          vnum EQ '8' OR vnum EQ '9' OR vnum EQ '10': BEGIN
            denom=26520
;            numer=[26842,26851]
            num_i=[3,3]
            num_j=[24,25]
            quick=1
          END 
          vnum EQ '7.1': BEGIN
            denom=39127          ; 202.04 line
            numer=[39479,39485]   ; 203.80+203.83 line
            quick=1
          END 
          vnum EQ '7': BEGIN
            denom=1054          ; 202.04 line
            numer=[1080,1081]   ; 203.80+203.83 line
          END 
          ELSE: BEGIN
            denom=9             ; 202.04 line
            numer=[11,12]       ; 203.80+203.83 line
          END 
        ENDCASE 
      END 
     ;
      2: BEGIN
        ionname='fe_12'
        CASE 1 OF
          vnum EQ '9' OR vnum EQ '10' OR vnum EQ '10.1': BEGIN
            denom=30051
            num_i=[2,3]
            num_j=[36,39]
            quick=1
          END 
          vnum EQ '8': BEGIN
            denom=30433
            num_i=[2,3]
            num_j=[36,39]
            quick=1
          END 
          vnum EQ '7.1': BEGIN
            denom=23391           ; 195.12 line
            numer=[21736,21747]   ; 186.85+186.89 line
            quick=1
          END 
          ELSE: BEGIN 
            denom=873           ; 195.12 line
            numer=[749,750]     ; 186.85+186.89 line
          END 
        ENDCASE 
      END
     ;
      3: BEGIN
        ionname='mg_7'
        CASE 1 OF
          vnum EQ '10' OR vnum EQ '10.1': BEGIN 
            denom=5182           ; 278.40 line
            numer=5266           ; 280.75 line
          END 
          vnum EQ '7.1' OR vnum EQ '7' OR vnum EQ '8' OR vnum EQ '9': BEGIN
            denom=231           ; 278.40 line
            numer=232           ; 280.75 line
          END 
          ELSE: BEGIN
            denom=195           ; 278.40 line
            numer=196           ; 280.75 line
          END 
        ENDCASE 
      END
     ;
      4: BEGIN
        ionname='si_10'
        CASE 1 OF
          vnum EQ '8' OR vnum EQ '9' OR vnum EQ '10' OR vnum EQ '10.1': BEGIN
            denom=2976           ; 261.04 line
            numer=2948           ; 258.37 line
          END 
          vnum EQ '7' OR vnum EQ '7.1': BEGIN
            denom=2403           ; 261.04 line
            numer=2378           ; 258.37 line
          END 
          ELSE: BEGIN
            denom=620           ; 261.04 line
            numer=619           ; 258.37 line
          END 
        ENDCASE 
      END 
     ;
      5: BEGIN
        ionname='fe_14'
        CASE 1 OF
          vnum EQ '7.1' OR vnum EQ '8' OR vnum EQ '9' OR vnum EQ '10' OR vnum EQ '10.1': BEGIN
            denom=33276         ; 274.20 line
            numer=32875         ; 264.79 line
            quick=1
         END 
          vnum EQ '7': BEGIN
            denom=11741         ; 274.20 line
            numer=11494         ; 264.79 line
         END 
          ELSE: BEGIN
            denom=115           ; 274.20 line
            numer=106           ; 264.79 line
          END 
        ENDCASE 
      END 
     ;
      6: BEGIN
        ionname='fe_13'
        CASE 1 OF
          vnum EQ '10.1': BEGIN
            den_i=1
            den_j=20
            num_i=4
            num_j=26
;            quick=1
          END 
          vnum EQ '8' OR vnum EQ '9' OR vnum EQ '10': BEGIN
            denom=26520          ; 202.04 line
;            numer=25571          ; 196.54 line
            num_i=4
            num_j=26
            quick=1
          END 
          vnum EQ '7.1': BEGIN
            denom=39127          ; 202.04 line
            numer=38021          ; 196.54 line
            quick=1
          END 
          vnum EQ '7': BEGIN
            denom=1054          ; 202.04 line
            numer=974           ; 196.54 line
          END 
          ELSE: BEGIN
            denom=9             ; 202.04 line
            numer=5             ; 196.54 line
          END 
        ENDCASE 
      END
     ;
      7: BEGIN
        ionname='fe_12'
        CASE 1 OF
           vnum EQ '10' OR vnum EQ '10.1': BEGIN
              den_i=1 & den_j=27
              num_i=3 & num_j=34
              quick=1
           END
          vnum EQ '9': BEGIN
            denom=30443           ; 195.12 line
;            numer=30969           ; 196.64 line
            num_i=3
            num_j=34
            quick=1
          END 
          vnum EQ '8': BEGIN
            denom=30051           ; 195.12 line
;            numer=30969           ; 196.64 line
            num_i=3
            num_j=34
            quick=1
          END 
          vnum EQ '7.1': BEGIN 
            denom=23391           ; 195.12 line
            numer=23707           ; 196.64 line
            quick=1
          END 
          ELSE: BEGIN
            denom=873           ; 195.12 line
            numer=897           ; 196.64 line
          ENDELSE 
        ENDCASE 
      END
     ;
      8: BEGIN
        ionname='fe_10'
        CASE 1 OF
           vnum EQ '10' OR vnum EQ '10.1': BEGIN
              den_i=[1,1] & den_j=[4,5]
              num_i=1 & num_j=27
              quick=1
           END 
          vnum EQ '8' OR vnum EQ '9': BEGIN
            denom=[13179,13180]   ; 257.259+257.263
            numer=6035           ; 184.53
            quick=1
          END 
          vnum EQ '7.1': BEGIN
            denom=[25846,25849]   ; 257.259+257.263
            numer=13453           ; 184.53
            quick=1
          END 
          ELSE: BEGIN 
            denom=[2081,2082]   ; 257.259+257.263
            numer=1505          ; 184.53
          END 
        ENDCASE 
      END 
     ;
      9: BEGIN
        ionname='fe_10'
        CASE 1 OF
           vnum EQ '10' OR vnum EQ '10.1': BEGIN
              den_i=1 & den_j=30
              num_i=2 & num_j=31
              quick=1
           END 
          vnum EQ '8' OR vnum EQ '9': BEGIN
            denom=5100          ; 174.53 line
            numer=5154          ; 175.27 line
            quick=1
          END 
          vnum EQ '7.1': BEGIN
            denom=11446          ; 174.53 line
            numer=11568          ; 175.27 line
            quick=1
          END 
          ELSE: BEGIN 
            denom=1381          ; 174.53 line
            numer=1392          ; 175.27 line
          END 
        ENDCASE 
      END 
     ;
      10: BEGIN
        ionname='fe_11'
        CASE 1 OF
           vnum EQ '10' OR vnum EQ '10.1': BEGIN
              num_i=2
              num_j=43
              den_i=1
              den_j=38
              quick=1
           END 
          vnum EQ '8' OR vnum EQ '9': BEGIN
            numer=20740         ; 182.17 line
            denom=22438         ; 188.22 line
            quick=1
          END 
          '7.1': BEGIN
            numer=18769         ; 182.17 line
            denom=20298         ; 188.22 line
            quick=1
          END 
        ENDCASE 
      END 
     ;
      11: BEGIN
        ionname='ar_14'
        CASE 1 OF
           vnum EQ '10' OR vnum EQ '10.1': BEGIN
              num_i=2 & num_j=10
              den_i=1 & den_j=8
           END 
          vnum EQ '8' OR vnum EQ '9': BEGIN
            numer=3868           ; 187.97 line
            denom=3895           ; 194.40 line
          END 
          ELSE: BEGIN 
            numer=619           ; 187.97 line
            denom=621           ; 194.40 line
          END 
        ENDCASE
      END 
     ;
      12: BEGIN
         ionname='ca_15'
        CASE 1 OF
           vnum EQ '10' OR vnum EQ '10.1': BEGIN
              num_i=3 & num_j=12
              den_i=3 & den_j=11
           END 
           ELSE: BEGIN 
              numer=290           ; 181.90 line
              denom=292           ; 182.87 line
           END
        ENDCASE 
      END 
      ;;**IF YOU ADD EXTRA RATIOS TO THE CASE STATEMENT, PLEASE ALSO
      ;;UPDATE THE IF STATEMENT (ANS LE 11).
    ENDCASE 
  ENDIF ELSE BEGIN
    print,''
    print,'%EIS_CHIANTI_DENS_RATIO: Exit selected. Returning...'
    return,-1
  ENDELSE 
ENDELSE 

dd=0.2    ; this is the density interval in dex (log density)
convertname,ionname,iz,ion
nd=(drange[1]-drange[0])/dd+1
dens=findgen(nd)*dd+drange[0]
;
IF n_elements(temp) EQ 0 THEN temp=alog10(get_tmax(ionname))
;
em=emiss_calc(iz,ion,dens=dens,temp=temp,/quiet)

;
; The following deals with the case when DEN_I and DEN_J have been
; specified instead of DENOM.
;
IF n_elements(denom) EQ 0 THEN BEGIN
  n=n_elements(den_i)
  denom=lonarr(n)
  FOR i=0,n-1 DO BEGIN
    k=where(em.level1 EQ den_i[i] AND em.level2 EQ den_j[i],nk)
    denom[i]=k[0]
  ENDFOR 
ENDIF 
em_denom=em[denom].em
em_denom=reform(em_denom)
IF n_elements(denom) GT 1 THEN em_denom=total(em_denom,2)

;
; The following deals with the case when NUM_I and NUM_J have been
; specified instead of NUMER.
;
IF n_elements(numer) EQ 0 THEN BEGIN
  n=n_elements(num_i)
  numer=lonarr(n)
  FOR i=0,n-1 DO BEGIN
    k=where(em.level1 EQ num_i[i] AND em.level2 EQ num_j[i],nk)
    numer[i]=k[0]
  ENDFOR 
ENDIF 
em_numer=em[numer].em
em_numer=reform(em_numer)
IF n_elements(numer) GT 1 THEN em_numer=total(em_numer,2)

ratio=em_numer/em_denom

;
; Create a string containing the numerator wavelengths (separated by +)
;
wvl=em[numer].lambda
num_wvl=''
FOR i=0,n_elements(wvl)-1 DO BEGIN
  wvls=trim(string(format='(f15.3)',wvl[i]))
  IF num_wvl EQ '' THEN num_wvl=wvls ELSE num_wvl=num_wvl+'+'+wvls
ENDFOR
 
;
; Create a string containing the denominator wavelengths (separated by +)
;
wvl=em[denom].lambda
den_wvl=''
FOR i=0,n_elements(wvl)-1 DO BEGIN
  wvls=trim(string(format='(f15.3)',wvl[i]))
  IF den_wvl EQ '' THEN den_wvl=wvls ELSE den_wvl=den_wvl+'+'+wvls
ENDFOR 

;
; The sum of the line emissivities is needed to compute the column depth
;
em_sum=em_numer+em_denom

;
; This is output structure that is intended to be sent to eis_density.pro.
;
ratio={dens: dens, ratio: em_numer/em_denom, num_wvl: num_wvl, den_wvl: den_wvl, $
       ion: ionname, em_sum: em_sum, temp: temp, num_ind: numer, den_ind: denom, $
       chianti_ver: vnum, time_stamp: systime(), quick: quick }

return,ratio

END