function gaincalib, images, x0, y0, additive=additive, shift=shift, adj=adj 
;+
;  NAME:  GAINCALIB
;  PURPOSE:    Produce either a gain table or an offset table
;              from a set of images
;              taken shifted with respect to one another
;  CALIING SEQUENCE:
;          Result = gaincalib(images, additive=additive)
; INPUT:
;          images    a three-dimensional array representing
;                    a sequence of two-dimensional images
;                    images(*,*,k) ( k=0, 1,.., N-1).
;                    It is strongly recommended that the number
;                    of elements in each image is a power of 2
;                    for the efficent FFT processing.
; OPTIONAL INPUT:
;          x0        an array of x-shift
;          y0        an array of y-shift
;                        
; KEYWORD:  
;         additive   If set,  the program produces
;                    the offset table rather than the gain table.
;                    If not set, the program produces the gain
;                    table. In this case, the input images should
;                    be positive-definite.
;         shift     if set, the program calculates the amount of shift
;                   in each image. 
;         adj       if set, use only the pairs two of which were $
;                   taken succesively
; OUTPUT:
;        Result     the gain table if the keyword  ADDITIVE is not set
;                   or the offset table if the keyword is set.
; OPTIONAL OUTPUT:
;          x0        an array of x-shift
;          y0        an array of y-shift
; COMMENTS:
;           This uses a Fouier method to determine the flat-pattern of
;           of a detector. Therefore it works best when the following
;           conditions are satisfied:
;            
;               1) negligible low frequency gain pattern. 
;
;                  If this pattern is non-negligible, 
;                  then obtain  first  the low
;                  frequency pattern based on surface fittings, and 
;                  devide (or  subtract in the case of additive pattern)
;                  the raw image by this low frequency pattern.
;                  Then it can be multiplied back to the output of
;                  this routine to yield the practical gain pattern.
;                
;              2) pixel number which is equal to a power of 2 for
;                  the efficient FFT.    
;-                                              

if keyword_set(additive) then additive=1 else additive=0
s=size(images)
x = indgen(s(1))#replicate(1, s(2))
x = x - s(1)*(x ge s(1)/2)
kx = 2.*!pi*x/s(1)
y = replicate(1, s(1))#indgen(s(2))
y = y - s(2)*(y ge s(2)/2)
ky=2.*!pi*y/s(2)
cimages = complexarr(s(1), s(2), s(3))

for k=0, s(3)-1 do if additive $
    then cimages(*,*,k) = fft(images(*,*,k),1) $
    else cimages(*,*,k) = fft(alog(images(*,*,k)),1)
if n_elements(x0) eq 0 then shift=1
if keyword_set(shift) then begin
x0 =fltarr(s(3))
y0 = fltarr(s(3))
w=bell_window(images(*,*,0), 0.1)
cref = fft((images(*,*,0)-median(images(*,*,0)))*w, -1)
for k=1, s(3)-1 do begin
cor = float(fft(fft((images(*,*,k)-median(images(*,*,k)))*w,1)*cref -1))
correl=max(cor)
ss=where(cor ge 1.0*correl)
x0(k)  = total(x(ss)*cor(ss))/total(cor(ss))
y0(k)  = total(y(ss)*cor(ss))/total(cor(ss))
print, k, x0(k), y0(k)
endfor
endif 

nom =0.
den =0.

if keyword_set(adj) then adj=1 else adj=0

for i=1, s(3)-1 do for j=0+adj*(i-1), i-1 do  begin

aa =  exp(complex(0., -(x0(i)*kx+y0(i)*ky)))
bb =  exp(complex(0., -(x0(j)*kx+y0(j)*ky)))
den = den + abs(aa-bb)^2
nom = nom + (cimages(*,*,i)*aa - cimages(*,*,j)*bb)*conj(aa-bb)
end
verysmall = den(1, 0)*0.01
den = den > verysmall

result = float(fft(nom/den, -1))
if not additive then result=exp(result) else result=result
return, result
end