PRO TWODPOLY,DATA,X,Y,NORDER,A0,ACOF,DATAFIT,SIGMA
; Fits 2-d  polynomial to 1-D array DATA, with X & Y locations
; Works for  NORDER = 1-3; assumes equal weights
; Uses IDL User Library Routine REGRESS:  see it for description of output
npts=n_elements(data)
print,npts,' pixels used in fitting'
w=replicate(1.,npts)

case norder of
1:   begin
     m = fltarr(2,npts)  &  nc=2
     end
2:   begin
     m = fltarr(5,npts)  &  nc=5
     x2 = x^2    &    y2 = y^2
     m(2,*) = x2  &  m(3,*) = x*y  &  m(4,*) = y2
     end
3:   begin
     m = fltarr(9,npts)  &  nc=9
     x2 = x^2    &    y2 = y^2 
    m(2,*) = x2  &  m(3,*) = x*y  &  m(4,*) = y2
     m(5,*) = x*x2  &  m(6,*) = y*x2  &  m(7,*) = x*y2  &  m(8,*) = y*y2
     end
else:  begin
     print,'Works only for NORDER = 1-3'
     return
     end
endcase
m(0,*) = x  &  m(1,*) = y  &  d = data

acof=regress(m,d,w,dum,a0,sigma,Ftest,R,Rmul,Chisq,/relative_weight) 
print,'Constant term:',a0,'    F-test:',Ftest,'    Chi-squared:',Chisq
print,'   Coefficient     Sigma  Correlation Coefficient'
for i=0,nc-1 do print,acof(i),sigma(i),r(i)

datafit=replicate(a0,npts)
datafit = datafit + acof(0)*x + acof(1)*y
if (norder ge 2) then datafit = datafit + acof(2)*x2 + acof(3)*x*y + acof(4)*y2
if (norder eq 3) then datafit = datafit + x2*(acof(5)*x + acof(6)*y) $
 + y2*(acof(7)*x + acof(8)*y)
end


PRO POLYIMAGE,DATA,NORDER,A0,ACOF,XS=X,YS=Y,FIT=DF,SIGMA=SIG,MASK=MSK
; Fits 2-D  polynomial to 2-D array DATA, with 1-D arrays X & Y of coordinates
; Works for  NORDER = 1-3; assumes equal weights
; Uses IDL User Library Routine REGRESS:  see it for description of output

nx=n_elements(data(*,0))
ny=n_elements(data(0,*))
if (not keyword_set(x)) then x=(findgen(nx)-float(nx-1)/2.)/(float(nx)/2.)
if (not keyword_set(y)) then y=(findgen(ny)-float(ny-1)/2.)/(float(ny)/2.)
x2d=fltarr(nx,ny)  &  y2d=fltarr(nx,ny)
for i=0,ny-1 do x2d(*,i) = x
for i=0,nx-1 do y2d(i,*) = y
if (not keyword_set(msk)) then msk=replicate(1,nx,ny)
tvscl,x2d,0  &  tvscl,y2d,1  &  tvscl,msk,2
w = where(msk ne 0)
d=data(w)  &  xx=x2d(w)  &  yy=y2d(w)
TWODPOLY,D,XX,YY,NORDER,A0,ACOF,DF1,SIG
df=replicate(a0,nx,ny)
df(W) = df1
RETURN
END


FUNCTION POLYEVAL,NX,A0,ACOF,NYS=NY,XS=X,YS=Y,MASK=MSK
; Evaluates 2-D  polynomial to 2-D array DATA, with 1-D arrays X & Y of coordinates

if (not keyword_set(ny)) then ny=nx
if (not keyword_set(x)) then x=(findgen(nx)-float(nx-1)/2.)/(float(nx)/2.)
if (not keyword_set(y)) then y=(findgen(ny)-float(ny-1)/2.)/(float(ny)/2.)
x2d=fltarr(nx,ny)  &  y2d=fltarr(nx,ny)
for i=0,ny-1 do x2d(*,i) = x
for i=0,nx-1 do y2d(i,*) = y
if (not keyword_set(msk)) then msk=replicate(1,nx,ny)
; tvscl,x2d,0  &  tvscl,y2d,1  &  tvscl,msk,2
w = where(msk ne 0)
xx=x2d(w)  &  yy=y2d(w)  &  x2d=0  &  y2d=0
x2 = xx^2    &  y2 = yy^2

datafit = replicate(a0,nx,ny)
nc = n_elements(acof)
datafit(w) = datafit(w) + acof(0)*xx + acof(1)*yy
if (nc ge 2) then begin
    x2 = xx^2    &  y2 = yy^2
    datafit(w) = datafit(w) + acof(2)*x2 + acof(3)*xx*yy + acof(4)*y2
    end
if (nc gt 5) then datafit(w) = datafit(w) + x2*(acof(5)*xx + acof(6)*yy) $
 + y2*(acof(7)*xx + acof(8)*yy)
return,datafit
end