% routine to calculate gradients in velocity given 
% a time series and geographic (lat/lon) position for each record

%ph=menu('What phase?','calulate gradients','calulate means',...
%          'plot gradient vectors','calculate overall stats');
ph=4;

if ph==1;
  % load time series with matrix "t" and "y"
  load /data5/jmanning/moor/moor82/wcrpl33
  stan=[1 1 2 2 3 3 3 3 3 3 5 5 6 6 6 6 6 6 7 7];
  depn=[8 8 8 8 40 40 100 100 305 305 8 8 8 8 35 35 84 84 8 8] 

  % load positions
  load /home3/ocn/jmanning/wcr/position/position_all.dat
  lat=position_all(:,7);
  lon=position_all(:,8);
  sta=position_all(:,1);
  dep=position_all(:,5);

  % loop through combinations of matrix y & calulate gradients
  % with index as follows i=vel.comp, j=var#1, and k=var#2
  [nr,nc]=size(y);

  % first do "u" component and then "v"
  % (note: we still call matrix "dudx" and "dudy" w/i=2 eventhough it is dy)

  for i=1:2;
    jj=0;
    for j=i:2:nc;
      jj=jj+1;
      kk=0;
      for k=i:2:nc;
        kk=kk+1;
        id1=find(sta==stan(j));
        id2=find(sta==stan(k));

        % calculate vertical gradients if same site & depn(k)>depn(j)
        if id1(1)==id2(1);
          distz=depn(j)-depn(k);
           % the sign should be positive when comparing shallow - deep since
           % up is positive
          distz=-1*distz;
          dudz(jj,kk)=mean((y(:,j)-y(:,k))./(distz*100*ones(nr,1)));
        end;

        %     calculate dist between stations
        distx=dist(lat(id1(1)),lon(id1(1)),lat(id1(1)),lon(id2(1)))*100000.;
        %     if 1st pos'n west of 2nd, reverse sign 
        if lon(id1(1))<lon(id2(1));distx=-1*distx;end;
        disty=dist(lat(id1(1)),lon(id1(1)),lat(id2(1)),lon(id1(1)))*100000.;
        %     if 1st pos'n south of 2nd, reverse sign 
        if lat(id1(1))<lat(id2(1));dist=-1*dist;end;
        if distx~=0;dudx(jj,kk)=mean((y(:,j)-y(:,k))./(distx*ones(nr,1)));end;
        if disty~=0;dudy(jj,kk)=mean((y(:,j)-y(:,k))./(disty*ones(nr,1)));end;
      end;
    end;
    if i==1;
      fploopy(dudx,'dxdy.out','o')
      fploopy(dudy,'dxdy.out','a')
      fploopy(dudz,'dxdy.out','a')
    else;
      fploopy(dudx,'dxdy.out','a')
      fploopy(dudy,'dxdy.out','a')
      fploopy(dudz,'dxdy.out','a')
    end;
  end;
end;

if ph==2;
  % phase 2: calculating means
  load /home3/ocn/jmanning/wcr/dxdy/dxdydz.dat
  d=dxdydz;
  for k=1:6;
    for i=1:10;
      dd=d(k*10-9:k*10,i);
      id=find(dd~=0&isnan(dd)~=1);
      m(k,i)=meannnan(dd);
      mi(k,i)=minnnan(dd);
      ma(k,i)=maxnnan(dd);
    end;      
  end;
  y=[m mi ma]'
  save  /home3/ocn/jmanning/wcr/dxdy/dxdydzm.out y -ascii
end;

if ph==3;
  %  plot gradients

  % load gradients
  load /home3/ocn/jmanning/wcr/dxdy/dxdydz.dat
  d=dxdydz;
  % load positions
  load /home3/ocn/jmanning/wcr/position/position_all.dat
  lat=position_all(:,7);
  lon=position_all(:,8);
  sta=position_all(:,1);
  dep=position_all(:,5);
  stan=[1 1 2 2 3 3 3 3 3 3 5 5 6 6 6 6 6 6 7 7];
  depn=[8 8 8 8 40 40 100 100 305 305 8 8 8 8 35 35 84 84 8 8] 

  %du/dxdy
  figure(1)
  for i=1:10;
    % find mid-point
    for j=1:10;
      latm(i,j)=mean([lat(i) lat(j)]);
      lonm(i,j)=mean([lon(i) lon(j)]);
    end;
    subplot(5,2,i);
    % start plot
%    plot([-71 -65],[38 42],'y+')
    pause(1);
%    hold on;
    pltwcr(2)
    pause(1);
    vecplot2(lonm(i,:),latm(i,:),d(1:10,i),d(11:20,i),.00001,...
        '/s',-75,35,'g')
    title([num2str(stan(i*2)),'@',num2str(depn(i*2)),'m'])
    axis([-69 -67 40 41.5])
    if rem(i,2)==0;set(gca,'yticklabels','');end;
    if i<9;set(gca,'xticklabels','');end;
  end;
  
  %dv/dxdy
  figure(2)
  for i=21:30;
    % find mid-point
    for j=1:10;
      latm(i-20,j)=mean([lat(i-20) lat(j)]);
      lonm(i-20,j)=mean([lon(i-20) lon(j)]);
    end;
    subplot(5,2,i-20);
    % start plot
    pltwcr(2)
    pause(1);
    vecplot2(lonm(i-20,:),latm(i-20,:),d(21:30,i-20),d(31:40,i-20),.001,...
        '/s',-75,35,'g')
    title([num2str(stan((i-20)*2)),'@',num2str(depn((i-20)*2)),'m'])
    axis([-69 -67 40 41.5])
    if rem(i,2)==0;set(gca,'yticklabels','');end;
    if i<29;set(gca,'xticklabels','');end;
  end;
end;
  
if ph==4; % calculate overall stats
  load /home3/ocn/jmanning/wcr/dxdy/dxdydzm.out
  d=dxdydzm;
  for k=1:6;
    meanshr(k)=mean(d(1:3:length(d),k));
    stdshr(k)=std(d(1:3:length(d),k));
    meanabs(k)=mean(abs(d(1:3:length(d),k)));
    minshr(k)=min(d(2:3:length(d),k));
    maxshr(k)=max(d(3:3:length(d),k));
  end;
  y=[meanshr; stdshr; meanabs; minshr; maxshr];
  save /home3/ocn/jmanning/wcr/dxdy/dxdydzst.dat y -ascii
end;
