pro dmg_flat, files,  vflat, flati, darkfile, x0, y0,  subpixel=subpixel,   $
                 ndark_integ=ndark_integ, ndata_integ=ndata_integ , shift=shift
;+
; NAME: DMG_FLAT
; PURPOSE: Produce the flat patterns for DMG observations using a new algorithm
; CALLING SEQUENCE:
;         DMG_flat, files, vflat, flati, darkfile
; INPUT:
;      files        an array of  polarization data files 
;
; OPTIONAL INPUT:
;      darkfile    an polarization dark frame file 
; INPUT KEYWORDS:
;   subpixel    if set, the determnaition is done with subpixel-accuracy
;   shift          if set, the given shift is used.
; OUTPUT:
;      vflat        a V flat image
;      flati        an intensity flat image     
; MODIFICATIONS:
;   J. Varsik   1999-11-10 Changed I/O to be similar to vmg_flat.  Added
;                          intensity flat field.
;-
          
          ; Read Dark Frames

if n_elements(ndata_integ) eq 0 then ndata_integ=1.
if n_elements(ndark_integ) eq 0 then ndark_integ=1. 
if n_elements(darkfile) eq 0 then begin
   dark=0.
endif else begin
 dark=float(readfits(darkfile, /sil))/ndark_integ
   
endelse

          ; Read Data and subtract dark levels
 
nf = n_elements(files)
ie= fltarr(512, 512, nf)
io= fltarr(512, 512, nf)
for k=0, nf-1 do begin
    tmp=(float(readfits(files(k), /sil))/ndata_integ -dark)
    ie(*,*,k) = tmp(*,*,0)
    io(*,*,k) = tmp(*,*,1)
endfor
        
           ; Determine Relative Offsets
if not keyword_set(shift) or n_elements(x0) ne nf then begin  
x0=fltarr(nf)
y0=fltarr(nf)
ss=nf/2

reference = (ie(*,*,ss)-io(*,*,ss))/(ie(*,*,ss)+io(*,*,ss))
w=1.; bell_window(reference, 0.1)

for k=0, nf-1 do begin  
sh =  alignoffset(w*(ie(*,*,k)-io(*,*,k))/(ie(*,*,k)+io(*,*,k)), w*reference)
x0(k) = sh(0)
y0(k) = sh(1)
print, sh
endfor   
endif         
            ; Determine Flat Images

;flate = gaincalib(ie, x0, y0, maxiter=10)
;flato = gaincalib(io, x0, y0,  maxiter=10)
flati = gaincalib(ie+io, x0, y0, maxiter=10)
print,'flati computed'
flatv = gaincalib((ie-io)/(ie+io)>(-0.02)<0.02, x0, y0, maxiter=10, additive=1)
print,'flatv computed'
flate=flati*(double(flatv)+1.d)*0.5d
flato=flati*(-double(flatv)+1.d )*0.5d
flato = float(flato*median(io)/median(ie))
flate = float(flate)

            ; Save Flat Images

; apply calibration scaling
;vflat = 14500.0 * [[[flate]], [[flato]]]
vflat = 14500.0 * flatv
                          

end