FUNCTION MINFUNC,X,Y,XMIN
; Returns minimum of Y(X), at XMIN; uses spline interpolation
; TDT  7/17/93
ITMAX = 2
N = N_ELEMENTS(X)
YMIN = MIN(Y,IT)
IF (IT EQ 0) THEN BEGIN 
    YMIN = MIN(Y(1:N-1),IT) 
    IT = IT+1
ENDIF
IF (IT EQ (N-1)) THEN XS = X(N-2)+FINDGEN(FIX(X(N-1)-X(N-2))+1) 
IF ((IT GT 0) AND (IT LT (N-1))) THEN XS = X(IT-1) +FINDGEN(FIX(X(IT+1)-X(IT-1))+1) 
YS = SPLINE(X,Y,XS)
YMIN=MIN(YS,IT)
XMIN=XS(IT)
; plot,xs,ys & oplot,[xmin,xmin,xmin],[0,ymin,0] &  wait,2
RETURN,YMIN
END

FUNCTION MAXFUNC,X,Y,XMAX
T = -MINFUNC(X,-Y,XM)
XMAX = XM
RETURN, T
END

PRO UNWRAP,PH
; Unwraps phases PH(x,y,i) by adding or subtracting 360 OR 720 to eliminate
;  to eliminate big jumps as a function of i
; Then adjusts phase of middle to be near 0
; TDT  7/17/93
NP = N_ELEMENTS(PH(0,0,*))
TCOUNT=0.*INTARR(4)
FOR I=1,NP-1 DO BEGIN
    T = PH(*,*,I)
    DP0 = ABS(T-PH(*,*,I-1))
    DP = ABS(T+360-PH(*,*,I-1))
    W = WHERE(DP LT DP0,COUNT)
    TCOUNT(0)=TCOUNT(0)+COUNT
    IF (COUNT GT 0) THEN BEGIN
      T(W) = T(W) + 360.
      DP0 = ABS(T-PH(*,*,I-1))
      DP = ABS(T+360-PH(*,*,I-1))
      W = WHERE(DP LT DP0,KOUNT)
      IF (KOUNT GT 0) THEN T(W) = T(W) + 360.
      TCOUNT(1)=TCOUNT(1)+KOUNT
    ENDIF
    DP0 = ABS(T-PH(*,*,I-1))
    DP = ABS(T-360-PH(*,*,I-1))
    W = WHERE(DP LT DP0,COUNT)
    TCOUNT(2)=TCOUNT(2)+COUNT
    IF (COUNT GT 0) THEN BEGIN
       T(W) = T(W) - 360.
       DP0 = ABS(T-PH(*,*,I-1))
       DP = ABS(T-360-PH(*,*,I-1))
       W = WHERE(DP LT DP0,KOUNT)
       TCOUNT(3)=TCOUNT(3)+KOUNT
       IF (KOUNT GT 0) THEN  T(W) = T(W) - 360.
    ENDIF
    PH(*,*,I) = T
ENDFOR
PRINT,'Phase Unwrapping pixel counts: ',TCOUNT
NP2 = NP/2
T = (PH(*,*,NP2-1)+PH(*,*,NP2)+PH(*,*,NP2+1))/1080.
W = WHERE(T GE 0.,COUNT) & IF (COUNT GT 0) THEN T(W) = T(W)+0.5
W = WHERE(T LT 0.,COUNT) & IF (COUNT GT 0) THEN T(W) = T(W)-0.5
T = FIX(T)*360.
FOR I=0,NP-1 DO PH(*,*,I)=PH(*,*,I)-T
END

PRO DTPROF,L,AMP,PHAS,OFF,IMIN,IMAX,LMIN,PHMIN,LMAX,PHMAX
; Analyzes a set of trig fits to a detune sequence
; L = array of wavelength settings of M2
; AMP, PHAS, OFF = results of SINFIT for M1 detuning at each of
;   the positions of M2
; Returns IMIN (IMAX) = minimum (maximum) intensities for the
;   entire 2-D set of detunings
; LMIN (LMAX) = M2 wavelength at minimum (maximum) intensity
; PHMIN (PHMAX) = phase in degrees of M1 with respect to M2 at the
;   minimum (maximum) intensity
; TDT  7/17/93
NX = N_ELEMENTS(AMP(*,0,0))  &  NY = N_ELEMENTS(AMP(0,*,0))
NL = N_ELEMENTS(AMP(0,0,*))
IMIN = FLTARR(NX,NY)
IMAX = FLTARR(NX,NY)
PHMIN = FLTARR(NX,NY)
PHMAX = FLTARR(NX,NY)
LMIN = FLTARR(NX,NY)
LMAX = FLTARR(NX,NY)
FOR J=0,NY-1 DO BEGIN
    FOR I=0,NX-1 DO BEGIN
    IMIN(I,J) = MINFUNC(L,OFF(I,J,*)-AMP(I,J,*),XMIN)
    LMIN(I,J) = XMIN
    PHMIN(I,J) = SPLINE(L,PHAS(I,J,*),XMIN)
    IMAX(I,J) = MAXFUNC(L,OFF(I,J,*)+AMP(I,J,*),XMAX)
    LMAX(I,J) = XMAX
    PHMAX(I,J) = SPLINE(L,PHAS(I,J,*),XMAX)+180.
    ENDFOR
if (((j/8)*8 eq j) and (j gt 0)) then begin & plot,lmin(*,j) & end	; pause & end
ENDFOR
RETURN
END

PRO DTMAP,L,AMP,PHAS,OFF,M1,M2,LCTR,CONT,M2CONT,noreorder=noreorder,reorder=reorder
;
;+
;NAME:
;	DTMAP
;PURPOSE:
;	Derives Michelson maps and corrected images from a detune sequence
;INPUT:
;	L 	- array of wavelength settings of M2
;	    findgen(10)/9*400 - 200
;	AMP, PHAS, OFF = results of SINFIT for M1 detuning at each of
;		the positions of M2
;OUTPUT:
;	M1,M2	- Michelson tunings at solar line center, in mA
; 	LCTR (CONT) - intensity at line center (continuum) for
;		properly tuned Michelsons
;	M2CONT	- Michelson 2 tuning at continuum
;OPTIONAL KEYWORD INPUT:
;	noreorder- If set, then do not shift the data around 
;		  Normally the data must be shifted so that the
;		  minimum intensities are centered in the array
;		  so that splines will work
;	reorder	- The order of the images that should be sorted
;		  The default is: [1,2,3,4,5,6,7,8, 9, 1,2]
;		  (9 and 0 are the same)
;HISTORY:
;	Written by Ted Tarbell 17-Jul-1993
;	16-Feb-94 (MDM) - Modified to generalize
;	26-May-94 (MDM) - Added REORDER keyword option
;-

NX = N_ELEMENTS(AMP(*,0,0))
NY = N_ELEMENTS(AMP(0,*,0))
NL = N_ELEMENTS(AMP(0,0,*))	;is generally 10
FSR1 = 200. 
FSR2 = 400.

; The following does a cyclic shift and extension of the data as a
;    function of M2 tuning--sometimes the amount and/or direction
;    of shift have to be changed, depending on where the solar line
;    center falls in the tuning scan of M2
;
if (keyword_set(noreorder)) then begin
    A = AMP
    P = PHAS
    O = OFF
    LL = L
end else begin
    if (n_elements(reorder) eq 1) then begin
	ss = [1,2,3,4,5,6,7,8, 9, 1,2]	;9 and 0 are the same
    end

    LL = [L(1:NL-2),L(NL-1),L(1:2)+FSR2]
    A = FLTARR(NX,NY,NL+1)
    P = FLTARR(NX,NY,NL+1)
    O = FLTARR(NX,NY,NL+1)

    A(*,*,0:NL-3) = AMP(*,*,1:NL-2)			;1,2,3,4,5,6,7,8
    A(*,*,NL-2) = 0.5*(AMP(*,*,0)+AMP(*,*,NL-1))	;0+9
    A(*,*,NL-1:NL) = AMP(*,*,1:2)			;1,2

    P(*,*,0:NL-3) = PHAS(*,*,1:NL-2)
    P(*,*,NL-2) = 0.5*(PHAS(*,*,0)+PHAS(*,*,NL-1))
    P(*,*,NL-1:NL) = PHAS(*,*,1:2)

    O(*,*,0:NL-3) = OFF(*,*,1:NL-2)
    O(*,*,NL-2) = 0.5*(OFF(*,*,0)+OFF(*,*,NL-1))
    O(*,*,NL-1:NL) = OFF(*,*,1:2)
end
;
print,'Unwrapping raw phases versus wavelength'
UNWRAP,P
nx8 = nx/8  &  nx2 =nx/2  &  ny2 = ny/2
plot,ll,p(nx2,ny2,*),title='Phase',yrange=[min(p),max(p)]
for i=nx8,nx-nx8,8 do oplot,ll,p(i,ny2,*)
pause
plot,ll,a(nx2,ny2,*),title='Amplitude',yrange=[min(a),max(a)]
for i=nx8,nx-nx8,8 do oplot,ll,a(i,ny2,*)
pause
plot,ll,o(nx/2,ny/2,*)-a(nx/2,ny/2,*),title='MINIMUM INTENSITY',yrange=[min(o-a),max(o-a)]
for i=nx8,nx-nx8,8 do oplot,ll,o(i,ny/2,*)-a(i,ny/2,*)
pause
DTPROF,LL,A,P,O,LCTR,CONT,M2,PHMIN,M2CONT,PHMAX
P = 0  &  A = 0  &  O = 0
PH1 = FLTARR(NX,2,NY) & PH1(*,0,*)=PHMIN & PH1(*,1,*)=PHMIN 
print,'Unwrapping M1 spatially'
UNWRAP,PH1
PH2 = FLTARR(NX,2,NY) & PH2(*,0,*)=M2*(360./FSR2) & PH2(*,1,*)=PH2(*,0,*) 
print,'Unwrapping M2 spatially'
UNWRAP,PH2
M2(*,*) = PH2(*,0,*)*(FSR2/360.)
;;M1(*,*) = PH1(*,0,*)*(FSR1/360.)
M1 = reform(PH1(*,0,*)*(FSR1/360.))
RETURN
END