;+
; NAME: 
;	WAVECAL.pro
;
; PURPOSE: 
;	To read a list of (lambda, x, y) and determine the polynomial coefs of 
;	the form:
;	
;	lambda(x,y) = a + bx + cy + dxy + ex^2 + fy^2        (if order = 2)
;   lambda(x,y) = a + bx + cy + dxy + ex^2 + fy^2 + 
;				  gx^2y + hxy^2 + ix^3 + jy^3.           (if order = 3)
;
;	The function supports up to order 4.
;
; DESCRIPTION:
; 	This procedure reads in a list of (x,y) points and their associated 
;	wavelengths in microns.  The (lambda, x, y) triplet list is generated 
;	by the procedure GetClix.  This procedure takes the (x,y) data and 
;	constructs the matrix A where Ax=b.  In our case b is the vector of 
;	wavelengths and x are the parameters that we are solving for.   
;
;	The parameter array is [a,b,c,d,e,f, etc].  Its length depends on the order
;	of the polynomial being fit.
;
; 	After constructing the matrix A, we call SVDC to do a singular value 
;	decomposition of the matrix into U, V, W.  Then the system of equations 
;	is solved with SVDSOL.  The parameter array is returned.  
; 
; 	Next the residual array is calculated along with the chi-squared.
;
; CATEGORY:
;
; CALLING SEQUENCE:
;  WAVECAL pts, order, dimX, dimY, par, resid, chiSq, lammap
;
; INPUTS:
;	pts   - (x,y, lambda) points corresponding the the wavlelengths in lambdas
;	order - the order of the polynomial to be fit
;   dimX  - the dimension in the x-direction of the output array, lammap
;   dimY  - the dimension in the y-direction of the output array, lammap
;
; OUTPUTS:
;   par  - The fitted parameters
;   resid - The fit residuals
;   chiSq - The chi-squared value of the fit
;   lammap -  a 2-D array of the wavelength map
;
; KEYWORD OUTPUT PARAMETERS:
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
; 19-APR-2004, CBO, SwRI
;-
pro	WAVECAL, pts, order, dimX, dimY, par, resid, chiSq, lammap


lambdas = REFORM(pts[*,0])

xys = pts[*,1:2]

npts = N_ELEMENTS(lambdas)

x = xys[*,0]
y = xys[*,1]

xv = REPLICATE(1.0, dimY)##FINDGEN(dimX)
yv=FINDGEN(dimY)##REPLICATE(1.0,dimX)

if (order EQ 2) then begin
	a=FLTARR( 6, npts)

	a[0,*]=1.0
	a[1,*]=x
	a[2,*]=y
	a[3,*]=x*y
	a[4,*]=x^2
	a[5,*]=y^2


	SVDC, a, w, u, v, /double
	par = SVSOL(u, w, v, lambdas, /double)

	calculated = par[0] + par[1]*x + par[2]*y + $
					par[3]*x*y + par[4]*x^2 + par[5]*y^2

	lammap = par[0]+par[1]*xv+par[2]*yv+par[3]*xv*yv+par[4]*xv^2+par[5]*yv^2

endif

if (order EQ 3) then begin
;  lambda(x,y) = a + bx + cy + dxy + ex^2 + fy^2 + gx^2y + hxy^2 + ix^3 + jy^3.  
;  The parameter array is [a,b,c,d,e,f, g, h, i, j].
	a=FLTARR( 10, npts)

	a[0,*]=1.0
	a[1,*]=x
	a[2,*]=y
	a[3,*]=x*y
	a[4,*]=x^2
	a[5,*]=y^2
	a[6,*]=y*x^2
	a[7,*]=x*y^2
	a[8,*]=x^3
	a[9,*]=y^3


	SVDC, a, w, u, v, /double
	par = SVSOL(u, w, v, lambdas, /double)

	calculated = par[0] + par[1]*x + par[2]*y + $
					par[3]*x*y + par[4]*x^2 + par[5]*y^2 + $
					par[6]*(x^2)*y + par[7]*(y^2)*x + par[8]*x^3 + par[9]*y^3

	lammap = par[0]+par[1]*xv+par[2]*yv+par[3]*xv*yv+par[4]*xv^2+par[5]*yv^2 + $
		 par[6]*(xv^2)*yv + par[7]*(yv^2)*xv + par[8]*xv^3 + par[9]*yv^3

endif


if (order EQ 4) then begin
	a=FLTARR( 15, npts)

	a[0,*]=1.0
	a[1,*]=x
	a[2,*]=y
	a[3,*]=x*y
	a[4,*]=x^2
	a[5,*]=y^2
	a[6,*]=y*x^2
	a[7,*]=x*y^2
	a[8,*]=x^3
	a[9,*]=y^3
	a[10,*]=(x^2)*(y^2)
	a[11,*]=(x^3)*y
	a[12,*]=x*(y^3)
	a[13,*]=x^4
	a[14,*]=y^4

	SVDC, a, w, u, v, /double
	par = SVSOL(u, w, v, lambdas, /double)

	calculated = par[0] + par[1]*x + par[2]*y + $
					par[3]*x*y + par[4]*x^2 + par[5]*y^2 + $
					par[6]*(x^2)*y + par[7]*(y^2)*x + par[8]*x^3 + par[9]*y^3 + $
					par[10]*(x^2)*(y^2) + par[11]*(x^3)*(y) + par[12]*(x)*(y^3) + $
					par[13]*x^4 + par[14]*y^4
					

	lammap = par[0]+par[1]*xv+par[2]*yv+par[3]*xv*yv+par[4]*xv^2+par[5]*yv^2 + $
		 par[6]*(xv^2)*yv + par[7]*(yv^2)*xv + par[8]*xv^3 + par[9]*yv^3 + $
		 par[10]*(xv^2)*(yv^2) + par[11]*(xv^3)*(yv) + par[12]*(xv)*(yv^3) + $
		 par[13]*xv^4 + par[14]*yv^4

endif

resid = lambdas - calculated

chiSq=TOTAL(resid^2)



end