;+
; NAME
;
;    EIS_GET_FITDATA()
;
; PROJECT
;
;    Hinode/EIS
;
; EXPLANATION
;
;    Extracts a line parameter array from a FITDATA Gaussian fit
;    structure (produced by EIS_AUTO_FIT). For example intensity,
;    velocity and line width arrays can be specified by using the
;    keyword inputs /INT, /VEL or /WID.
;
;    Pixels where a Gaussian fit was not possible are flagged with the
;    value FITDATA.MISSING in the output array.
;
;    The returned line width is automatically corrected for the EIS
;    instrumental width. The THERMAL_WID keyword allows a thermal
;    width to also be subtracted. Note that a number of pixels will
;    likely have widths smaller than the combined instrumental and
;    thermal widths due to noise, for example. The widths are set to
;    zero for these pixels.
;
; INPUTS
;
;    FITDATA   A structure returned by the Gaussian fitting routine
;              EIS_AUTO_FIT. 
;
; OPTIONAL INPUTS
;
;    LINE  Choose emission line within fit structure. Only applies if
;          FITDATA contains more than one emission line.
;
;    THERMAL_WID  This is a two element array of the form
;                 [LOG10_TEMP,MASS_NUM] where LOG10_TEMP is the
;                 logarithm of the temperature where the ion species
;                 is formed, and MASS_NUM is the mass number of the
;                 ion (e.g., 56 for iron). These parameters allow the
;                 routine to compute the thermal width of the emission
;                 line which is then subtracted from the measured
;                 width. If THERMAL_WID is not specified, then the
;                 thermal width is not subtracted. For Fe XII, the
;                 input would be THERMAL_WID=[6.20,56]. Note that this
;                 keyword is only active if the /WID keyword is set.
;
;     YSHIFT_WVL  This keyword allows the user to shift the output
;                 image up or down (a Y-shift) within the array by an
;                 amount that depends on the wavelength
;                 YSHIFT_WVL. This essentially allows the user to
;                 align the image relative to an array obtained from
;                 another emission line. The example given below shows
;                 how the keyword is used.
;
;     CALIB  This allows the EIS calibration to be modified, changing
;            the output intensity and error arrays. One of the
;            following values should be set:

;            0 - no effect
;            1 - Del Zanna (2013) calibration is implemented.
;            2 - Del Zanna (2013) calibration, and the
;                /correct_sensitivity correction is "undone".
;            3 - Warren et al. (2014) calibration is implemented.
;            4 - Warren et al. (2014) calibration, and the
;                /correct_sensitivity correction is "undone". 
;            5 - Del Zanna et al. (2025) calibration is implemented.
;            6 - Del Zanna et al. (2025)  calibration, and the
;                /correct_sensitivity correction is "undone". 
;
;            WARNING: before using this keyword, the user should know
;            what calibration options were applied in the original
;            call to eis_prep. See the EIS wiki for more details
;            (http://solarb.mssl.ucl.ac.uk:8080/eiswiki/Wiki.jsp?page=EISCalibration)
;
; KEYWORDS
;
;    INT   Return line intensity. Units depend on original call to
;          eis_prep. Typically they will be erg/cm^2/s/sr.
;
;    VEL   Return line velocity. Units: km/s.
;
;    WID   Return FWHM of line, with the instrumental width
;          subtracted. Units: angstroms.
;
;    VWID  Returns width in km/s units. In the case that the non-thermal width
;          is returned, this is the non-thermal velocity (usually assigned Greek
;          letter xi).
;
;    CEN   Return line centroid. Units: angstroms.
;
;    BG1   Return 1st background parameter. Units depend on original call to
;          eis_prep. Typically they will be erg/cm^2/s/sr/angstrom.
;
;    BG2   Return 2nd background parameter. Units depend on original call to
;          eis_prep. Typically they will be erg/cm^2/s/sr/angstrom.
;
;    CHI   Return reduced chi^2 array.
;
;    KEEP_INST_WIDTH  By default the routine removes the instrumental
;                     line width from the measured width. Setting this
;                     keyword stops this correction. If the routine is being used
;                     for data from other instruments (not EIS) then the
;                     instrumental width will not be subtracted so there is no
;                     need to set keep_inst_width.
;
;    QUIET If set, then the routine will not print information
;          messages to the IDL window.
;
;    MAP   If set, then the 2D array is returned within an IDL map
;          structure. 
;
; OPTIONAL OUTPUTS
;
;    ERROR   Contains the error array that matches the user-selected
;            fit parameter.
;
;    UNITS   A string containing the units of the output quantity.
;
;    AV_INST_WIDTH  This contains the average instrumental width for
;                   the spatial pixels that have good fits. It will
;                   only contain a non-zero value if both /WID and
;                   /KEEP_INST_WIDTH have been set. Care should be
;                   taken in using the output value due to the
;                   variation of instrumental width with slit height,
;                   but the value can be useful in some
;                   circumstances. 
;
; OUTPUT
;
;    Returns a 2D array containing the user-selected line
;    parameter. Pixels where a Gaussian fit was not possible are
;    flagged with the value FITDATA.MISSING in the output array.
;
;    If the /MAP keyword was given, then the array is returned within
;    an IDL map structure.
;
; EXAMPLES
;
;    To return a velocity array and the 1-sigma errors on the
;    velocity, do:
;
;    IDL> vel=eis_get_fitdata(fitdata,/vel,error=velerr)
;
;    If FITDATA contains fits for three Gaussians, then the intensity
;    array for the 2nd Gaussian is returned by doing:
;
;    IDL> int=eis_get_fitdata(fitdata,/int,line=1)
;
;    (the first line is Gaussian 0, the second Gaussian 1, etc.).
;
;    ---
;    The YSHIFT_WVL keyword works as follows. Consider fit structures
;    fit200 and fit280 for two emission lines at 200 and 280
;    angstroms. The user wants to form an image ratio from the two
;    lines and this is done as follows:
;
;    IDL> int1=eis_get_fitdata(fit200,/int)
;    IDL> int2=eis_get_fitdata(fit280,/int,yshift_wvl=fit200.refwvl)
;    IDL> plot_image,int2/int1
;
;    There is a spatial offset in the Y-direction between the two
;    wavelengths which can found from the routine EIS_CCD_OFFSET. By
;    specifying yshift_wvl to be the wavelength of the denominator
;    line, the routine applies the pixel shift between 200 and 280
;    angstroms to shift the fit280 image array in order to spatially
;    match the image from the fit200 data. Note that the routine
;    EIS_DENSITY makes use of the YSHIFT_WVL keyword to match the
;    intensity images from two density sensitive lines.
;
; CALLS
;
;    LAMB2V(), EIS_FIT_Y_SHIFT, EIS_LTDS, EIS_RECALIBRATE_INTENSITY
;    EIS_RECALIBRATE_INTENSITY_NEW
;
; INTERNAL ROUTINES
;
;    EGF_SHIFT_ARRAY(), EGF_MAKE_MAP()
;
; HISTORY
;
;    Ver.1, 21-Dec-2009, Peter Young
;    Ver.2, 21-Jan-2010, Peter Young
;      points where the fit parameter has a zero error are now set to
;      be missing.
;    Ver.3, 14-Apr-2010, Peter Young
;      velocity error is now correctly calculated
;    Ver.4, 28-Apr-2010, Peter Young
;      added /CHI keyword
;    Ver.5, 3-May-2010, Peter Young
;      if FITDATA contains 2 or more lines the routine now complains
;      if LINE has not been specified (instead of assuming LINE=0).
;    Ver.6, 24-Mar-2011, Peter Young
;      routine now subtracts the instrumental line width (from
;      EIS_SLIT_WIDTH) from the measured line width.
;    Ver.7, 4-Apr-2011, Peter Young
;      added THERMAL_WID= and /QUIET keywords; modified computation of
;      corrected width to handle thermal width subtraction.
;    Ver.8, 5-Apr-2011, Peter Young
;      update information only printed if keyword /WID is set.
;    Ver.9, 6-Apr-2011, Peter Young
;      warning message is no longer printed when /quiet is set.
;    Ver.10, 7-Apr-2011, Peter Young
;      line= keyword now works again after previous changes
;    Ver.11, 5-May-2011, Peter Young
;      added yshift_wvl= keyword and introduced egf_shift_array
;      internal routine
;    Ver.12, 16-May-2011, Peter Young
;      added /MAP keyword; fixed some problems with setting missing
;      data in the arrays.
;    Ver.13, 9-Mar-2012, Peter Young
;      I mistakenly referred to the atom number in the ion in the
;      description of the thermal width (header) instead of the mass
;      number. This has been corrected. There is no change to the
;      code. 
;    Ver.14, 5-May-2012, Peter Young
;      Added UNITS= optional output; now prints out a help message if
;      no parameters are specified; added AV_INST_WIDTH= optional
;      output. 
;    Ver.15, 4-Apr-2013, Peter Young
;      I'm attempting to make the auto_fit suite of routines
;      work with CDS data, so I now check the INSTRUME tag when
;      computing the line width..
;    Ver.16, 9-Apr-2013, Peter Young
;      Corrected bug introduced with version 15.
;    Ver.17, 22-Jan-2014, Peter Young
;      If yshift_wvl and /map are set, then the routine now adjusts
;      the yc tag of the output structure to match the Y-cen value of
;      the reference wavelength, yshift_wvl.
;    Ver.18, 24-Feb-2014, Peter Young
;      Added the CALIB= optional input; missing data were not
;      correctly treated for some of the parameters when yshift_wvl
;      was specified, so this has been fixed (note that intensity was
;      OK). 
;    Ver.19, 29-Jul-2014, Peter Young
;      Now apply EIS-AIA offset to the output map coordinates.
;    Ver.20, 18-Aug-2014, Peter Young
;      ...but don't apply offset to IRIS data!
;    Ver.21, 22-Aug-2014, Peter Young
;      Caught a bug if the 'instrume' doesn't exist for fitdata.
;    Ver.22, 24-Sep-2014, Peter Young
;      Now doesn't require the line= input if /chi, /bg1 or
;      /bg2 are set. 
;    Ver.23, 2-Mar-2015, Peter Young
;      Modified the calls to 'reform' in order to return 2D arrays if
;      there is only one X element.
;    Ver.24, 03-Mar-2022, Peter Young
;      Added slit_ind to the map output structure.
;    Ver.25, 07-Oct-2022, Peter Young
;      Updated the numerical constant for computing the thermal width.
;      The previous constant had used the mass of the proton, when it
;      should be the atomic mass unit; added keyword /VWID for computing
;      the width in velocity units.
;    Ver.26, 26-mar-2025, Peter Young
;      Updated header; no change to code.
;    Ver.27, 09-May-2025, Peter Young
;      Implemented the Del Zanna et al. (2025) calibration.
;-

FUNCTION egf_make_map, fitdata, array, missing=missing, quiet=quiet,instrument=instrument

CASE instrument OF
  'IRIS': BEGIN
    id='IRIS'
    xy=[0.0,0.0]
  END 
  'CDS': BEGIN
    id='CDS'
    xy=[0.0,0.0]
  END 
  ELSE: BEGIN
    id='Hinode/EIS'
    xy=eis_aia_offsets(fitdata.date_obs)
    slit_ind=fix(fitdata.slit_ind)
    IF NOT keyword_set(quiet) THEN BEGIN
      print,'% EGF_MAKE_MAP:  The map coordinates have been corrected for the EIS-AIA spatial offset.'
    ENDIF 
  END 
ENDCASE 


IF tag_exist(fitdata,'xcen') THEN xc=fitdata.xcen+xy[0] ELSE xc=0.
IF tag_exist(fitdata,'ycen') THEN yc=fitdata.ycen+xy[1] ELSE yc=0.
IF tag_exist(fitdata,'scale') THEN BEGIN
  dx=fitdata.scale[0]
  dy=fitdata.scale[1]
ENDIF ELSE BEGIN
  IF fitdata.slit_ind EQ 2 THEN dx=2.0 ELSE dx=1.0
  dy=fitdata.ybin
ENDELSE 

outmap=make_map(array,xc=xc, yc=yc, dx=dx, dy=dy, time=fitdata.date_obs, id=id, $
                missing=missing, slit_ind=slit_ind)

return,outmap

END


;-----
FUNCTION egf_shift_array, int, err, wvl1, wvl2, missing=missing, outerr=outerr, ybin=ybin, $
                          delta_y=delta_y

IF n_elements(missing) EQ 0 THEN missing=-100.

siz=size(int)
nx=siz[1]
ny=siz[2]

int_new=fltarr(nx,ny)+missing
err_new=fltarr(nx,ny)+missing

off1=eis_ccd_offset(wvl1)
off1=off1[0]
off2=eis_ccd_offset(wvl2)
off2=off2[0]

delta_y=off2-off1

n=fix((off2-off1)/ybin)

CASE 1 OF
  n GT 0: BEGIN
    int_new[*,n:*]=int[*,0:ny-n-1]
    err_new[*,n:*]=err[*,0:ny-n-1]
  END
  n LT 0: BEGIN
    np=-n
    int_new[*,0:ny-np-1]=int[*,np:*]
    err_new[*,0:ny-np-1]=err[*,np:*]
  END 
  n EQ 0: BEGIN
    int_new=int
    err_new=err
  END 
ENDCASE

offset=(off2-off1)/ybin - n   ; note:  -1 < offset < 1

eis_fit_y_shift,int_new,err_new,offset,outint,outerr,missing=missing

return,outint

END


;-----------------------
FUNCTION eis_get_fitdata, fitdata, line=line, int=int, vel=vel, $
                          wid=wid, cen=cen, error=error, $
                          bg1=bg1, bg2=bg2, chi=chi, $
                          keep_inst_width=keep_inst_width, $
                          thermal_wid=thermal_wid, $
                          quiet=quiet, yshift_wvl=yshift_wvl, map=map, $
                          units=units, av_inst_width=av_inst_width, $
                          calib=calib, iris=iris, vwid=vwid


IF n_params() EQ 0 THEN BEGIN
  print,''
  print,'Use:  IDL> arr=eis_get_fitdata( fitdata [, /int, /vel, /wid, /cen, /bg1, /bg2, /chi,'
  print,'                                /keep_inst_width, thermal_wid=, /quiet, /map'
  print,'                                yshift_wvl=, units=, av_inst_width=, calib=)'
  print,''
  print,' Note: thermal_wid specified as [LOG_TMAX, MASS_NUM]'
  print,' Calib:  1 - Del Zanna (2013)'
  print,'         2 - Del Zanna (2013), with /undecay'
  print,'         3 - Warren et al. (2014)'
  print,'         4 - Warren et al. (2014), with /undecay'
  print,'         5 - Del Zanna et al. (2025)'
  print,'         6 - Del Zanna et al. (2025), with /undecay'
  print,''
  return,-1
ENDIF 


;
; Here I establish what instrument the fit structure belongs to. Note
; that if the 'instrume' doesn't exist then I assume it's an EIS file.
;
CASE 1 OF
  keyword_set(iris): instrument='IRIS'
  keyword_set(cds): instrument='CDS'
  ELSE: BEGIN
    IF tag_exist(fitdata,'instrume') EQ 1 THEN instrument=fitdata.instrume ELSE instrument='EIS'
  END
ENDCASE

;
; For these keywords the LINE= input is not important, so I have a
; flag to identify if one of these is set.
;
swtch=keyword_set(bg1) OR keyword_set(bg2) OR keyword_set(chi)

;
; The input LINE determines which of the emission lines (in the case
; of a multi-Gaussian fit) will be used to determine the orbit
; variation. The default is the first line.
;
ng=fitdata.ngauss
IF n_elements(line) EQ 0 AND ng GT 1 AND swtch EQ 0 THEN BEGIN
  print,''
  print,'%EIS_GET_FITDATA: The fit structure contains '+trim(ng)+' emission lines. Please select'
  print,'                  which line to use by specifying the keyword LINE=.'
  form1=trim(ng)+'i9'
  form2=trim(ng)+'f9.3'
  print,''
  print,format='("   LINE: ",'+form1+')',indgen(ng)
  print,format='(" REFWVL: ",'+form2+')',fitdata.refwvl
  print,''
  return,-1
ENDIF ELSE BEGIN
  IF n_elements(line) EQ 0 THEN line=0
ENDELSE 


refwvl=fitdata.refwvl[line]


;
; Intensity calibration adjustment (CALIB keyword).
;
IF n_elements(calib) NE 0 THEN BEGIN
  CASE calib OF 
    1: calib_factor=eis_ltds(fitdata.date_obs,refwvl)
    2: calib_factor=eis_ltds(fitdata.date_obs,refwvl,/undecay)
    3: calib_factor=eis_recalibrate_intensity(fitdata.date_obs,refwvl,1)
    4: calib_factor=eis_recalibrate_intensity(fitdata.date_obs,refwvl,1,/undecay)
    5: calib_factor=eis_recalibrate_intensity_new(fitdata.date_obs,refwvl,1)
    6: calib_factor=eis_recalibrate_intensity_new(fitdata.date_obs,refwvl,1,/undecay)
    ELSE: calib_factor=1.0
  ENDCASE 
 ;
 ; Output of eis_ltds is an array, so need to convert to scalar
  calib_factor=calib_factor[0]
ENDIF ELSE BEGIN
  calib_factor=1.0
ENDELSE 


th_wid=0.    ; set thermal width to zero by default
IF n_elements(thermal_wid) GT 0 THEN BEGIN
  IF n_elements(thermal_wid) NE 2 THEN BEGIN
    print,''
    print,'%EIS_GET_FITDATA: the optional input THERMAL_WID must be specified as a two element'
    print,'                  array of the form THERMAL_WID=[LOG10_TEMP,MASS_NUM] - see the routine'
    print,'                  header for more details.'
    print,'                  Returning...'
    return,-1
  ENDIF ELSE BEGIN
   ;
    th_wavel=fitdata.refwvl[line]
    th_temp=10.^float(thermal_wid[0])
    th_massnum=float(thermal_wid[1])
   ;
   ; compute thermal width
   ;  (7-Oct-2022: updated constant)
    th_wid=5.130e-13*th_temp*th_wavel^2/th_massnum
;    th_wid=5.093d-13*th_temp*th_wavel^2/th_massnum
    th_wid=sqrt(th_wid)
    IF NOT keyword_set(quiet) THEN BEGIN
      print,''
      print,'%EIS_GET_FITDATA: Subtracting thermal width from the measured width.'
     ;
      print,format='("   Thermal width (A): ",f7.4)',th_wid

      print,'Computed with the following parameters:'
      print,format='("   Wavelength (A): ",f10.3)',th_wavel
      print,format='("   Temperature (K): ",e10.2)',th_temp
      print,format='("   Mass number: ",i3)',th_massnum
    ENDIF 
  ENDELSE 
ENDIF 


;
; Background parameter 1
; ----------------------
IF n_elements(bg1) NE 0 THEN BEGIN
  out_array=reform(fitdata.aa[fitdata.ngauss*3,*,*],fitdata.nx,fitdata.ny)
  error=reform(fitdata.sigmaa[fitdata.ngauss*3,*,*],fitdata.nx,fitdata.ny)
 ;
 ; Flag missing data
  i=where(error EQ 0.,ni)
 ;
 ; Perform calibration correction (if necessary)
  out_array=out_array*calib_factor
  error=out_array*calib_factor
 ;
 ; Set missing values in arrays
  IF ni GT 0 THEN BEGIN
    out_array[i]=fitdata.missing
    error[i]=fitdata.missing
  ENDIF 
 ;
  IF n_elements(yshift_wvl) NE 0 THEN BEGIN
    out_array=egf_shift_array(out_array,error,refwvl,yshift_wvl,missing=fitdata.missing,outerr=error,ybin=fitdata.ybin,delta_y=delta_y)
  ENDIF 
 ;
  IF keyword_set(map) THEN output=egf_make_map(fitdata,out_array,missing=fitdata.missing,quiet=quiet,instrument=instrument) ELSE output=out_array
 ;
  units='erg cm^-2 s^-1 sr^-1 A^-1'
 ;
  return,output
ENDIF 


;
; Background parameter 2
; ----------------------
IF n_elements(bg2) NE 0 THEN BEGIN
  IF fitdata.nback EQ 2 THEN BEGIN
    out_array=reform(fitdata.aa[fitdata.ngauss*3+1,*,*],fitdata.nx,fitdata.ny)
    error=reform(fitdata.sigmaa[fitdata.ngauss*3+1,*,*],fitdata.nx,fitdata.ny)
   ;
   ; Flag missing data
    i=where(error EQ 0.,ni)
   ;
   ; Perform calibration correction (if necessary)
    out_array=out_array*calib_factor
    error=out_array*calib_factor
   ;
   ; Set missing values in arrays
    IF ni GT 0 THEN BEGIN
      out_array[i]=fitdata.missing
      error[i]=fitdata.missing
    ENDIF 
   ;
    IF n_elements(yshift_wvl) NE 0 THEN BEGIN
      out_array=egf_shift_array(out_array,error,refwvl,yshift_wvl,missing=fitdata.missing,outerr=error,ybin=fitdata.ybin,delta_y=delta_y)
    ENDIF 
   ;
    IF keyword_set(map) THEN output=egf_make_map(fitdata,out_array,missing=fitdata.missing,quiet=quiet,instrument=instrument) ELSE output=out_array
   ;
    units='erg cm^-2 s^-1 sr^-1 A^-1'
   ;
    return,output
  ENDIF ELSE BEGIN
    print,'%EIS_GET_FITDATA: fit only has one background parameter.'
    return,-1
  ENDELSE 
ENDIF


;
; Reduced Chi^2
; -------------
IF keyword_set(chi) THEN BEGIN
  out_array=fitdata.chi2
  IF n_elements(yshift_wvl) NE 0 THEN BEGIN
    error=fltarr(fitdata.nx,fitdata.ny)   ; use an empty error array for call to 
                                          ; egf_shift_array
    out_array=egf_shift_array(out_array,error,refwvl,yshift_wvl,missing=fitdata.missing,outerr=error,ybin=fitdata.ybin,delta_y=delta_y)
  ENDIF 
 ;
  IF keyword_set(map) THEN output=egf_make_map(fitdata,out_array,missing=fitdata.missing,quiet=quiet,instrument=instrument) ELSE output=out_array
 ;
  units='None'
 ;
  return,output
ENDIF 

;
; Velocity
; --------
IF keyword_set(vel) THEN BEGIN
  cen=reform(fitdata.aa[line*3+1,*,*],fitdata.nx,fitdata.ny)
  error=reform(fitdata.sigmaa[line*3+1,*,*],fitdata.nx,fitdata.ny)
  i=where(error EQ 0.,ni)     ; missing data
 ;
 ; convert to velocity
  out_array=lamb2v(cen-refwvl,refwvl)
  error=lamb2v(error,refwvl)
 ;
 ; set missing data 
  IF ni GT 0 THEN BEGIN
    out_array[i]=fitdata.missing
    error[i]=fitdata.missing
  ENDIF 
 ;
  IF n_elements(yshift_wvl) NE 0 THEN BEGIN
    out_array=egf_shift_array(out_array,error,refwvl,yshift_wvl,missing=fitdata.missing,outerr=error,ybin=fitdata.ybin,delta_y=delta_y)
  ENDIF 
 ;
  IF keyword_set(map) THEN output=egf_make_map(fitdata,out_array,missing=fitdata.missing,quiet=quiet,instrument=instrument) ELSE output=out_array
 ;
  units='km s^-1'
 ; 
  return,output
ENDIF


;
; Width
; -----
; The width stored in the FITDATA structure is the Gaussian width. I
; convert this to FWHM by multiplying by 2*sqrt(2*ln(2))=2.35.
;
; The instrumental width is automatically removed from the measured
; line widths (use /KEEP_INST_WIDTH to switch this off). 
;
; When subtracting the instrumental and/or thermal widths, the
; width may end up disappearing (see the 'chck' and 'ichck' parameters
; below). In this case I just set the output width to zero, but
; keep the error as it is. I do not flag these pixels as missing. 
;
; I do not correct the error when performing the instrumental and
; thermal width corrections.
;

IF keyword_set(vwid) THEN wid=1b
;
IF keyword_set(wid) THEN BEGIN
  out_array=reform(fitdata.aa[line*3+2,*,*],fitdata.nx,fitdata.ny)*2.*sqrt(2.*alog(2.))
  error=reform(fitdata.sigmaa[line*3+2,*,*],fitdata.nx,fitdata.ny)*2.*sqrt(2.*alog(2.))
 ;
 ; Get instrumental width for the Y-positions of FITDATA
 ;
  inst_width_arr=fltarr(fitdata.nx,fitdata.ny)
  ybin=round(fitdata.ybin)
  yip=fitdata.yip + indgen(fitdata.ny)*ybin + ybin/2
 ;
 ; 4-Apr-2013. I've added the check below in case CDS data are being used.
 ;
  IF instrument EQ 'EIS' THEN BEGIN 
    inst_width=eis_slit_width(yip,slit_ind=fitdata.slit_ind)
  ENDIF ELSE BEGIN
    inst_width=0.
  ENDELSE 
 ;
  FOR j=0,fitdata.nx-1 DO inst_width_arr[j,*]=inst_width
 ;
  i=where(error EQ 0.,ni)  ; these are missing data points
  j=where(error NE 0.,nj)
 ;
 ; This gives the average instrumental width for the pixels where the
 ; fit is good
 ;
  IF  nj NE 0 THEN BEGIN
    av_inst_width=mean(inst_width_arr[j])
  ENDIF 
 ;
  IF keyword_set(keep_inst_width) THEN BEGIN 
    chck=out_array[j]^2  - th_wid[j]^2
  ENDIF ELSE BEGIN
    chck=out_array[j]^2 - inst_width_arr[j]^2 - th_wid[j]^2
  ENDELSE 
 ;
  ichck=where(chck LT 0.,nchck)
  perc_chck=float(nchck)/float(n_elements(out_array))*100.
  perc_chck_str=trim(string(format='(f10.2,"%")',perc_chck))
  IF nchck GE 1 AND NOT keyword_set(quiet) THEN BEGIN
    print,''
    print,'WARNING: the instrumental width and/or thermal width is larger than the measured width '
    print,'for '+trim(nchck)+' pixels ('+perc_chck_str+') . The widths for these pixels have been set to zero.'
  ENDIF
  chck=chck>0.   ; Force negative chck values to be zero.
 ;
 ; If /vwid was set, then the FWHM in angstroms is converted to the 1/e
 ; width in km/s units. This is the standard non-thermal velocity that
 ; is usually assigned the Greek letter xi. See for example, Equation 1 in
 ; Young et al. (2013, ApJ, 766, 127).
 ;
  IF keyword_set(vwid) THEN BEGIN
    out_wid=sqrt(chck)
    out_wid=out_wid/5.554e-6/fitdata.refwvl[line]
  ENDIF ELSE BEGIN
    out_wid=sqrt(chck)
  ENDELSE 
 ;
  IF nj NE 0 THEN out_array[j]=out_wid
  IF ni NE 0 THEN BEGIN
    out_array[i]=fitdata.missing
    error[i]=fitdata.missing
  ENDIF 
 ;
  IF n_elements(yshift_wvl) NE 0 THEN BEGIN
    out_array=egf_shift_array(out_array,error,refwvl,yshift_wvl,missing=fitdata.missing,outerr=error,ybin=fitdata.ybin,delta_y=delta_y)
  ENDIF 
 ;
  IF keyword_set(map) THEN output=egf_make_map(fitdata,out_array,missing=fitdata.missing,quiet=quiet,instrument=instrument) ELSE output=out_array
 ;
  IF keyword_set(vwid) THEN units='km s^-1' ELSE units='A'
 ;
  return,output
ENDIF 


;
; Centroid
; --------
IF keyword_set(cen) THEN BEGIN
  out_array=reform(fitdata.aa[line*3+1,*,*],fitdata.nx,fitdata.ny)
  error=reform(fitdata.sigmaa[line*3+1,*,*],fitdata.nx,fitdata.ny)
 ;
  i=where(error EQ 0.,ni)     ; missing data
  IF ni GT 0 THEN BEGIN
    out_array[i]=fitdata.missing
    error[i]=fitdata.missing
  ENDIF 
 ;
  IF n_elements(yshift_wvl) NE 0 THEN BEGIN
    out_array=egf_shift_array(out_array,error,refwvl,yshift_wvl,missing=fitdata.missing,outerr=error,ybin=fitdata.ybin,delta_y=delta_y)
  ENDIF 
 ;
  IF keyword_set(map) THEN output=egf_make_map(fitdata,out_array,missing=fitdata.missing,quiet=quiet,instrument=instrument) ELSE output=out_array
 ;
  units='A'
 ;
  return,output
ENDIF 


;
; Intensity
; ---------
; Note that the fit.int and fit.interr tags have missing values set to fitdata.missing
;
out_array=reform(fitdata.int[line,*,*],fitdata.nx,fitdata.ny)
error=reform(fitdata.interr[line,*,*],fitdata.nx,fitdata.ny)

; Perform calibration correction (if necessary)
i=where(error NE fitdata.missing,ni)
out_array[i]=out_array[i]*calib_factor
error[i]=error[i]*calib_factor

IF n_elements(yshift_wvl) NE 0 THEN BEGIN
  out_array=egf_shift_array(out_array,error,refwvl,yshift_wvl,missing=fitdata.missing,outerr=error,ybin=fitdata.ybin,delta_y=delta_y)
ENDIF 

IF keyword_set(map) THEN output=egf_make_map(fitdata,out_array,missing=fitdata.missing,quiet=quiet,instrument=instrument) ELSE output=out_array

units='erg cm^-2 s^-1 sr^-1'

IF keyword_set(map) AND n_elements(yshift_wvl) NE 0 THEN output.yc=output.yc-delta_y

return,output


END