function sc_inverse,n,diag, below, above
;+
;$Id: sc_inverse.pro,v 1.3 2007/05/18 11:13:45 crothers Exp $
;
; Project   : STEREO SECCHI
;
; Name      : sc_inverse
;
; Purpose   : Computes the inverse of a matrix of size n where the major diagonal
;             contains the value diag, the lower triangle contains below and the
;             upper triangle contains the value above.
;                i.e. X[i,i]=diag, X[0:i-1,i]=below, X[0,i+1:n-1]=above

; Explanation: This kind of matrix is generated by an exposure on a CCD without
;              a shutter where diag is the time the image is still on the CCD
;              below is the time between line shifts in the clear of the CCD and
;              above is the time between line shitfs in the readout.  The
;              observed values can be matrix multiplied to remove the exposure
;              properties from the image.
;
; Use      : IDL> correction=sc_inverse(nrows,still_time,downtime,uptime)
;
; Inputs   : nrows - the number of rows/columns in the output matrix
;            still_time - the time the image is static
;            downtime - the shift time down the image
;            uptime - the shift time up the image
;
; Outputs   : correction - the matrix to perform correction
;
; Keywords  :
;
; Common    :
;
; Calls     :
;
; Category  : Calibration
;
; Prev. Hist. : None.
;
; Written     : Steve Crothers, Chris Johnson, RAL
;
;$Log: sc_inverse.pro,v $
;Revision 1.3  2007/05/18 11:13:45  crothers
;Commented routine to show intended purpose and removed "correction" that altered this.
;
;
; Revision 1.2  2007/03/30 00:00:00  newmark
;    fixed bug to not divide by exposure time
;
;
; Revision 1.1  2007/01/22 19:35:08  colaninn
; new from RAL

compile_opt idl2

wt_above=double(above)/diag
wt_below=double(below)/diag

wt_above_1=wt_above-1
wt_below_1=wt_below-1

power_above=dblarr(n-1)
power_below=dblarr(n-1)

power_above[0]=1
power_below[0]=1

for row=1,n-2 do begin
    power_above[row]=power_above[row-1]*wt_above_1
    power_below[row]=power_below[row-1]*wt_below_1
end

v= [ 0 , wt_below *(power_below*reverse(power_above))]
u= [ 0 , wt_above *(power_above*reverse(power_below))]

d = -u[1]/wt_above - (total(v)-v[N-1])
f = 1/(diag*(d+wt_above*total(v)))
;f = 1/(d+wt_above*total(v))     ;bug fix - no change of function do this elsewhere

u[0]=d
v[0]=d
u=u*f
v=reverse(v)*f
p=dblarr(n,n,/nozero)
p[0,0]=u
for row=1,n-2 do begin
p[0,row]=v[n-row-1:n-2]
p[row,row]=u[0:n-row-1]
end
p[0,n-1]=v

return,p

end