%
%             "foi"
%
%    ...  Fisheries Oceanography Investigation's
%
%         MATLAB DRIVER PROGRAM
%         TO CONDUCT ANALYSIS
%         OF TIME SERIES DATA
%
%
% on user's first pass through the routine, display a message
%
e=exist('first');
if e==0;
 disp(' ');
 disp(' !! PLEASE TAKE NOTES OF PROBLEMS ENCOUNTERED !!');
 disp(' ');
end;%
first=1;
%
% This routine often calls routines in Jim's "/atsee" directory
%

cp=computer;
path(path,'/home3/ocn/jmanning/atsee/stick');
%
%
%                    Main Menu

todo=menu(' Main Menu ',...
          '-1) quit   ',...
          ' 0) load data into MATLAB from other formats', ...
          ' 1) time series plots    ', ...
         '  2) filter or rotate    ', ...
         '  3) autospec  ', ...
         '  4) cross-spec   ', ...
         '  5) eofs  ', ...
         '  6) stats [general (mean,min,etc.), correlation matrix, histogram]  ', ...
         '  7) enter any MATLAB command (keyboard)', ...
         '  8) view modeled flow fields (on SUN machine only)', ...
         '  9) process CTD data',...
         '  10) principal axis');
%
% alter the value of "todo".
%
todo=todo-2;
if todo==7;
  todo=99;
end;
%
if todo==-1;
  return;
end;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%% LOAD DATA  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if todo == 0

%  note: when finished this phase of foi, one should have two variables:
%        t=yrday0_gmt
%        y=matrix with multiple colunms of variables
%  we get the data with a subroutine "gettsd.m"

   gettsd;
end
%
%%%%%%%%%%%%%%%%%%%%%   TIME SERIES PLOT   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if todo == 1
%
  

% Next, pickout particular yearday desired
  pick=menu('Which do you want?','all the data','particular subset');
  if pick==2;
    st=input('Enter starting yearday (note: noon Jan1 = 0.5): ');
    et=input('Enter ending yearday   (note: noon Jan1 = 0.5): ');
  else;
    st=0;
    et=365;
  end;

  exy=exist('yy');
  if exy==0

     yy=input('Enter a two-digit starting year (88, for example, for 1988): ');
  end;
           
  ff=find(t>=st&t<=et);


% Next choose between stick plot or progressive vector
  wtp=menu('Which type of plot?','stickplot','progressive vector')
  
  if wtp==1;
% 
    nsp=input(' Enter # of stick plots: ');
    tnp=input(' Enter total # of plots: ');
    cya=menu('How do you want y-axis scale set?',...
             'same size for each plot',...
             'default constant',...
             'user specified (prompted for each plot)',...
             'read from a yax.dat file');

    cyl=menu('How do you want specify labels?',...
             'user specified (prompted for each plot)',...
             'read from a yaxl.dat file');
    if cyl==1;
      for kk=1:tnp;
        label(kk,:)=input(['Enter label #',num2str(kk),': '],'s');
      end;
    elseif cyl==2;
      flid=fopen('yaxl.dat');
      for kk=1:tnp;
        label(kk,1:14)=fgetl(flid);
      end;
    end;

    anom=input('Enter level at which to plot red to indicate the strongest: ');

    nth=input('plotting every nth point ... enter nth: ');
    stkcombo(yy,t(ff),y(ff,:),nsp,tnp,cya,label,anom,nth)

  elseif wtp==2;
    wop=menu('Annotation?','about 5 times','mthly','user specied','none');
    if wop==1;
      disp('Hardcoded to annotate about 5 times')
     tint2=range(t(ff))/5*3600*24;
    elseif wop==2;
      tint2=30*24*3600;
    elseif wop==3;
      tint2=input('Enter annotation interval in seconds: ');
    else;
      tint2=0;
    end;
%   decide scaling function to get cm/sec
%    if ek==2;
      sc=.1;
%    else;
%      sc=100;
%    end;
    disp('Warning: "provecxy" assumes cm/s values');
    disp('Warning: scaling factor set ');

    
%   make subplots if #col in y >2
    [nr,nc]=size(y);
    sp=0;

%   get positions of these series

    cp=computer;

    if cp(1:3)=='SGI';
       load /home3/ocn/jmanning/wcr/position/position_all.dat

    else;

       load e:\wcr\position\position_all.dat

    end;   
    p=position_all;
 %   ps=[1 2 6 7 10]

     lcol=['r','g','k','m','b']
   ps=[3 4 5 8 9] 
    for i=1:2:nc-1;
      sp=sp+1;
%      if nc>2;
%        subplot(nc/2,1,sp)   
%      end;
      if sp==1;georbat(1,0,0,2,1);end;
       provecxy(y(ff,i)*sc,y(ff,i+1)*sc,t(ff),tint2,p(ps(sp),7),p(ps(sp),8),lcol(sp));
      hold on;
    end;
%   axis('normal')
%   axis([-1000 200 -400 100]);
    axis([-69 -67.5 40 41]);

    if cp(1:3)=='PCW';

       cd c:\matlab\toolbox\matlab\graph2d

    end;

%     legend('11','21','51','61','71',-1)
    legend('32','33','34','62','63',-1)
%      
    lllabel(30,30,40,42,-70,-67)
%    title('Surface Advection (10m)','fontsize',12);
    title('Deeper Advection (>30m)','fontsize',12);
    if sp~=nc/2;
      set(gca,'xticklabels',' ');xlabel('');ylabel('');
    end;
  end;
  orient tall

end
%
%
%%%%%%%%%%%%%%%%%%%%%   FILTERS            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if todo == 2
   typefilt=menu(' What do you want?','rotate U & V',...
                                              'pl33',...
                                       'butterworth',...
                                       'running ave',...
                                       'subsample',...
                                       ' '    );
%%%%%%%%   rotate
%
   if typefilt==1
     angl=input('Enter the angle to rotate (positive clockwise from north): ');
     [y(:,1),y(:,2)]=rot(angl,y(:,1),y(:,2));
   end
%
%%%%%%%    PL33
%
   if typefilt==2

%   note: this option will filter all coluumns of matrix "y"
%         and output a new matrix yf

      
%   now filter

      for kk=1:hm;
        yf(:,kk)=plfilt(y(:,kk));
        plot(t,y(:,kk),t,yf(:,kk));
        xlabel('days');
        disp(' ');
        disp(['Variable # ',num2str(kk),'.  Press return to continue.']);
        pause
      end
   end;
%
%%%%%%%%  Butterworth filter
%
   if typefilt==3
      hm=input(' How many variables do you want to filter? ');
      minp=input(' Enter the minimum period of interest in hours: ');
%
% if sampling interval is not known it is asked for
%
     if ti==[];
       ti=input('What is your sampling interval in hours?');
     end;
      nyqf=(24/ti)/2;
      cof=(24/minp)/nyqf;
%      ord=menu(' Order Butterworth filter?','second','third','forth','none');
%      ord=ord+1;
      ord=3;
      for ih=1:hm;
      wv=input(' Which variable? (integer order in which they were loaded) ');
         [b,a]=butter(ord,cof);
         yf(:,wv)=filtfilt(b,a,y(:,wv));
         plot(t,y(:,wv),t,yf(:,wv));
         disp(' ');
         disp(' press return to continue')
         pause;
         clf
       end;
    end
%
%%%%%%%%  Running average filter
%
  if typefilt==4
       n=input('# of points to average?: ');
       disp('')
       disp('Processing..');
       [yf,tf]=runave(t,y,n);
  end;

%
%%%%%%%%  Subsample
%
  if typefilt==5
%   check to make sure we have hourly data
%   and if not interpolate to such
%    [rg,mi,ma]=rangejm(diff(t));
   
%    if mi~=ma|rg~=0; 
%     make the first value "on-the-hour"
      st=(t(1)*24-rem(t(1)*24,1))/24+1/24;
      tf=[st:1/24:max(t)];

%     interp to hourly values from running ave file
      [hr,hm]=size(y);
      disp('interpolating to one hour...')
      for kk=1:hm;
        yf(:,kk)=interp1(t,y(:,kk),tf);
      end;
      y=yf;
      t=tf;
%    end;
  end;
%

%
  if typefilt > 1
     pltyn=menu('What Next? ',...
           ' hardcopy',...
           ' ascii store filtered set with "t" as 1st col',...
           ' save tf & yf as new t & y in mat file',...
           'continue');
     if pltyn == 1
       print -deps yfplt
     elseif pltyn == 2
       y=[tf yf];
       fploopy(y);
     elseif pltyn == 3
       fn=input('Enter filename (include path but no spec): ','s');
       y=yf;
       t=tf';
       eval(['save ',fn,' t y']);
     end
  end
     hold off
end
%


%%%%%%%%%%%%%%%%%%%%%   AUTO SPECTRA       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if todo == 3
%
disp(' ');
disp(' We will step through spectra using different segment lengths' )
disp(' pausing five seconds for each case.')
disp(' ');
   pause(2); 
   %
 which=input('Which variable do you want to examine? 1? 2? ..n? ');
   %
   endpt=length(y(:,which));
   %
   % find the power of two that forms "endpt"
   %
   j=4;
   while 2^j < endpt, j=j+1; end
   %
   j=j-1;
   %
 wtype=menu('Which type of spectra?','power spectra','variance preserving');
   sampd=24/ti;
%
 numfb=input('Number of frequency bands to average?');
%
   for i=6:j
     npt=2^i
     ovl=npt/2;
       p=spectrum(y(:,which),npt,ovl);
       [n,m]=size(p);
       f=(1:n)/n*sampd/2;
       if numfb>1;
         [p2,f2]=runave(f,p,numfb);
       end;
%
       if wtype==1;
       spcpltjm2(p2,f2,1);
       elseif wtype==2;
       f=[f;f];
       p=p.*f';
       spcpltjm(p2,sampd/numfb);
       end;
       pause(1);
     %
   end
   cont=menu('Now what?','hardcopy best result','continue');
     if cont == 1;
   %
     npt=input(' what npt value? ')
     ovl=npt/2;
   %
     p=spectrum(y(:,which),npt,ovl);
       [n,m]=size(p);
       f=(1:n)/n*sampd/2;
       if numfb>1;
         [p2,f2]=runave(f,p,numfb);
       end;
%
       if wtype==1;
       spcpltjm2(p2,f2,1);
       elseif wtype==2;
       [n,m]=size(p);
       f=(1:n)/n*sampd/2;
       f=[f;f];
       p=p.*f';
       spcpltjm(p2,sampd/numfb);
       end;
       aspec=int2str(npt)
disp(' ');
disp(' Enter appropriate text and print commands, then type "return".')
disp(' ');
keyboard; 
      end
end
%
%%%%%%%%%%%%%%%%%%%%%   CROSS - SPECTRA    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if todo==4

hmp=input('How many pairs do you have to analyze?');
maxcypd=input('Enter maximum cycles per day to consider: ');
disp(['Maximum cycle per day: ',num2str(maxcypd)]);
disp(' ');
numfb=input('Number of frequency bands to average?');
skip=menu('What would you like?','try different segment lengths',...
                                 'user specified value');
if skip==2;
    npt=input(' What segment length (#pts) do you want? ')
end;

% set-up instr codes
icode=[001 002 111 112 211 212 321 322 331 332 341 342 ...
       511 512 611 612 621 622 631 632 711 712]

% get time in fractional days
ti=mean(diff(t))*24;

% loop through all combinations of pairs
for kp=1;
if hmp==1;hmp=2;end;
for ip=2:hmp;

% if hmp>3 assume an automatic generation of pairs
  if hmp<=3;
    disp(' ');
    disp('Two variables are needed for cross-spectra.');
    wv1=input('  Which one in your input? ');
    wv2=input('  Which one is your output? ');
  else;
%    wv1=kp*2-1;
    wv1=1;
    wv2=ip*2;
  end;
  if wv1~=wv2;
  cont=0;
  if skip==1;
   %
     endpt=length(y(:,wv1))
   %
   % find the power of two that forms "endpt"
   %
     j=4;
     while 2^j < endpt, j=j+1; end
   %
     j=j-1;
   %
     sampd=24/ti;
%
%    Try calculating spectrum for various powers of two
%
     for i=6:j
       npt=2^i
       ovl=npt/2;
%
%      Calculate the # of pieces (including overlap) used in this case
%
       nump=floor((length(y)/npt)*2-1)
%
%      Calculate spectrum
%
       p=spectrum(y(:,wv1),y(:,wv2),npt,ovl);
       [n,m]=size(p);
       f=(1:n)/n*sampd/2;
%
%      Average over "numfb" frequencies using MATLAB's mean
%
       if numfb>1;
         lf2=floor(n/numfb);
         f2=ones(lf2,1);
         for jr=1:lf2;
           [f2(jr)]=mean(f(jr*numfb-numfb+1:jr*numfb));
         end; 
 
         for jc=1:m;
           for jr=1:lf2;
             [p2(jr,jc)]=mean(p(jr*numfb-numfb+1:jr*numfb,jc));
           end; 
         end;
       else;
         p2=p;
         f2=f;
       end;
%
% The following statistical quantities are taken from a writeup
% by Marlene Noble for the Buoy Processing Programs
%
% Calculate phase confidence intervals (using the "stats toolbox")
%
       degf=floor(nump*numfb-2);
       st=tinv(.975,degf);
%
       [nr,nc]=size(p2);
       for j=1:nr;
       ciph(j)=180/pi*asin(st*sqrt((1-real(p2(j,5))^2)/(degf*real(p2(j,5))^2)));
       end;
%
% calculate the coherence C.I.
%
       degf=nump*numfb/2
       cic=sqrt(1-.05^(1/(degf-1)));

%
       spcpltjm2(p2,f2,maxcypd,ciph,cic,99);
       pause(3)
    end;
%
%

    cont=menu('Now what?','view cross-spec plots',...
                        'output results to ascii file',...
                        'both ',...
                        'continue');
  end;
%

  if cont < 4| skip==2;
    clear p p2 f2 ciph;
    sampd=24/ti;
    if hmp<3;
      npt=input(' What segment length (#pts) do you want? ')
    end;
%
%   Calculate the # of pieces (including overlap) used in this case
%
    nump=floor((length(y)/npt)*2-1)
    ovl=npt/2;
    p=spectrum(y(:,wv1),y(:,wv2),npt,ovl);
    [n,m]=size(p);
    f=(1:n)/n*sampd/2;
    if numfb>1;
      lf2=floor(n/numfb);
      f2=ones(lf2,1);
      for jr=1:lf2;
        f2(jr)=mean(f(jr*numfb-numfb+1:jr*numfb));
      end; 
 
      for jc=1:m;
        for jr=1:lf2;
          [p2(jr,jc)]=mean(p(jr*numfb-numfb+1:jr*numfb,jc));
        end; 
      end;
    else;
      p2=p;
      f2=f;
    end
%
%   Calculate phase confidence intervals (using the "stats toolbox")
%
    degf=floor(nump*numfb-2);
    st=tinv(.975,degf);
%
    [nr,nc]=size(p2);
    for j=1:nr;
      ciph(j)=180/pi*asin(st*sqrt((1-p2(j,5)^2)/(degf*real(p2(j,5))^2)));
    end;
%
%   calculate the coherence C.I.
%
    degf=floor(nump*numfb/2)
    cic=sqrt(1-.05^(1/(degf-1)));
  end;
%
  if skip==2;
    if hmp<3;
      cont=menu('Now what?','view cross-spec plot',...
                        'output results to ascii file',...
                        'both ',...
                        'continue');
    else;
      cont=3;
    end;
  end;


  if cont==1|cont==3;
    spcpltjm2(p2,f2,maxcypd,ciph,cic,99);
%   cspecp=int2str(npt)
%    print cspecp
%      !lpr cspecp.ps
%      title([num2str(icode(wv1)),' vs ',num2str(icode(wv2))]);
      title(['eastward wind vs ',num2str(icode(wv2))]);
  end

  if cont==2|cont==3;
%
%   the columns output are:
%
%   var#1 var#2   freq  days     gain          phase     phase confidence   coher. 95%
%
%
%   Note also that we output only periods greater than one day
%
%   (NOTE: ONE MAY NEED A TRANSPOSE ON f2 in the FOLLOWING LINE)

     p3=[ones(nr,1)*icode(wv1) ones(nr,1)*icode(wv2)  f2*100. 1./f2 abs(p2(:,4)) ...
         180/pi*angle(p2(:,4)) abs(ciph)' p2(:,5) ones(nr,1)*cic];
     idd=find(p3(:,4)>1&p3(:,8)>=p3(1,9));
     eval(['fploopy(p3(idd,:),''spec.out'');']);
%    pr=menu('Print?','y','n');
    if idd~=[];
      orient tall
%      pps([num2str(icode(wv1)),'vs',num2str(icode(wv2))]);
      pps(['tauyvs',num2str(icode(wv2))]);
    end;

  end
  end;
% start another variable 2
end;
% start another variable 1
end;

% end cross-spec phase of foi
end
%
%%%%%%%%%%%%%%%%%%%%%       EOFs           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if todo == 5

   disp('Calculating EOFS.... hold on.')

   % Note we are using the filtered version of the time series, "yf".
   %
   % each row in "yf" should represent one observation
   %
   % remove the mean from yf; let "m" = #observation, "n" = #pts
   %
   [m,n]=size(yf);
   ybar=mean(yf);
   %
   % need to make "ybar" a m*n matrix with each columns of constant values
   % so that we can subtract it from "y". So far "ybar" is only a 1 x n matrix
   %
   ybar=ones(m,1)*ybar;
   ymybar=yf-ybar;
   %
   % Now, do the eof 2-liner!
   %
   c=cov(ymybar);
   [v,d]=eig(c);
   % 
   % write out results of eofs
   %
     trace=sum(d); 
     perctvar=d/trace*100.;
     disp('   mode#     percent variance')
     ld=length(d);
     od=[ld:-1:1];
     result=[od;perctvar'];
     disp(result');
     disp('Press Return to Continue.');
     pause
   %
     plot(v(:,n));
     title(' first mode');
     disp('Press Return to Continue.');
     pause;
     hold off
     plot(v(:,n-1));
     title(' second mode');
     disp('Press Return to Continue.');
     pause;
     hold off
   %
   % calculate eof coefficients (mode time dependence)
   %
   % in Fortran we would write:
   %
   % for l=1:m;
   %   for j=1:n;
   %    for i=1:n;
   %    a = a + ymybar(l,i)*v(i,j);
   %    end
   %   end
   %  end
   %
   % But in MATLAB, we simply write:
   %
     a=ymybar*v;
   %
     plot(t,a(:,n));
     title(' first mode time dependence')
     disp('Press Return to Continue.');
     pause;
     plot(t,a(:,n-1));
     title(' second mode time dependence')
     disp('Press Return to Continue.');

     save eofs result t a v;
end
%
%%%%%%%%%%%%%%%%%%%%%   STATISTICS         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if todo==6
 statsout(y,t)
end
%
%%%%%%%%%%%%%%%%%%%%%   KEYBOARD SESSION   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if todo == 99
% 
% Note: Your variables are:
who
%
%  You may enter any MATLAB command at the prompt below. When finished
%  type RETURN.
%
% 
  keyboard
  else
end
%
%%%%%%%%%%%%%%%%%%%%%  VIEW FLOW FIELDS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if todo==8;
  ff=menu('Choose a flow field','Cod season (no density)',...
        'Fall85 wind event',...
        'Georges Bank Rotation');
  axis('off');
  if ff==1;
%
%%%%%%%% CODNODEN
%
  !gzip -d /home3/ocn/jmanning/model/d3/codnoden.mat
  load /home3/ocn/jmanning/model/d3/codnoden
  movie(egg,3,12,[0 0 600 600]);
  !gzip  /home3/ocn/jmanning/model/d3/codnoden.mat
  end;
%
 if ff==2;
%
%%%%%% FALL 1985 WIND EVENT
%
  !gzip -d /home3/ocn/jmanning/gbr/movies/fall85.mat
  load /home3/ocn/jmanning/gbr/movies/fall85
  movie(fall85,2,12,[0 0 600 600]);
  clear fall85
  !gzip    /home3/ocn/jmanning/gbr/movies/fall85.mat
  end;
%
  if ff==3;
%
%%%%%  ROTATE GEORGES
%
%  !gzip -d /home3/ocn/jmanning/gbr/movies/rotbank.mat
%  load /home3/ocn/jmanning/gbr/movies/rotbank
%  movie(rotbank,3,12,[0 0 600 600]);
%  !gzip  /home3/ocn/jmanning/gbr/movies/rotbank.mat
  end;
end;
%
%
%%%%%%%%%%%%%%%%%%%%   PROCESS CTD DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if todo==9;
  lookat
end;
%
%%%%%%%%%%%%%%%%%%%%   Principal Axis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if todo==10;
  [deg xs ys]=prinaxis
end;
%
%%%%%%%%%%%%%%%%%%%%   RESTART FOI PROGRAM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
foi

