function eis_int_err_fit, windata, plot=plot, quiet=quiet, g=g, h=h, $
                          whole_array=whole_array, $
                          noabs=noabs, photon=photon

;+
; NAME:
;     EIS_INT_ERR_FIT
;
; PURPOSE:
;     Returns the coefficients F for the straight line fit that relates
;     the function H to the function G (defined below). These
;     coefficients are are used for estimating uncertainties for
;     interpolated EIS data. Further details are given in EIS Software
;     Note No.13.
;
; CATEGORY:
;     Hinode; EIS; interpolation; errors.
;
; CALLING SEQUENCE:
;     Result = EIS_INT_ERR_FIT( Windata )
;
; INPUTS:
;     Windata:  The structure output by EIS_GETWINDATA.
;
; KEYWORD PARAMETERS:
;     PLOT:    If set, then two plots are created showing the fit. The
;              left-hand plot shows all the data, while the right-hand
;              plot shows the same data but with the X-range set to a
;              factor 10 smaller. Dots are used to mark the data
;              points, and a thick line is used to show the fit.
;
;     WHOLE_ARRAY: If set then the whole array is processed even if it
;                  has more than 100,000 elements.
;
; OUTPUTS:
;     A two element array describing the straight line fit to G and H,
;     i.e.,
;
;      H = F(1) * G + F(0)
;
; OPTIONAL OUTPUTS:
;     G:  An array of same size as WINDATA.INT containing the
;         calculated G values.
;
;     H:  An array of same size as WINDATA.INT containing the
;         calculated H values.
;
;     NOABS: Set to 1 if the data were calibrated with the /noabs option.
;            This is determined from the windata.hdr.calstat.
;
;     PHOTON: Set to 1 if the data were calibrated with the /photon option.
;             This is determined from the windata.hdr.calstat.
;
; PROGRAMMING NOTES:
;     For data in calibrated units of erg/cm2/s/sr/A, the G and H
;     functions are defined as:
;
;      G = INT * EFF_AREA / WVL
;      H = ERR^2 * EFF_AREA^2
;
;     These definitions mean that there is a linear relationship between
;     G and H such that H = G * F(1) + F(0). The F parameters are derived
;     by fitting a straight line to G and H for the good data in the EIS
;     dataset. The parameters can then be used to infer uncertainties
;     for interpolated EIS data.
;
;     For data in DN units (when the /noabs keyword is given to eis_prep),
;     G and H are defined as:
;
;      G = INT / WVL
;      H = ERR^2
;
;     For data in photon units (when /photon keyword given to eis_prep),
;     G and H are defined as
;
;      G = INT / WVL^2
;      H = ERR^2 / WVL^2
;
; EXAMPLE:
;     IDL> wd=eis_getwindata(file,195.12)
;     IDL> f=eis_int_err_fit(wd)
;
; MODIFICATION HISTORY:
;     Ver.1, 16-Apr-2009, Peter Young
;     Ver.2, 15-Jan-2010, Peter Young
;        tried to improve memory usage
;     Ver.3, 15-Apr-2010, Peter Young
;        added check in case interpolation method has not worked well
;     Ver.4, 16-Apr-2010, Peter Young
;        to improve memory usage the routine now checks the size of
;        the input array; if too big then only a sub-array centered on
;        the brightest pixel is processed. (Note: this behavior is
;        controlled by the keyword WHOLE_ARRAY.)
;     Ver.5, 13-Apr-2011, Peter Young
;        I now check for Infs and NaNs in the h array and exclude
;        these from the line fit; added status= output
;     Ver.6, 13-Nov-2012, Peter Young
;        Expanded text of error message; no change to code.
;     Ver.7, 3-Mar-2016, Peter Young
;        Updated my e-mail address in the error message. No change to
;        code.
;     Ver.8, 13-May-2025, Peter Young
;        Added /noabs keyword; updated my email
;        address.
;     Ver.9, 06-Aug-2025, Peter Young
;        Changed expressions for G and H after finding problems with
;        spectral atlas data with the /noabs keyword; implemented g and
;        h definitions for /photon data ; remove status=
;        output as it was not being used; updated header; g and h
;        pixels not used in the fit are now set to missing.
;-

int=windata.int
err=windata.err
wvl=windata.wvl

;
; The following checks the size of the data array. If it's too
; big, then only a sub-array is processed.
;
IF n_elements(int) GT 100000 AND NOT keyword_set(whole_array) THEN BEGIN
  getmax=max(int,imax)
  ind=array_indices(int,imax)
  siz=size(int)
  xmin=max([0,ind[1]-50])
  xmax=min([siz[2]-1,ind[1]+50])
  ymin=max([0,ind[2]-50])
  ymax=min([siz[3]-1,ind[2]+50])
  int=int[*,xmin:xmax,ymin:ymax]
  err=err[*,xmin:xmax,ymin:ymax]
 ;
  IF NOT keyword_set(quiet) THEN BEGIN
    message,/info,/cont,'Interpolation performed on sub-array [*,'+trim(xmin)+':'+trim(xmax)+','+trim(ymin)+':'+trim(ymax)+']'
  ENDIF 
ENDIF


siz=size(int)
n1=siz[1]
n2=siz[2]
n3=siz[3]
IF siz[0] EQ 4 THEN n4=siz[4]

ea=eis_eff_area(wvl,/quiet)

noabs=0b
chck=strpos(windata.hdr.calstat,'ABS')
IF chck LT 0 THEN BEGIN
  noabs=1b
  IF NOT keyword_set(quiet) THEN BEGIN 

    message,/info,/cont,'ABS was not found in windata.hdr.calstat. The NOABS option has been set. Set NOABS=0 to switch this off.'
  ENDIF 
ENDIF 

photon=0b
chck=strpos(windata.hdr.calstat,'PHOT')
IF chck GE 0 THEN BEGIN
  photon=1b
  IF NOT keyword_set(quiet) THEN BEGIN 
    message,/info,/cont,'PHOT has been found in windata.hdr.calstat, so assuming units of phot/cm2/s/sr/Ang.'
    ENDIF 
ENDIF

;
; Create G and H arrays
;
g=int  ; Set g and h to have same size as int; values will be overwritten
h=int  ;
FOR ii=0,n1-1 DO BEGIN
  CASE siz[0] OF
    3: BEGIN
      CASE 1 OF
        keyword_set(noabs): BEGIN
          g[ii,*,*]=int[ii,*,*]/wvl[ii]
          h[ii,*,*]=err[ii,*,*]^2
        END
        ;
        keyword_set(photon): BEGIN
          g[ii,*,*]=int[ii,*,*]/wvl[ii]^2
          h[ii,*,*]=err[ii,*,*]^2/wvl[ii]^2
        END
        ;
        ELSE: BEGIN
          g[ii,*,*]=int[ii,*,*]*ea[ii]/wvl[ii]
          h[ii,*,*]=err[ii,*,*]^2*ea[ii]^2
        END 
      ENDCASE 
    END
    ;
    4: BEGIN
      CASE 1 OF
        keyword_set(noabs): BEGIN
          g[ii,*,*,*]=int[ii,*,*,*]/wvl[ii]
          h[ii,*,*,*]=err[ii,*,*,*]^2
        END
        ;
        keyword_set(photon): BEGIN
          g[ii,*,*,*]=int[ii,*,*,*]/wvl[ii]^2
          h[ii,*,*,*]=err[ii,*,*,*]^2/wvl[ii]^2
        END
        ;
        ELSE: BEGIN
          g[ii,*,*,*]=int[ii,*,*,*]*ea[ii]/wvl[ii]
          h[ii,*,*,*]=err[ii,*,*,*]^2*ea[ii]^2
        END 
      ENDCASE 
    END
  ENDCASE
ENDFOR

;
; Very large values of err (i.e., > 1e19) will generate Infs when
; squared. The following checks for Infs (or NaNs).
; 
finite_chck=finite(h)

;
; Only perform fit to non-missing data, those intensities greater
; than zero, and finite values of h.
; 
ind_fit=where(err NE windata.missing AND int GT 0. AND finite_chck EQ 1,n_ind)

err=0  ; free up disk space
int=0  ;

f=poly_fit(g[ind_fit],h[ind_fit],1)

hx=f[0]+f[1]*g[ind_fit]

;
; The following checks the number of points that are more than 20%
; from the fit. If the number exceeds 1% a warning message is
; produced. 
;
chck=where((abs(hx-h[ind_fit])/h[ind_fit]) GT 0.20,nchck)
chck=float(nchck)/float(n_ind)
IF chck GT 0.01 THEN BEGIN
  print,''
  print,'% EIS_INT_ERR_FIT:  ***PROBLEM WITH /REPLACE METHOD***'
  print,''
  print,format='("Percentage of points more than 20% from fit:",f5.1,"%   (> 1%)")',chck*100.
  print,''
  print,'The method used to interpolate intensities and errors for missing pixels has '
  print,'not worked well for this data-set, and so the data should not be used.'
  print,'Please run the routine again without the /replace (or /refill) keyword.'
  print,''
  print,'Please also send an e-mail to peter.r.young@nasa.gov giving the filename and'
  print,'the calling sequence you used so that the problem can be diagnosed.'
  print,''
ENDIF ELSE BEGIN
  IF NOT keyword_set(quiet) THEN BEGIN 
    print,format='("Percentage of points more than 20% from fit:",f6.2,"%   (< 1%)")',chck*100.
  ENDIF 
ENDELSE 

IF NOT keyword_set(quiet) THEN BEGIN 
  print,format='("Percentage of pixels used in fit: ",f5.1," %")', $
        float(n_ind)/float(n_elements(h))*100.
ENDIF

IF keyword_set(plot) THEN BEGIN
  !p.multi=[0,2,1]
  xmin=min(g[ind_fit])
  xmax=max(g[ind_fit])
 ;
  plot,psym=3,g[ind_fit],h[ind_fit], xtit='G', ytit='H', $
       /xsty,charsiz=2
  x=findgen(101)/100.*(xmax-xmin)+xmin
  oplot,th=2,x,f[1]*x+f[0]
 ;
  plot,psym=3,g[ind_fit],h[ind_fit], xtit='G', ytit='H', $
       xra=[xmin,xmax/10.],charsiz=2
  x=findgen(101)/100.*(xmax/10.-xmin)+xmin
  oplot,th=2,x,f[1]*x+f[0]
  
  !p.multi=0
ENDIF

;
; Set pixels not used in fit to be missing.
;
mask=bytarr(size(g,/dim))
mask[ind_fit]=1b
k=where(mask EQ 0,nk)
IF nk NE 0 THEN BEGIN
  g[k]=windata.missing
  h[k]=windata.missing
ENDIF 

return,f

END