%
% Script file to extract the tidal and steady current fields
% from a ship-board ADCP survey of a given region.
%
% for details see: Candela et al., 1992, "Separation of tidal and ..."
%                  JGR, 97 C1, 769-788.
%
% for MATLAB 4.0,      Julio Candela, VII/1993. (WHOI).
%

% if metros=1, no scaling is used
% prt = 1, for hardcopy plots.
% vel = 1, to fit vertically averaged velocity instead of transport.
% nf = # of tidal constituents to consider. You can use up to 2 (M2 and S2)
%      or 4 (M2,S2,K1 and O1), however you'll have to augment arrays "vu, fan"
%      with the proper numbers for these diurnal constituents using program
%      gphase.f
% pm = polynomial degree for steady part.
% pt = polynomial degree for each tidal constituent to consider.
% icontran = 1, to draw contours of tidal ellipse components instead of axes.

 metros=0;
 iprt=1; 
 vel=1;
 nf=1; pm = 1; pt=[2,0]; nt(1:4)=zeros(1,4); %
 icolor=0; 
tree='/vol/mar/d/georges/sjohnson/';
cruises=['sj9503';'sj9504';'sj9505';'sj9506';'sj9507';'sj9508'];
 cruise='sj9505';

%for ncru=1:6
%cruise = cruises(ncru,:)

%
% phase of the tidal potential analysed constituents w.r.t. the
% Greenwich meridian. vu = vu /(2*pi),
% include also amplitude modulation fan,
% Use the fortran program "gphase.f" to calculate these numbers
% for a different time reference. If you use the numbers given below
% your data should have as time reference 0h , 0d, Jan, 1994.
%
% (M2,S2) for time reference (0,0,1,93) [i.e., (h,d,m,y) UT]
% 
% vu = [-.4383293, -.000327218];% vu = [0,0] % for no Greewich phase correction.
%
% (M2,S2) for time reference (0,0,1,94) [i.e., (h,d,m,y) UT]
% vu = [-1.159125, -.0002829144]; 
%
% (M2,S2) for time reference (0,1,1,94) [i.e., (h,d,m,y) UT]
vu = [-1.122682, -.0002636539];
% 
%  vu = [0,0]; for no Greewich phase correction.
  

fan = [1. , 1.];

% constants
% define scales and labels for plots
% scales 

if vel,
 scale= 7.e-1; stext='50 cm/s'; svec=.50; 
 seli= 3.e-1; stexeli='1.0 m/s'; sveli=.50; 
 sscale= 3.e-0; sstext='20 cm/s'; ssvec=.20;
else 
 scale= 8.e-3; stext='50 m^2/s'; svec=50.; 
 seli= 5.e-3; stexeli='100 m^2/s'; sveli=50.; 
 sscale= 4.e-2; sstext='20 m^2/s'; ssvec=20.;
end

% xymax = longitude and latitude range of plots.
% xloref = reference longitude.
% ylaref = reference latitude.
% xloesc = length in 1.e-6 m of a minute in longitude at the mean latitude
%          of the surveyed area.
% ylaesc = length in 1.e-6 m of a minute in latitude at ...
%

  xymax=[-71 -65. 39 45];
  xloref = -71. ; ylaref = 39;
  ylaesc = .111052 ; xloesc = .084137 ;

% load ADCP measured transports

t=[];x=[];y=[];h=[];u=[];v=[];

%% load data 
%
% data file should consist of the columns [t,x,y,u,v,h];
% where:
% t = time (in hours) of ADCP averaged observation according to the time 
%     reference chosen.
% x = scaled x coordinate in 1.e-6 meters from "xloref".
% y = scaled y coordinate in 1.e-6 meters from "ylaref".
% u = East component of observed transport in [m^2/s].
% v = North component of transport in [m^2/s].
% h = local water depth recorded by ADCP and averaged the same as the current
%     profiles to obtain u and v.

% file = ['tran',num2str(60),'c'];	
% eval(['load ',file,'.dat']) ; D = eval([file]); eval(['clear ',file]);
% t=[t;D(:,2)]; x=[x;D(:,3)]; y=[y;D(:,4)];
% u=[u;D(:,8)]; v=[v;D(:,9)]; h=[h;D(:,7)]; clear D;

 file = [tree,cruise,'/tranh',cruise];	
 eval(['load ',file]) ;
 t=[t;th']; x=[x;lonh']; y=[y;lath'];
 u=[u;uth']; v=[v;vth']; h=[h;deph']; 

xg=x; yg=y;

 ntran=length(t);

% plot data locations
% georbat is a script file that plots coastlines and bathymetric contours
% of the surveyed area and holds the screen.

 figure (1);clf, gbmap,%georbat;
 set(gca,'xlim',[-71 -64.5],'ylim',[39.5 43.2],'position',[.15 .13 .7 .6806])
 hp=plot(xg,yg,'g.');set(hp,'markersize',11)
 title(['ADCP hourly data locations from Seward Johnson Cruise ',cruise])
 xlabel('Longitude'),  ylabel(['Latitude, #data = ',int2str(ntran)]),grid
 hold off; pause(3)
  if iprt, orient landscape,  
 if icolor, eval(['print -dpsc ',cruise,'dat']),
     else, eval(['print ',cruise,'dat']), end, end

% plot raw data vectors 
 
figure(2),clf, gbmap,%georbat;
if vel,title(['Vertically averaged hourly current vectors, Cruise ',cruise]),
 else,
title(['Observed hourly average transport vectors']), end
  xlabel('Longitude'),  ylabel(['Latitude, #data = ',int2str(ntran)])
 set(gca,'xlim',[-71 -64.5],'ylim',[39.5 43.2],'position',[.15 .13 .7 .6806])
 xyaxis=[-71 -64.5 39.5 43.2]; gsca = [11 8.5 .7 .6806];
 hold on
 pltvecgb(xg,yg,u,v,scale,xyaxis,'g',gsca)
 scavecgb(-65.5, 40,svec,stext,scale,xyaxis,'r',gsca)
 hold off 
 pause(3)
 if iprt, orient landscape,  
 if icolor, eval(['print -dpsc ',cruise,'raw']),
     else, eval(['print ',cruise,'raw']), end, end

% Construct model matrix elements

 Mm = matpold(x,y,pm,0);  [mm,nm]=size(Mm);
 
   for k=1:nf
      s = ['Mt',num2str(k)];
      eval([s,'= matpold(x,y,pt(k),0);']);
      ss = ['size(',s,')']; [mt(k),nt(k)]=eval(ss); 
   end

%
% construct matrix chosing the number (nf)
% of tidal frequencies (f) (in cyc/h) to include in the analysis.
% freq = [m2=.08051140, s2=.08333333, k1=.04178075, o1=.03873085]
%
 freq = [.08051140 .08333333 .04178075 .03873085];
 Tt = zeros(ntran,nf);

 for k=1:nf, w(k) = 2*pi*freq(k); 
  eval(['tt',num2str(k),'c=cos(w(k).*t);']);
  eval(['tt',num2str(k),'s=sin(w(k).*t);']);

  eval(['tt',num2str(k),'c=tt',num2str(k),'c(:,ones(1,nt(',num2str(k),')));']);
  eval(['tt',num2str(k),'s=tt',num2str(k),'s(:,ones(1,nt(',num2str(k),')));']);

  eval(['z',num2str(k),'=zeros(size(Mt',num2str(k),'));']);

% convert frequency to rad/sec
% w(k) = w(k)/3600;
 end

% Construct model matrix

 zm = zeros(size(Mm)); 

% u part

 a = [Mm;zm];
for k=1:nf
 eval(['tema=[Mt',num2str(k),'.*tt',num2str(k),'c,Mt',num2str(k), ...
'.*tt',num2str(k),'s;z',num2str(k),',z',num2str(k),'];']);
 a = [a,tema];
end

% v part

tema = [zm;Mm]; a = [a,tema];
for k=1:nf
 eval(['tema=[z',num2str(k),',z',num2str(k),';Mt',num2str(k), ...
'.*tt',num2str(k),'c,Mt',num2str(k),'.*tt',num2str(k),'s];']);
 a = [a,tema];
end


[ma,na]=size(a),

% construct right hand side vector b
 
   b = [ u; v];

% solve system and plot results

 tol = 0.; [c,C,bm,cor,X,r,cond] = solsys(a,b,tol);cor
 res = b - bm;
 figure(3), clf 
 plot(b,bm,'o'); grid; title('b vs. bm'); pause(3)

% define indexes for coefficients
% "you can learn more about the model matrix structure by trying to write
%  these indexes into a more compact form"

 mum=1:nm; 
 mvm=nm+sum(2*nt(1:nf))+1:2*nm+sum(2*nt(1:nf));

 mut1c=nm+1:nm+nt(1);
 mut1s=nm+nt(1)+1:nm+2*nt(1);
 mut2c=nm+2*nt(1)+1:nm+2*nt(1)+nt(2);
 mut2s=nm+2*nt(1)+nt(2)+1:nm+2*nt(1)+2*nt(2);

 mut3c=nm+2*sum(nt(1:2))+1:nm+2.*sum(nt(1:2))+nt(3);
 mut3s=nm+2*sum(nt(1:2))+nt(3)+1:nm+2*sum(nt(1:3));
 mut4c=nm+2*sum(nt(1:3))+1:nm+2.*sum(nt(1:3))+nt(4);
 mut4s=nm+2*sum(nt(1:3))+nt(4)+1:nm+2*sum(nt(1:4));

 mvt1c=2*nm+2*nt(1)+2*nt(2)+1:2*nm+3*nt(1)+2*nt(2);
 mvt1s=2*nm+3*nt(1)+2*nt(2)+1:2*nm+4*nt(1)+2*nt(2); 
 mvt2c=2*nm+4*nt(1)+2*nt(2)+1:2*nm+4*nt(1)+3*nt(2);
 mvt2s=2*nm+4*nt(1)+3*nt(2)+1:2*nm+4*nt(1)+4*nt(2);

 mvt3c=2*nm+4*sum(nt(1:2))+2*nt(3)+1:2*nm+4*sum(nt(1:2))+3*nt(3);
 mvt3s=2*nm+4*sum(nt(1:2))+3*nt(3)+1:2*nm+4*sum(nt(1:3));
 mvt4c=2*nm+4*sum(nt(1:3))+2*nt(4)+1:2*nm+4*sum(nt(1:3))+3*nt(4);
 mvt4s=2*nm+4*sum(nt(1:3))+3*nt(4)+1:2*nm+4*sum(nt(1:4));

%
% code contributed by Lyn Harris

mum=1:nm;                                        % Steady coefficients.
mvm=mum+na/2.;

mut1c=(1:nt(1)) + mum(nm);                       % Tidal coefficients.
mut1s=(1:nt(1)) + mut1c(nt(1));
mvt1c=(1:nt(1)) + mvm(nm);
mvt1s=(1:nt(1)) + mvt1c(nt(1));

if nf > 1
  mut2c=(1:nt(2)) + mut1s(nt(1));
  mut2s=(1:nt(2)) + mut2c(nt(2));
  mvt2c=(1:nt(2)) + mvt1s(nt(1));
  mvt2s=(1:nt(2)) + mvt2c(nt(2));
  if nf > 2
    mut3c=(1:nt(3)) + mut2s(nt(2));
    mut3s=(1:nt(3)) + mut3c(nt(3));
    mvt3c=(1:nt(3)) + mvt2s(nt(2));
    mvt3s=(1:nt(3)) + mvt3c(nt(3));
    if nf > 3
      mut4c=(1:nt(4)) + mut3s(nt(3));
      mut4s=(1:nt(4)) + mut4c(nt(4));
      mvt4c=(1:nt(4)) + mvt3s(nt(3));
      mvt4s=(1:nt(4)) + mvt4c(nt(4));
    end
  end
end

% save coefficients for detiding each cruise separetly

cut = [c(mut1c);c(mut1s)]; cvt = [c(mvt1c);c(mvt1s)];
if nf==2,
cut =[cut;c(mut2c);c(mut2s)]; cvt = [cvt;c(mvt2c);c(mvt2s)];end
if nf==3,
cut =[cut;c(mut3c);c(mut3s)]; cvt = [cvt;c(mvt3c);c(mvt3s)];end
if nf==4,
cut =[cut;c(mut4c);c(mut4s)]; cvt = [cvt;c(mvt4c);c(mvt4s)];end

save ctotal cut cvt
%

% plot transport cotidal charts

% blanking polygon
% poligono.mat contains the coordinates (geographical) of the vertices
% [xpoly,ypoly] of a poligon enclosing the ADCP observation points for
% ploting purposes. It can be generated from a MATLAB window having a plot
% of the surveyed area with the comand: "[xpoly,ypoly]=ginput" ; chosing the
% vertices with the mouse and then saving them: "save poligono xpoly ypoly"

 ires=.2; xel=-70:ires:-65.4; yel=40:ires:42.8;
 [xxe,yye]=meshgrid(xel,yel); mme=length(xxe(:));
 zze=zeros(size(xxe));
load poligono
%xpoly=(xpoly-xloref.*ones(size(xpoly))).*xloesc;
% ypoly=(ypoly-ylaref.*ones(size(ypoly))).*ylaesc;
%if metros, xpoly=1.e6.*xpoly; ypoly=1.e6.*ypoly; end
Bnan = inside(xxe,yye,xpoly,ypoly);
% xelg = (xxe(:)./xloesc) + xloref.*ones(size(xxe(:)));
% yelg = (yye(:)./ylaesc) + ylaref.*ones(size(yye(:)));

for k=1:nf;

 DE=matpold(xxe(:),yye(:),pt(k),0);

% tidal ellipse chart
% plot tidal ellipse axis at specified points (xp,yp)

 mutc = eval(['mut',num2str(k),'c']); muts = eval(['mut',num2str(k),'s']);
 mvtc = eval(['mvt',num2str(k),'c']); mvts = eval(['mvt',num2str(k),'s']);

 cx =  DE*real(c(mutc)); sx = DE*real(c(muts));
 cy =  DE*real(c(mvtc)); sy = DE*real(c(mvts));
 au = sqrt(cx.^2+sx.^2); fu = atan2(sx,cx);
 av = sqrt(cy.^2+sy.^2); fv = atan2(sy,cy);
 wp = .5.*(au.*exp(-i.*fu) + i.*av.*exp(-i.*fv));
 wm = .5.*(au.*exp(i.*fu) + i.*av.*exp(i.*fv));
 awp = abs(wp) / fan(k);
 awm = abs(wm) / fan(k);
 gwp = (unwrap(angle(wp))./(2*pi)) - vu(k); gwp = 360*(gwp - fix(gwp));
 gwm = (unwrap(angle(wm))./(2*pi)) + vu(k); gwm = 360*(gwm - fix(gwm));
 orip = (gwm+gwp)/2;
 fase = (gwm-gwp)/2;

% calculate fase for flooding tide, i.e. into Mass bay.
 nori = find(orip<=30 & orip>-150); orip(nori)=orip(nori)+180;
 fase(nori) = fase(nori)+180.;
 ipha = find(fase<0); fase(ipha)=fase(ipha)+360;
 ipha = find(fase>360); fase(ipha)=fase(ipha)-360;

figure(4),clf, gbmap,%georbat;
%
if vel,
% title(['Vertically average tidal current ellipse constituent ', num2str(freq(k)) ,' (cyc/h)']),
title(['M2 barotropic current ellipses, Cruise ',cruise])
else
 title(['Tidal transport ellipse constituent ',num2str(freq(k)),' (cyc/h)']), end
%  xlabel([file,', pm=',num2str(pm), ...
%' pt(',num2str(pt(1)),',',num2str(pt(2)),'), #data = ',num2str(ntran)])
ylabel('Latitude'); xlabel('Longitude');% ylabel(['corr= ',num2str(cor)])
 set(gca,'xlim',[-71 -64.5],'ylim',[39.5 43.2],'position',[.15 .13 .7 .6806])
 xyaxis=[-71 -64.5 39.5 43.2]; 
 hold on
pltelip(xxe(:).*Bnan(:),yye(:),(awp+awm),awp-awm,seli,xyaxis,orip,'c',gsca)
 scapax4(-65.5,40,sveli,stexeli,seli,xymax,'r',gsca)
 hold off, pause(3)
  if iprt, orient landscape,  
 if icolor, eval(['print -dpsc ',cruise,'elip']),
     else, eval(['print ',cruise,'elip']), end, end

% plot phase of ellipses

ifase=1,
if ifase, vpha=0:30:360; zz=zeros(size(xxe));
 zz(:) = fase.*Bnan(:);
 figure(5), clf,georbat;
 set(gca,'xlim',[-71 -64.5],'ylim',[39.5 43.2],'position',[.15 .13 .7 .6806])
 cs=contour(xxe,yye,zz,vpha); clabel(cs); %; clabel(cs,'manual');
%contour(xxe,yye,zz,vpha); 
 title(['Phase of ellipse (UT) frequency ',num2str(freq(k)),' (cyc/h)']),
  xlabel([file,', pm=',num2str(pm), ...
' pt(',num2str(pt(1)),',',num2str(pt(2)),'), #data = ',num2str(ntran)])
 ylabel(['corr= ',num2str(cor)])
 hold off; pause (3)
  if iprt, orient landscape,  
 if icolor, eval(['print -dpsc ',cruise,'pha']),
     else, eval(['print ',cruise,'pha']), end, end

end

end ,%end for loop


% plot steady field

steady=0;
if steady

 DE=matpold(xxe(:),yye(:),pm,0);
 u0 = DE*c(mum); v0 = DE*c(mvm);
 figure(6), clf; gbmap;
if vel, title(['Steady vertically averaged currents, Cruise ',cruise]), else
title(['Steady transports']), end
ylabel('Latitude'); xlabel('Longitude');
%  xlabel([file,', pm=',num2str(pm), ...
%' pt(',num2str(pt(1)),',',num2str(pt(2)),'), #data = ',num2str(ntran)])
% ylabel(['corr= ',num2str(cor)])
 set(gca,'xlim',[-71 -64.5],'ylim',[39.5 43.2],'position',[.15 .13 .7 .6806])
 xyaxis=[-71 -64.5 39.5 43.2]; gsca = [11 8.5 .7 .6806];
 hold on
 set(gca,'xlim',[-71 -64.5],'ylim',[39.5 43.2],'position',[.15 .13 .7 .6806])
 xyaxis=[-71 -64.5 39.5 43.2];gsca = [11 8.5 .7 .6806];
 hold on
%
% gau = gau2d(xxe(:),yye(:));
% u0 = u0.*gau; v0=v0.*gau;
%
 pltvecgb(xxe(:).*Bnan(:),yye(:),u0,v0,sscale,xyaxis,'g',gsca)
 scavecgb(-65.5, 40,ssvec,sstext,sscale,xyaxis,'g',gsca)
 hold off 
 pause(1)
  if iprt, orient landscape,  
 if icolor, eval(['print -dpsc ',cruise,'sub']),
     else, eval(['print ',cruise,'sub']), end, end
end
%save variables for later use
save steady pm c mum mvm pt nf mutc mvtc muts mvts w

end ,%ncru
