function [bcs] = ampphase(bndele,meshnod,bathy);
% [bcs] = ampphase(bndele,meshnod,bathy);
% mesh.bcs.s2c 
%  Works on bndele elements which are flagged as open water  bndele(:,5)=5
% amplitude and phase of sea surface for M2 forcing
% 
%  Two crossing Kelvin waves working now.

disp(['meshtools/ampphase.m DEFAULT ampphase from meshtools directory WARNING'])
disp(['ampphase.m crashes if bndele and meshele are not good'])

% Kelvin wave solution for exp decay of height off shore
      TIDEP = 12.42*3600;
      FFREQ = 2.*pi/TIDEP;
      Ho = 20;
      OMEGA2 = 1.4584E-4 ;
      RAD = 6.3675E6 ;
      DEG = 31.9333;
      f = OMEGA2*sin(DEG/180.*pi); disp(['Coriolis = ',num2str(f),'  s-1'])
      Co = sqrt(9.8*Ho);
      Ro = Co/f  ;  disp(['Rossby Radius = ',num2str(Ro),'  m'])
%  ffreq = Co k  
      k = FFREQ/Co;

theta = 135 * pi/180  ;    % Coastline angle  135 is Wassaw island
disp(['Coastline angle is ',num2str(theta*180/pi)])
XO = 0;    %  Origin on coastline of wave
YO = 0;

theta = 90  * pi/180  ;    % Middle of nowhere
XO = 0;   
%  Origin on coastline of wave
YO = 0;  
% Amplitude of wave that decays off shore  (flow parallel coast):
  Ampoff   = 1.20;
% Amplitude of wave that decays along shore (flow perpendicular coast):
  Ampalong = .70;
%  Aoff=.7  along = .7  produces 1:2 ratio ellipses of -1 slope
%  Aoff=.7, along=1.2  produces lines of -1 slope
%  Aoff=1.2, along=.7  produces 1:2 ratio ellipses of -1 slope

% Redo with center of grid at 0,0
XO = -30000;
YO = -30000;
theta = 0;
Ampoff = 1.0;
Ampalong = 0;
Phase=90;


ct = cos(theta);
st = sin(theta);

bcs = [];
[nb,mb]=size(bndele);
for i = 1:nb
if bndele(i,5) == 5
% open boundary 

x1 = meshnod(bndele(i,2),2);
x2 = meshnod(bndele(i,3),2);
y1 = meshnod(bndele(i,2),3);
y2 = meshnod(bndele(i,3),3);

% Kelvin wave solution for exp decay of height off shore
      TIDEP = 12.42*3600;
      FFREQ = 2.*pi/TIDEP;
      Ho =  ( bathy(bndele(i,2),2)+ bathy(bndele(i,3),2) )/2.;
Ho = max(8,Ho);
      OMEGA2 = 1.4584E-4 ;
      RAD = 6.3675E6 ;
      DEG = 31.9333;
      f = OMEGA2*sin(DEG/180.*pi); 
% disp(['coriolis = ',num2str(f),'  s-1'])
      Co = sqrt(9.8*Ho);
      Ro = Co/f  ;  
%  disp(['Rossby Radius = ',num2str(Ro),'  m'])
%  ffreq = Co k  
      k = FFREQ/Co;


% Rotate the shore line of the Kelvin wave
dah1 = (x1-XO)*ct - (y1-YO)*st;     % Distance along shore
doff1 = (x1-XO)*st + (y1-YO)*ct;    % Distance off shore
dah2 = (x2-XO)*ct - (y2-YO)*st;     % Distance along shore
doff2 = (x2-XO)*st + (y2-YO)*ct;    % Distance off shore




A  = Ampoff*exp(-doff1/Ro);
ap = k*dah1 ;
B  = Ampalong*exp(-dah1/Ro);
bp = k*doff1 + Phase*pi/180;
phase1 = -atan2((A*sin(ap)+B*sin(bp)), (A*cos(ap)+B*cos(bp)))* 180/pi;
amp1   = sqrt( (A*cos(ap)+B*cos(bp))^2 + (A*sin(ap)+B*sin(bp))^2);
% disp([int2str(dah1),'  ',int2str( doff1),' : ',...
%                 num2str(amp1),'  ',num2str(phase1)])


A  = Ampoff*exp(-doff2/Ro);
ap = k*dah2 ;
B  = Ampalong*exp(-dah2/Ro);
bp = k*doff2 + Phase*pi/180;
phase2 = -atan2((A*sin(ap)+B*sin(bp)), (A*cos(ap)+B*cos(bp)))* 180/pi;
amp2   = sqrt( (A*cos(ap)+B*cos(bp))^2 + (A*sin(ap)+B*sin(bp))^2);



bcs = [bcs ; [bndele(i,2) amp1 phase1]];
bcs = [bcs ; [bndele(i,3) amp2 phase2]];
%disp([bndele(i,2) x1 bndele(i,3) x2 amp phase1]);
end

end

k = 1;
for i=2:length(bcs)
if bcs(i)~=bcs(k)
k=k+1;
bcs(k,:)=bcs(i,:);
end
end
bcs= bcs(1:k,:);

l = 1:k;
clf
subplot(2,2,1)
plot3prj(meshnod(bcs(l,1),2)/1000,meshnod(bcs(l,1),3)/1000,bcs(l,2))
title(' amplitude')
zlabel('meters')
xlabel(' W-E kilometers ')
ylabel(' S-N kilometers ')

subplot(2,2,2)
plot3prj(meshnod(bcs(l,1),2)/1000,meshnod(bcs(l,1),3)/1000,bcs(l,3))
title(' Phase')
zlabel('Degrees')
xlabel(' W-E kilometers ')
ylabel(' S-N kilometers ')

% Rotate the shore line of the Kelvin wave
dah = (meshnod(bcs(l,1),2)-XO)*ct - (meshnod(bcs(l,1),3)-YO)*st;     % Distance along shore
doff = (meshnod(bcs(l,1),2)-XO)*st + (meshnod(bcs(l,1),3)-YO)*ct;    % Distance off shore

subplot(2,2,3)
plot3prj(dah/1000,doff/1000,bcs(l,2))
title(['Ampoff=',num2str(Ampoff),'  Ampalong=',num2str(Ampalong),...
   ' + ',int2str(Phase),' deg'])
zlabel('meters')
xlabel(' alongshore km ')
ylabel(' offshore km ')

subplot(2,2,4)
plot3prj(dah/1000,doff/1000,bcs(l',3))
title(['Origin(',int2str(XO),',',int2str(YO),...
        ')  theta=',num2str(theta),' deg'] )
zlabel('Degrees')
xlabel(' alongshore km ')
ylabel(' offshore km ')

