PRO MAPLASER, IMAGE, M1, M2, trans, OBSMODE, pos, s1, s2, amp1, amp2, $
	phas1, phas2, offset1, offset2, m1sigma, m2sigma, uncert, sum_arr, $
	tsigma, $
	qfixed_dark=qfixed_dark, progver=progver, qstop=qstop, $
	m1img_avoid=m1img_avoid, m2img_avoid=m2img_avoid, nodefringe=nodefringe, $
	all_heads=all_heads
;
;NAME:
;	MAPLASER
;PURPOSE:
;	Produce the central wavelength and peak transmission maps
;	for one set of 27 or 29 images.
;INPUT:
;	image	- a NxNxM data array where M is 27 or 29 depending
;		  on the sequence
;OUTPUT:
;	m1	- The M1 central wavelength array
;	m2	- The M2 central wavelength array
;	trans	- The transmission
;	obsmode	- The image showing the location of the laser on the CCD
;	pos	- The location of the laser on the CCD [x0,y0]
;OPTIONAL KEYWORD INPUT:
;	nodefringe - If set, then do not perform defringing
;COMMENTS:
;	29 set: 1 obs, 1 dark, 10 m1 varies, 10 m2 varies, 
;	        1 obs, 1 dark, 4 "nominal", 1 "continuum"
;	27 set: 1 obs, 4 "nominal", 1 dust (?), 1 dark
;		10 m1 varies, 10 m2 varies
;	21 set: 1 dark, 10 m1 varies, 10 m2 varies
;	22 set: 1 dark, 1 obs, 10 m1 varies, 10 m2 varies
;HISTORY:
;	Written by Ted Tarbell (MTMTUNE and MAPLASER687)
;	19-Jan-94 (MDM) - Renamed to MAPLASER and generalized 
;V2.0	 3-Feb-94 (MDM) - Added /PROGVER and started tracking version
;			  numbers
;			- Added uncertainty and SUM_ARR outputs
;V2.01   7-Feb-94 (MDM) - Expanded to handle 24 image series 575-586
;V2.02	 8-Feb-94 (MDM) - Modified SINFIT to avoid the /0 case and to
;			  not create infinity phase values
;V2.03	10-Feb-94 (MDM) - Corrected another case of wrap around on
;			  the M1 and M2 result
;			- Added /QSTOP option
;V2.04	15-Feb-94 (MDM) - Added OFFSET1 and OFFSET2 as output params
;			- Modified to allow 21 image set (series 504
;			  solar dtune)
;V2.05	 4-Mar-94 (MDM) - Modified to handle 22 image set
;			- Expanded to return the sigma on the M1 and
;			  M2 fits
;V2.06	 7-Mar-94 (MDM) - Corrected SINFIT error in SIGMA calculation
;			- Added IMG_AVOID keyword input option
;V2.07	 8-Mar-94 (MDM) - Added average PHASE1 and PHASE2 to SUM_ARR
;V2.08	25-Mar-94 (MDM) - Added /nodefringe option
;V2.10	12-Apr-94 (MDM) - ** BIG CHANGE TO TRANSMISSION CALCULATION**
;			  Previously, the program assumed that the
;			  starting position of the M1 scan was the
;			  value that M1 sat at during the M2 scan
;			  (and like-wise for M2).  The new code will
;			  shift so that the transmission can be 
;			  calculated even if they were not set up
;			  properly.
;V2.11	13-Apr-94 (MDM) - Derive the location of the OBS image using
;			  SAT_OBS_LOC instead of CENTROIDW
;V2.2	27-Apr-94 (MDM) - Corrected error with SAT_OBS_LOC which was
;			  causing "Segmentation fault (core dumped)"
;			  AND "bus error (core dumped)"
;V2.3	15-Jun-94 (MDM) - Added call to MMAP_IMGNUM to figure out
;			  the images to use in the set
;			- Corrected error in sigma calculation for M1 & M2
;			  (see SINUS)
;			- Modified SINFIT to get alphabeta
;			- Added "Tsigma" as output image
;			- Adjusted how wrap around is fixed (170 vs 180)
;V2.4	17-Jun-94 (MDM) - Moved a bunch of "wrap" code into MMAP_UNWRAP
;V3.0	 1-Sep-94 (MDM) - Corrected error in dark image calculation.  When
;			  more than one dark image was specified with "idark"
;			  it ended up subtracting zero as the dark frame.
;-
;
if (keyword_set(progver)) then begin
    image = '3.000'
    return
end
;
; Needs a routine to locate the bright spot in OBSMODE, in CCD 
;   coordinates where x & y range from -1 to +1  
; TDT  1/18/94 
 
FSR1 = 200.  &  FSR2 = 400. 
DIM = SIZE(IMAGE) 
NX = DIM(1) 
NY = DIM(2) 
N  = DIM(3) 
NP = FLOAT(NX)*NY

mmap_imgnum, n, idark, iiobs, im1, im2, comment
if (im1 lt 0) then begin
    print, 'Returning...
    return
end
;
dark = image(*,*,idark)
;;if (n_elements(idark) gt 1) then dark = fix(total(dark,3) / n_elements(dark))
if (n_elements(idark) gt 1) then dark = fix(total(dark,3) / n_elements(idark))	;corrected error 1-Sep-94
if (keyword_set(qfixed_dark)) then begin
    dark = rfits('/home/mons/soft/mdi/data/dark_image.fits')
    print, 'Using fixed dark image "DARK_IMAGE.FITS"'
end
; 
FOR I=im1,im1+9 DO IMAGE(*,*,I)=IMAGE(*,*,I)-DARK 
FOR I=im2,im2+9 DO IMAGE(*,*,I)=IMAGE(*,*,I)-DARK 

MM1 = FLTARR(10)
MM2 = MM1 
FOR I=im1,im1+9 DO MM1(I-im1) = TOTAL(IMAGE(*,*,I))/NP
FOR I=im2,im2+9 DO MM2(I-im2) = TOTAL(IMAGE(*,*,I))/NP

ANGLE = FINDGEN(10)/9.*360.-180. 
if (getenv('ys_batch') eq "") then begin
    PLOT,ANGLE,MM1,LINE=0  & WAIT,.1
    oPLOT,ANGLE,MM2,LINE=1 &  WAIT,.1
end

; Fringing spatial frequencies set by trial & error
;FRINGE = [NX/2.,NY/2.]  &  FWHM = 1.
FRINGE = [NX-16.,4.]   &  FWHM = 12.  
if (keyword_set(nodefringe)) then fwhm = -1

SINFIT,ANGLE,IMAGE(*,*,im1:im1+9), PHAS1, AMP1, OFFSET1, m1sigma, RAD, FRINGE, alphabeta1, $
				FWHM=fwhm, img_avoid=m1img_avoid
SINFIT,ANGLE,IMAGE(*,*,im2:im2+9), PHAS2, AMP2, OFFSET2, m2sigma, RAD, FRINGE, alphabeta2, $
				FWHM=fwhm, img_avoid=m2img_avoid
;
PRINT,'Mean Phase, Amplitude, Offset for M1, M2:' 
PRINT,TOTAL(PHAS1)/NP,TOTAL(AMP1)/NP,TOTAL(OFFSET1)/NP 
PRINT,TOTAL(PHAS2)/NP,TOTAL(AMP2)/NP,TOTAL(OFFSET2)/NP 
;
if (keyword_set(qstop)) then stop
phas1 = mmap_unwrap(phas1)
phas2 = mmap_unwrap(phas2)
;
M1 = (PHAS1+180.)*FSR1/360. 
M2 = (PHAS2+180.)*FSR2/360.
;
phase1_shift = -180.
phase2_shift = -180.
if (sxpar(all_heads(*,im1), 'M1TUNE') ne 0) then begin
    m1_m1scan = sxpar(all_heads(*,im1), 'M1TUNE') ;almost always 26
    m2_m1scan = sxpar(all_heads(*,im1), 'M2TUNE') ;variable
    m1_m2scan = sxpar(all_heads(*,im2), 'M1TUNE') ;variable
    m2_m2scan = sxpar(all_heads(*,im2), 'M2TUNE') ;almost always 78
    phase1_shift = -180 + 360/45.*(m1_m2scan - m1_m1scan)
    phase2_shift = -180 + 360/45.*(m2_m1scan - m2_m2scan)
    print, 'PHASE1 and 2 SHIFT: ', phase1_shift, phase2_shift
end
S1 = OFFSET1 - AMP1*COS((phase1_shift - PHAS1)/!RADEG)
S2 = OFFSET2 - AMP2*COS((phase2_shift - PHAS2)/!RADEG)
I0av = (s1+s1)/2.
TRANS = (AMP1+OFFSET1)*(AMP2+OFFSET2) / I0av
PRINT,'Origin from M1, M2 scans, Mean Peak Transmission'
PRINT,TOTAL(S1)/NP,TOTAL(S2)/NP,TOTAL(TRANS)/NP
;
n = float(n_elements(angle))	;10 always
sigT1 = m1sigma^2 * ( ( 3. / (offset1+amp1)^2 ) + $
			( 9. / (4*I0av^2) ) - $
			( ( amp1 + 4*alphabeta1 ) / (amp1*(offset1+amp1)*I0av) ) )
sigT2 = m2sigma^2 * ( ( 3. / (offset2+amp2)^2 ) + $
			( 9. / (4*I0av^2) ) - $
			( ( amp2 + 4*alphabeta2 ) / (amp2*(offset2+amp2)*I0av) ) )
Tsigma = sqrt(1/n * (sigT1 + sigT2) )

;
if (iiobs ge 0) then begin
    OBSMODE = (IMAGE(*,*,iiobs) - dark)>0
    ;;centroidw, obsmode>0, x0, y0
    sat_obs_loc, obsmode, x0, y0
    pos = [x0,y0]
end else begin
    pos = [0,0]
end
;
;-- Added uncertainty 3-Feb-94

o1 = total(offset1)/n_elements(offset1)
o2 = total(offset2)/n_elements(offset2)
a1 = total(amp1)/n_elements(amp1)
a2 = total(amp2)/n_elements(amp2)
savg = total(s1+s2)/n_elements(s1)/2
davg = total(dark)/n_elements(dark)
p1 = total(phas1)/n_elements(phas1)
p2 = total(phas2)/n_elements(phas2)
;
n = float(n_elements(angle))
uncert = ( sqrt(3/(2*n)) * sqrt(1 + (.01*o1)^2) ) / (o1+a1) + $
	( sqrt(3/(2*n)) * sqrt(1 + (.01*o2)^2) ) / (o2+a2) - $
	( sqrt(3/(2*n)) * sqrt(1 + (.01*savg)^2) ) / savg
sum_arr = [o1,o2,a1,a2,savg,davg,p1,p2]
RETURN 
END