function [ZREG,Zreg]=polyregional(XBIG,YBIG,ZBIG,Xsmall,Ysmall,Zsmall)
%
% [ZREG,Zreg ]=polyregional(XBIG,YBIG,ZBIG,Xsmall,Ysmall,Zsmall);
%
% polyregional least squares fits a polynomial surface to a
% regional scale grid XBIG YBIG ZBIG and projects this regional (ZREG)
% onto a smaller grid Xsmall Ysmall Zsmall contained within the regional scale grid
% the regional projected onto the small grid is Zreg
% anomalous regions within the regional scale grid can be masked
% out so that they do not influence the polynomial fit
%
% [ZREG,Zreg ]=polyregional(XBIG,YBIG,ZBIG,Xsmall,Ysmall,Zsmall);

%***********************************************************************************
%
%                              Jim Merriam
%                       University of Saskatchewan
%                          jim.merriam@usask.ca
%                              (Dec, 1997)  
%                            (last Oct, 1999)
%***********************************************************************************
%
%
% modify this so that the regional also makes the perimeter as small as possible -
% or is this a task for LAPLACE? 
% ALso work on a Fourier function fitting version

clf
H1=axes('position',[.05 .05 .4 .55]); % axes for initial contour maps
contour(XBIG,YBIG,ZBIG,20),title('REGIONAL SCALE GRID'),colorbar('horiz')
hold on
[my,ny]=size(Ysmall);
[my,ny]=size(Xsmall);
  boxtop=Ysmall(my,1); % top of box around mask
  boxbot=Ysmall(1,1); % bottom of box around mask
  boxleft=Xsmall(1,1);% left side of mask
  boxright=Xsmall(1,ny);% right side of mask
  hb=plot([boxleft boxright boxright boxleft boxleft],[boxbot boxbot boxtop ...
  boxtop boxbot],'w-');
  set(hb,'LineWidth',1);
H2=axes('position',[.5 .39 .4 .55]);
contour(Xsmall,Ysmall,Zsmall,20),title('SMALL SCALE GRID'),colorbar('horiz')
disp('ENTER TO CONTINUE');pause;

ans=input(' MASK A SUBGRID? [y]/n ','s');
 if(isempty(ans) | ans=='Y' | ans=='y');
 ZMASK=gridmask(XBIG,YBIG,ZBIG); % call gridmask to mask out regional
 else
 ZMASK=ZBIG;
 end

pause(2)
  KeepAtIt='y';
  while KeepAtIt=='y'
    badchoice='y';
    while badchoice=='y'
    disp('         SELECT LINEAR, QUADRATIC, OR CUBIC 2D POLYNOMIAL')
    disp(' 1 FOR LINEAR')
    disp(' 2 FOR QUADRATIC')
    disp(' 3 FOR CUBIC')
    disp(' 4 FOR QUARTIC')
    response=input('SELECTION? ');
      if response==1;level='linear';badchoice='n';
      elseif response==2;level='quad  ';badchoice='n';
      elseif response==3;level='cubic ';badchoice='n';
      elseif response==4;level='quart ';badchoice='n';    
      else disp('THIS IS NOT A VALID CHOICE. TRY AGAIN')
      end
    end
  disp(level);    
  XX=XBIG;
  YY=YBIG;
  [n,m]=size(ZMASK);
  ZZ=reshape(ZMASK,n*m,1); % reshaped ZMASK etc
  XX=reshape(XBIG,n*m,1);
  YY=reshape(YBIG,n*m,1);
  i=isnan(ZZ); % find Nan's in masked grid
  ZZ(i)=[]; %clear grid points with NaN's
  XX(i)=[];
  YY(i)=[];
%       calculate the model matrix

    if(level=='linear')
    A=[ones(size(ZZ)) XX YY];
    elseif(level=='quad  ')
    A=[ones(size(ZZ)) XX YY XX.^2 YY.^2 XX.*YY];
    elseif(level=='cubic ')
    A=[ones(size(ZZ)) XX YY XX.^2 YY.^2 XX.*YY XX.*YY.^2 YY.*XX.^2 XX.^3 YY.^3];    
    elseif(level=='quart ')
    A=[ones(size(ZZ)) XX YY XX.^2 YY.^2 XX.*YY XX.*YY.^2 YY.*XX.^2 XX.^3 YY.^3 XX.*YY.*3 XX.^3.*YY (XX.^2).*(YY.^2) XX.^4 YY.^4 ];    
    end
    
  polynomialsurface=A\ZZ;
  disp(' ') 
    if(level=='linear')
    disp('    COEFFICIENTS OF LINEAR REGIONAL '); 
    disp(['  MEAN   X    Y']);
    S=sprintf(' %4.1f %4.1f %4.1f ',polynomialsurface);
    disp(S)
    elseif(level=='quad  ')
    disp('    COEFFICIENTS OF QUADRATIC REGIONAL');   
    disp(['  MEAN   X    Y     X^2     Y^2       XY']);
    S=sprintf(' %4.1f %4.1f %4.1f %4.1e %4.1e %4.1e',polynomialsurface);
    disp(S)
    elseif(level=='cubic ')
    disp('    COEFFICIENTS OF CUBIC REGIONAL ');
    disp(['  MEAN   X    Y     X^2     Y^2       XY      XY^2       YX^2      X^3      Y^3']);
    S=sprintf(' %4.1f %4.1f %4.1f %4.1e %4.1e %4.1e %4.1e %4.1e %4.1e %4.1e %4.1e ',polynomialsurface);
    disp(S)
    elseif(level=='quart ')
    disp('    COEFFICIENTS OF QUARTIC REGIONAL ')
        disp(['  MEAN   X    Y     X^2     Y^2       XY      XY^2       YX^2      X^3      Y^3    XY^3 YX^3 X^2Y^2  Z^4      Y^4' ]);
    S=sprintf(' %4.1f %4.1f %4.1f %4.1e %4.1e %4.1e %4.1e %4.1e %4.1e %4.1e %4.1e %4.1e %4.1e %4.1e %4.1e % 4.1e',polynomialsurface);
    disp(S)

    end
    disp(' ');
    
  [n,m]=size(XBIG);
  ZZ=reshape(ZBIG,n*m,1); % reshape ZREG to etc
  XX=reshape(XBIG,n*m,1); % tnew XX and YY
  YY=reshape(YBIG,n*m,1); % on unmasked regional scale grid

%        calculate the model matrix

    if(level=='linear')
    B=[ones(size(ZZ)) XX YY];
    elseif(level=='quad  ')
    B=[ones(size(ZZ)) XX YY XX.^2 YY.^2 XX.*YY];
    elseif(level=='cubic ')
    B=[ones(size(ZZ)) XX YY XX.^2 YY.^2 XX.*YY XX.*YY.^2 YY.*XX.^2 XX.^3 YY.^3];   
    elseif(level=='quart ')
    B=[ones(size(ZZ)) XX YY XX.^2 YY.^2 XX.*YY XX.*YY.^2 YY.*XX.^2 XX.^3 YY.^3 XX.*YY.*3 XX.^3.*YY (XX.^2).*(YY.^2) XX.^4 YY.^4];    
    end

  ZREG=B*polynomialsurface; 
  ZREG=reshape(ZREG,n,m); %polynomial surface on regional scale
  clf

  
  H1=axes('position',[.05 .05 .4 .55]); % axes for regional scale grid
  contour(XBIG,YBIG,ZMASK,20),title('(MASKED) REGIONAL SCALE GRID'),colorbar('horiz')
  H2=axes('position',[.5 .39 .4 .55]); % axes for regional on regional
%                                      scale grid
  contour(XBIG,YBIG,ZREG,20);title('REGIONAL ON REGIONAL SCALE GRID'),colorbar('horiz')
  disp('ENTER TO CONTINUE'),pause,clf

  H1=axes('position',[.05 .05 .4 .55]); % axes for regional scale grid
  contour(XBIG,YBIG,ZMASK,20),title('(MASKED) REGIONAL SCALE GRID'),colorbar('horiz')
  H2=axes('position',[.5 .39 .4 .55]); % axes for residual on regional
%                                      scale grid
  contour(XBIG,YBIG,ZBIG-ZREG,20);title('RESIDUAL ON REGIONAL SCALE GRID'),colorbar('horiz')
  disp('ENTER TO CONTINUE'),pause,clf


  [n,m]=size(Xsmall);
  XX=reshape(Xsmall,n*m,1); % new XX etc for small scale grid
  YY=reshape(Ysmall,n*m,1);
  ZZ=reshape(Zsmall,n*m,1);
  clf 

%         calculate the model matrix

    if(level=='linear')
    C=[ones(size(ZZ)) XX YY];
    elseif(level=='quad  ')
    C=[ones(size(ZZ)) XX YY XX.^2 YY.^2 XX.*YY];
    elseif(level=='cubic ')
    C=[ones(size(ZZ)) XX YY XX.^2 YY.^2 XX.*YY XX.*YY.^2 YY.*XX.^2 XX.^3 YY.^3]; 
    elseif(level=='quart ')
    C=[ones(size(ZZ)) XX YY XX.^2 YY.^2 XX.*YY XX.*YY.^2 YY.*XX.^2 XX.^3 YY.^3 XX.*YY.*3 XX.^3.*YY (XX.^2).*(YY.^2) XX.^4 YY.^4];    

    end

  Zreg=C*polynomialsurface; % regional projected on samll scale grid
  Zreg=reshape(Zreg,n,m);
  H1=axes('position',[.05 .05 .4 .55]); % axes for original 
%                                        small scale grid
  contour(Xsmall,Ysmall,Zsmall,20),title('ORIGINAL DATA GRID');colorbar('horiz')
  H2=axes('position',[.5 .39 .4 .55]); % axes for regional projected  
%                                      onto small scale grid
  contour(Xsmall,Ysmall,Zreg,20),title('REGIONAL PROJECTED ONTO DATA GRID'),colorbar('horiz')
  disp('ENTER TO CONTINUE'),pause,clf

  H1=axes('position',[.05 .05 .4 .55]);% axes for original on 
%                                       small scale grid
  contour(Xsmall,Ysmall,Zsmall,20),title('ORIGINAL DATA GRID');colorbar('horiz')
  H2=axes('position',[.5 .39 .4 .55]); % axes for residual on small
%                                     scale grid
  contour(Xsmall,Ysmall,Zsmall-Zreg,20),title('RESIDUAL DATA GRID'),colorbar('horiz')

    KeepAtIt=input('TRY A DIFFERENT POLYNOMIAL? [y]/n ','s');
    if isempty(KeepAtIt)
    KeepAtIt='y'
    end

  end