PRO fe12_density, fit195, fit186, dens, sum195=sum195, yoffset=yoffset, $
                  use192=use192, indens=indens

;+
; NAME
;
;     FE12_DENSITY
;
; PROJECT
;
;     Hinode/EIS
;
; EXPLANATION
;
;     Takes the Fe XII 195 and 186 fit structures produced by eis_auto_fit,
;     and creates a density structure containing a density map. The densities
;     are derived used the CHIANTI atomic database.
;
; INPUTS
;
;     FIT195  Fit structure for Fe XII 195 line.
;
;     FIT186  Fit structure for Fe XII 186 line.
;
; OPTIONAL INPUTS
;
;     YOFFSET Works out densities assuming a spatial offset between 195 and
;             186. YOFFSET has to be positive which indicates that the 186
;             image lies *above* the 195 image on the detector.
;
;     INDENS  Instead of using the computed density to derive the column
;             depth, the density contained in INDENS is used. E.g., the
;             density derived from Fe XIII 203/202 or Fe XII 186/195 could
;             be used as an input.
;
; KEYWORDS
;
;     SUM195  If set, then the 195 line is assumed to comprise of the 195.12
;             and 195.18 lines when computing the densities.
;
;     USE192  Instead of using the 195 line as the dominator, instead use the
;             192.39 line.
;
; OUTPUTS
;
;     DENS    Structure containing derived densities. Tags are:
;               .ldens   Log10 of density.
;               .dens    Density
;               .lo      Lower limit on density (log10 value)
;               .hi      Upper limit on density (log10 value)
;               .ratio   Ratio values
;               .hiperc  Percentage upper limit
;               .loperc  Percentage lower limit
;
; CALLS
;
;     ERRORS, EMISS_CALC
;
; HISTORY
;
;     Ver.1, 18-Oct-2007, Peter Young
;     Ver.2, 9-Nov-2007, Peter Young
;         updated header and tidied up routine
;     Ver.3, 12-Mar-2008, Peter Young
;         added /sum195 keyword
;     Ver.4, 13-Mar-2008, Peter Young
;         added YOFFSET keyword
;     Ver.5, 28-Mar-2008, Peter Young
;         added /sum195 and /use192 keywords
;     Ver.6, 17-Apr-2008, Peter Young
;         added indens= optional input
;     Ver.7, 2-Aug-2010, Peter Young
;         added a number of reform() commands to make the routine
;         compatible with the fit structures produced by the new
;         version of eis_auto_fit
;-

siz=size(reform(fit186.int))

IF n_elements(yoffset) EQ 0 THEN yoffset=0

int186=reform(fit186.int)
int195=reform(fit195.int)
err186=reform(fit186.interr)
err195=reform(fit195.interr)

IF yoffset NE 0 THEN BEGIN
  ny=siz[2]
  int186=int186[*,yoffset:*]
  int195=int195[*,0:ny-1-yoffset]
  err186=err186[*,yoffset:*]
  err195=err195[*,0:ny-1-yoffset]
ENDIF

siz=size(int186)

;
; Get missing data value
;
missing=fit195.missing

;
; Calculate ratio and 1-sigma errors
;
errors,int186,err186,int195,err195,rat,err

;
; Make sure missing data are flagged.
;
i_miss=where(int195 EQ missing OR int186 EQ missing,ni)
IF ni NE 0 THEN rat[i_miss]=missing

;
; The theoretical Fe XII ratio is calculated over a fine density grid.
;
dens=findgen(501)/100.+7.5

;
; Get line emissivities from CHIANTI and extract the values for the two
; lines under consideration.
;
; Note that the 186 line is a blend of two Fe XII lines.
;
; WARNING: if the line indices change in future versions of CHIANTI, then
; this routine will need to be modified!
;
em=emiss_calc(26,12,temp=6.1,dens=dens,/quiet)

IF keyword_set(sum195) THEN BEGIN
  em195=reform(em[[873,875]].em)
  em195=total(em195,2)
ENDIF ELSE BEGIN
  em195=reform(em[873].em)
ENDELSE
;
IF keyword_set(use192) THEN BEGIN
  em195=reform(em[827].em)
ENDIF
;
em186=em[749:750].em
em186=reform(em186)
em186=total(em186,2)
;
rat_theory=em186/em195

outdens=fltarr(siz[1],siz[2])
hidens=fltarr(siz[1],siz[2])
lodens=fltarr(siz[1],siz[2])
emiss=fltarr(siz[1],siz[2])

FOR i=0,siz[1]-1 DO BEGIN
  FOR k=0,siz[2]-1 DO BEGIN
    IF rat[i,k] GT 0. THEN BEGIN
      getmin=min(abs(rat_theory-rat[i,k]),j)
      outdens[i,k]=dens[j]
     ;
     ; re-compute emissivity if INDENS has been specified
     ;
      IF n_elements(indens) NE 0 THEN BEGIN
        getmin=min(abs(indens.ldens[i,k]-dens),j)
        emiss[i,k]=em186[j]+em195[j]
      ENDIF ELSE BEGIN
        emiss[i,k]=em186[j]+em195[j]
      ENDELSE
;      emiss[i,k]=em186[j]+em195[j]
   ;
      getmin=min(abs(rat_theory-(rat[i,k]+err[i,k])),j)
      hidens[i,k]=dens[j]
   ;
      getmin=min(abs(rat_theory-(rat[i,k]-err[i,k])),j)
      lodens[i,k]=dens[j]
    ENDIF ELSE BEGIN
      outdens[i,k]=-1.
    ENDELSE
  ENDFOR
ENDFOR

f_OF_t=0.294
ab_fe=1.26e-4

;
; for denserr I'm using only hidens. I should really the max of hidens and
; lodens...
;
interr=sqrt(reform(fit186.interr)^2+reform(fit195.interr)^2)

;
; Allow an alternative density profile to be input when computing column depth.
;
IF n_elements(indens) NE 0 THEN BEGIN
  dd=indens.ldens
  denserr=10.^(indens.hi)-indens.dens
ENDIF ELSE BEGIN
  dd=outdens
  denserr=10.^(hidens)-10.^(outdens)
ENDELSE


alpha=4.0*!pi/0.83/ab_fe/f_OF_t / emiss

errors,alpha*(reform(fit186.int)+reform(fit195.int)),alpha*interr,10.^dd,denserr,coldepth,coldepth_err


dens={ldens: outdens, $
      dens: 10.^outdens, $
      lo: lodens, $
      hi: hidens, $
      ratio: rat, $
      hiperc: hidens-hidens, $
      loperc: lodens-lodens, $
      coldepth: coldepth, $
      coldepth_err: coldepth_err}

dens.hiperc=abs(dens.dens-10.^(dens.hi))/dens.dens *100.
dens.loperc=abs(dens.dens-10.^(dens.lo))/dens.dens *100.

END