% routine to grid data
%
% function [zi,xi,yi]=gridopts(x,y,z,gd,ngp,[zto],[xysc],[opt],[xiyi]);
%
% gd=1 griddata
% gd=2 sagas mincurvi
% gd=3 sagas kriging
% gd=4 sagas objmap
%
% zto= some number betwe 0 and 1 (for objmap & krig)
% xysc= some 2 element vector with x & y scales 
% opt= see help for specified method (typically [32 8 500 1/3 1])
% xiyi= ngpX2 matrix defining desired x & y grid nodes 

function [zi,xi,yi]=gridopts(x,y,z,gd,ngp,zto,xysc,opt,xiyi);

if nargin==3;
    gd=menu('Grid data?','use griddata',...
                         'use sagas mincurvi',...
                         'use sagas kriging',...
                         'use sagas objmap',...
                          'none');       
end;

%if gd>=2&gd<=4;

cp=computer;

if cp(1:3)=='SGI';
   path(path,'/home3/ocn/jmanning/atsee/saga');

else;

   path(path,'e:\atsee\saga');

end;   

%end;

  cgd=['griddata';'mincurvi';' kriging';'  objmap'];

  if gd<=4;
    disp(' ');
    disp(' ');
    disp('Hardcoded to grid data');

%   process differently if x & y are only vectors & z is matrix
    if length(x)~=length(y);
      xi=x;
      yi=y';

    else;
%     find distinct values of x & y
      disp('Finding distinct values...');

      [xy,nr,c]=unique([x y]);
      x=xy(:,1);
      y=xy(:,2);
      snr=sort(nr);
      id=find(diff(snr)==1);
      z=z(snr([id;id(length(id))+1]));
      if nargin<9;
        % calculated grid nodes
        disp(['Hardcoded to ',num2str(ngp),' grid pts in x & y both']);
        xi=[min(x):(max(x)-min(x))/ngp:max(x)];
        yi=[min(y):(max(y)-min(y))/ngp:max(y)]';
      else;
        % user specified grid nodes
        disp(['Hardcoded to ',num2str(length(xiyi)),' grid pts in x & y both']);
        xi=xiyi(:,1)';
        yi=xiyi(:,2);
      end;
    end;

%     define grid
    if nargin<5;
      ngp=25;
    end;
    [X,Y]=meshgrid(xi,yi);
      
    disp('Gridding in progress ...');

    if (gd==3|gd==4)&nargin<4;
      zto=input('Enter some number betw 0 and 1 (typically .5): ');
      xysc=input('Enter vector w/x&y scales (depends on your data units): ');
    end;
      
    if nargin<8;
      opt=[32 8 500 1/3 1];
    end;
    if gd==1;
      zi=griddata(x,y,z,X,Y);
    elseif gd==2;
      zi=mincurvi(x,y,z,X,Y);
    elseif gd==3;
      zi=kriging(x,y,z,X,Y,xysc,zto,opt);
    elseif gd==4;
      zi=objmap(x,y,z,X,Y,xysc,zto,opt);
    end;
  end;
