FUNCTION lsback_away,image_in,index_in,model

; =========================================================================
;+
; PROJECT:
;       Solar-B / XRT
;
; NAME:
;       LSBACK_AWAY   3f
;  (contains ests of T dep of base, elims separate orb  dep)
;  (outputs just model)
;
; CATEGORY:
;       Data calibration
;
; PURPOSE:
;       Remove Large-scale background (pedestal+dark+readnoise) 
;       introduced to the data by thermal noise and read-out electronics.
;
; CALLING SEQUENCE:
;       image_out = LSBACK_AWAY (image_in , index_in, [,model]) 
;
; INPUTS:
;       IMAGE_IN  - [Mandatory] (2-dim float array, [Nx,Ny])
;                   The image containing read-out signals.
;       INDEX_IN  - [Mandatory] (Index structure corresponding to image_in)
;
; OUTPUTS:
;       IMAGE_OUT - [Mandatory] (2-dim float array, [Nx,Ny])
;                   The image with the fit to read-out signals removed.
;       MODEL     - [Optional] (2-dim float array, [Nx,Ny])
;                   The model of the large scale background
;
; EXAMPLES:
;
;       Basic usage:
;       IDL> image_out = lsback_away(image_in, index_in)
;
;       You want the dark model for some reason:
;       IDL> image_out = lsback_away(image_in, index_in, model)
;
;
;  COMMON BLOCKS:
;       None.
;
; NOTES:
;
;       1) a set of 2129 darks of sizes >=2048^2 (~97.8% of those available
;          up to 2/08) processed by this routine produced an average result 
;          <dark> = -0.005 +/- 0.231 DN.
;       2) parameters of the fit may vary in time and need updating
;          in the future
;       3) dependence on binning, exposure time, and compression (none?)
;          are less well deteremined
;       4) Does no error checking! (this is done by xrt_clean_ro 
;          which calls it)
;
; MODIFICATION HISTORY:
;
progver = 'v2007.Mar.06' ; --- (SSaar, MWeber) Written by SSaar. Cleaned 
;                               up by MWeber.
progver = 'v2008.May.18' ; --- (SSaar) Fixed bug relating to binned data.
;                               Includes empirical estimates of the effects
;                               of CCD temperature, orbital variation and 
;                               exposure time.
progver = 'v2008.Jun.18' ; --- (SSaar) Merged of CCD temperature and orbital 
;                               variation in a revised calibration.  
;                  
;
;
;-
; =========================================================================

;===== Get dimensions and parameters of image.
  nx = index_in.naxis1
  ny = index_in.naxis2
  chip_sum = index_in.chip_sum   
  t_exp = index_in.exccdex/1e6  
  tmx = 1677721500L
  if t_exp lt 0 then t_exp =(tmx -index_in.e_sh_ope +index_in.e_sh_clo)/1e6
  t_ccd = index_in.ccd_tmpc      
  time0 = anytim2jd(index_in.date_obs)
  t0i = 2454035d0                       ;  JD time of 1st dark image
  time = time0.int - t0i + time0.frac

;===== use best fits from 2048^2 darks to compute ramp parameters 
;=====  amplitude of exponential 
  c_amp = [4.01,5,4.29,6.6]
  tx = alog10(t_exp)+6.
  a_slope = (c_amp(2)-c_amp(0))/(c_amp(3)-c_amp(1))
  amp =  c_amp(0)*(tx lt c_amp(1)) + $
        (a_slope*(tx-c_amp(1))<(c_amp(2)-c_amp(0)) + c_amp(0))*(tx ge c_amp(1)) 
;=====  e-folding width of exponential
  wid = 188.199 - 8.42898*chip_sum  
;=====  slope of ramp 
  c_slope=[0.0004556465,2.522937e-06]
  slope = poly(t_ccd,c_slope)

;===== construct 1-D ramp in y
  y = findgen(ny)
  ramp = amp*exp(-y/wid) + slope*y 

;===== orbital dependence of base level
;  c_orb=[0.067944796,0.70718654,1.071]
;  a_orb=10.^c_orb(0)*chip_sum^c_orb(1) - c_orb(2) ; orb term amp (sig=0.03264)
;  p_orb=0.0684052                ; orbit in days  
;  phi_orb = 0.025                ; phase shift of orbital term
;  phi_obs=time/p_orb mod 1.
;  orb_term  = a_orb*sin(2*!pi*(phi_obs + phi_orb))
   orb_term = 0.

;===== compute base level 
;===== dependence on binning
  c_bin=[2.0165648,1.0042904,62.597]
  back0=10.^c_bin(0)*chip_sum^c_bin(1) - c_bin(2)   ; binning term (sig=0.1603)

;===== dependence on CCD temperature
  c0=44.502073d0
  c_tccd1=[c0,0.16950656,0.0019546097]   ; predicts 41.7338 @ -68.8979
  c_tccd2=[c0,2.4585699,0.023489056]  ; sig=0.53  c0=207.93645
  c_tccd4=[c0,2.4585699*2*0.9,0.023489056*2*0.81] ; scaled (from 8 refit) 
                                                  ; extrap of bin=2
  c_tccd8=[c0,2.4585699*4-0.93596664,0.023489056*4-0.017488685] ; refit of 
                                                            ; extrap from bin=2
                                                            ; sig=2.9836
  case 1 of 
      (chip_sum eq 1): c_tccd=c_tccd1
      (chip_sum eq 2): c_tccd=c_tccd2
      (chip_sum eq 4): c_tccd=c_tccd4
      else: c_tccd=c_tccd8
  end
;===== dependence on exposure time
  time_term =  .04454/31*chip_sum^2*t_exp
  case 1 of
     (chip_sum eq 1): offset = -0.024      ; small correction to overall fit
     (chip_sum eq 2): offset = 58.111956-0.034165056 ; correction to overall fit
     (chip_sum eq 4): offset = 119.715    ; correction to overall fit
     else: offset = 247.79423+3.2745614   ; correction to overall fit
  end
  base = poly(t_ccd,c_tccd) + offset - (10.^c_bin(0) - c_bin(2)) 

;===== add up base level terms, construct full 1-D background
  back1 = back0 +  base + time_term + orb_term
  yback = ramp + back1

;===== Make 2D.  Rotate into y direction
  back = rotate(rebin(yback,ny,nx),1)         

;===== Subtract, and you're done!
  image_out = image_in - back 
  model = back
  
return, image_out
 
end
;
;
;