;------------------------------------------------------------------
;+
; NAME:
;	LCSETUPFF
; PURPOSE:
;	This procedure accepts two BBSO VMG FITS files, one V
;	image (magnetogram) and one J or I image (direct image).
;	A flat-field magnetogram generated from Jongchul Chae's
;	vmg_flat procedure is subtracted from the magentogram
;	to remove large-scale non-uniformities in the magnetogram.
;
;	Then, both the V and the J (or I) images are rebinned from
;	512 x 481 to 512 x 384 to make the pixels square.
;	The resulting images are returned in the viss and jiss
;       variables.
;
;	Note that the old IRAF version deleted the original FITS
;	files.  This version does not.
;
;	Note that the old version (using the boxcar smoothing)
;	did not scale the VMG image when it was loaded.  This
;	version does.
;
;       This version is for use with old Datacube files without
;       flatfield frames.  A pseudo-flatfield is made using
;       boxcar smoothing.
;
;       This version also makes a fake flatfield for the J image
;       and applies it.
; CATEGORY:
;	POLAR
; CALLING SEQUENCE:
;	LCSETUPFF,vfile,jfile,vih,jih,viss,jiss
; INPUTS:
;	vfile = file name for input magnetogram (FITS format).
;	jfile = file name for input direct image (FITS format).
;       imscl = image scale (pixels per arcsec).
; KEYWORD PARAMETERS:
; OUTPUTS:
;	viss = IDL variable containing rebinned and "flattened" magnetogram
;	jiss = IDL variable containing rebinned, flattened direct image
; COMMON BLOCKS:
; NOTES:
; MODIFICATION HISTORY:
;	J. Varsik,	13 Oct. 1995 --- Port from IRAF task
;	J. Varsik,	08 Jan. 1999 --- Change to use flatfield image.
;       J. Varsik,      24 Feb. 1999 --- Add J flatfield.
;                                        J. Varsik,      03 May  1999
;                                        --- Modified from lcsetup,
;                                            remove flats.
;
;-
;-------------------------------------------------------------------

function lcsetupff,vfinam,jfinam,vih,jih,viss,jiss,imscl

    ; Display IDL header if help is required.
    if (keyword_set(help)) then begin
      get_idlhdr,'lcsetup.pro'
      goto,finishup
    endif

    ; Read V file
    vi = READFITS(vfinam,vih)
    vixa = make_array(640,640,/float,value=0.0)
    vixa(64,64) = vi
    ; Generate boxcar smoothed file and subtract.
    s = SMOOTH(vixa,50)
    vo = vixa - s
    vos = vo(64:575,64:544)
    so = s(64:575,64:544)

    ; Use this (vos) as initial guess
    vbig = where(abs(vos) gt 10.0)
    vib = vi
    vib(vbig) = so(vbig) ;clip large magnetic field elements

    ; Now regenerate boxcar smoothing and subtract
    vixb = make_array(640,640,/float,value=0.0)
    vixb(64,64) = vib
    ; Generate boxcar smoothed file and subtract.
    s = SMOOTH(vixb,50)
    vo = vixa - s
    vos = vo(64:575,64:544)
    

    ; Read J file
    ji = READFITS(jfinam,jih,/NOSCALE)

    ; Construct fake solar image
    ew = fxpar(jih,'CRVAL1')
    ns = fxpar(jih,'CRVAL2')
    ; Find center, radius, p angle, b angle, etc.
    datec = FXPAR(jih,'DATE-OBS')
    timec = FXPAR(jih,'TIME-OBS')
    IF (STRMID(datec, 2, 1) EQ '/') THEN BEGIN
        READS,datec,day,mon,yrs,FORMAT='(i2,1x,i2,1x,i2)'
        yrs = yrs + 1900
    ENDIF ELSE BEGIN
        READS,datec,yrs,mon,day,FORMAT='(i4,1x,i2,1x,i2)'
    ENDELSE
    READS,timec,hour,min,sec,FORMAT='(i2,1x,i2,1x,i2)'
    hr = hour + ((min + (sec / 60.0)) / 60.0)
;   Use solar ephemeris from Johns Hopkins library
    SUN,yrs,mon,day,hr,pa=p,lat0=b,sd=rsun,long0=lnaught

    h = 256.0 + (imscl * ew)
    k = 240.0 - (imscl * ns)
    rad = imscl * rsun
    rad2 = rad * rad

    imout = MAKE_ARRAY(512,481,/byte,value=0)

    FOR i = 0, 511 DO BEGIN
        FOR j = 0, 480 DO BEGIN
           r1 = (i+1-h)^2 + (j+1-k)^2
           IF (r1 LE rad2) THEN BEGIN
              cost2 = 1.0 - (r1 / rad2)
              cost = SQRT(cost2)
              ldfunc = 1.00 - 0.88 + 0.23 + (0.88 * cost) - (0.23 * cost2)
              imout(i,j) = byte(255.0 * ldfunc)
           ENDIF
        ENDFOR
    ENDFOR
    ; imout is our fake Sun.

    ; Use fake sun to obtain fake flatfield image.
    jf = fltarr(512,481)
    nz = where(imout ne 0)
    z = where(imout eq 0, numzeros)
    IF (numzeros gt 0) THEN BEGIN
        jf(nz) = FLOAT(ji(nz))/imout(nz)
        jf(z) = 1.0
    ENDIF ELSE BEGIN
        jf = FLOAT(ji)/imout
    ENDELSE
    ; Replace bad columns
    ; For some reason, up to col 25 is bad sometimes!
    FOR i = 0, 25 DO jf(i,*) = jf(26,*)
    FOR i = 504, 511 DO jf(i,*) = jf(503,*)
    jmean = mean(jf)
    jflat = jf /jmean
    IF (numzeros gt 0) THEN BEGIN
        jflat(z) = 1.0
    ENDIF
    ; Smooth the fake flat
    jfs = SMOOTH(jflat,50,/EDGE_TRUNCATE)
    ; Renormalize
    jmean = mean(jfs)
    jflat = jfs / jmean

    ; Correct real solar image with fake flatfield.
    nz = where(jflat ne 0.)
    z = where(jflat eq 0., numzeros)
    jo = ji
    IF (numzeros gt 0) THEN BEGIN
        jo(nz) = ji(nz) / jflat(nz)
        jo(z) = 0.
    ENDIF ELSE BEGIN
        jo = ji / jflat
    ENDELSE

    ; Rebin images
    viss = CONGRID(vos,512,384)
    jiss = CONGRID(jo,512,384)

    finishup:

    return, jflat
end