%function Z=bathy(x,y,params,[make_plot])
%
%  Usage example: Z=bathy32(0:.5:100,0:-.5:-100,[17.5,27.5,60,4,2],1);
%
%  model coordinates are right-handed with X north, Y east and Z down
%  params=[B1,B2,B3,Z1,Z2]
%    where: B1=distance from coast to 1st seafloor bend in X or Y direction
%           B2=distance from coast to 2nd seafloor bend in X or Y direction
%           B3=distance from coast to 3rd seafloor bend in X or Y direction
%           Z1=constant depth of seafloor between 1st and 2nd bends
%           Z2=constant depth of seafloor after 3rd bend
%           NOTE: Z1 must be >=0
%                 Z2 must also be >=0 but can be bigger or smaller
%                    than Z1. 
%                 If Z2<Z1 you have a trench; if Z2>Z1 you have a ridge
%  x and y are vectors that define grid for Z(x,y) 
%  x or y can be positive or negative
%  returns Z(x,y) >= 0 (ie right-handed as defined above)
%  ALL distances are in km
%
%  if optional argument make_plot > 0
%    plots a pcolor image with y east-west and x north-south
%    pauses
%    Then plots a depth y transect at user-selected x
%       
%  To create sites bathymetry you have to set x (North-south orientation), y
%  (East-West)
%   x   y
% -50	50  
% -20	20  
% -10	10  
% -6	6   
% -3	3	
% 3	-3  
% 6	-6  
% 10	-10 
% 15	-15 
% 20	-20 
% 25	-25 
% 35	-35 
% 50	-50 
% 75	-75 
% 100	-100
%
% The bathymetry output is written to bathydtm32_out.dat file
%


function Z=bathy(x,y,params,varargin)

  if nargin<3
     error('Needs at least 3 args: x,y and params')
  end
  if  nargin>3
    make_plot=varargin{:};
  else
    make_plot=0;
  end
  
%  params=[B1,B2,Z0,Z1]
  B1=params(1);
  B2=params(2);
  B3=params(3);
  Z1=params(4);
  Z2=params(5);

% compute slopes of seafloor
  slope1= Z1/B1;
  slope2= (Z2-Z1)/(B3-B2);

  yp=y;
  xp=x;
  
  [X,Y]=meshgrid(xp,yp);
  % fill all Z's with Z2
  Z=Z2*ones(size(X));
  if size(Z,1) == 1
    make_plot=0;  
  end    
  
  kkk=find( Y<(-1)*B2 & Y>=(-1)*B3 | X>B2 & X<=B3);
  if ~isempty(kkk)
      Z(kkk)=Z1+slope2*min(abs(Y(kkk))-B2,X(kkk)-B2);
  end

  kk=find(Y<(-1)*B1 & Y>=(-1)*B2 | X>B1 & X<=B2);
  if ~isempty(kk)
      Z(kk)=Z1*ones(size(kk));
  end

  k=find(Y>=(-1)*B1 | X<=B1);
  if ~isempty(k)
     Z(k)=slope1*min(abs(Y(k)),X(k));
  end

  k=find(Y>=0 | X<=0);
  if ~isempty(k)
     Z(k)=0;
  end
  
  Z=Z';

  i=1;
  for xx = 1 : length(x)
    for yy = 1 : length(y)
     Zout(i,1) = y(yy);
     Zout(i,2) = x(xx);
     Zout(i,3) = 1000*Z(xx,yy);
     i = i + 1;
    end
  end

  dlmwrite('bathydtm32_out.dat',Zout);
  
  if make_plot > 0
    close all
    figure
      pcolor(y,x,Z)
      shading interp
      axis equal
      ylabel('X km')
      xlabel('Y km')
      hcb=colorbar;
      set(hcb,'ydir','reverse')
      ncolors=round(4*Z1);
      colormap(flip(jet(ncolors)));
      title(['Seafloor bends at: ',...
                 num2str(B1),' ',num2str(B2),' ',num2str(B3),' km'])
      pause
    figure
      while 1
          x0=input('Enter position of x transect (0 to quit): ')
          if x0==0
            break
          end
          j=find(x==x0);
          if ~isempty(j)
            plot(y,Z(j,:),'linewidth',3)
            axis([-inf inf, 0,max(Z1,Z2)])
            set(gca,'ydir','reverse')
            title(['Seafloor transect at x = ',num2str(x0),' km'])
          else
            disp(['Your model does not have x0 = ',num2str(x0)])
          end  
      end
  end
end