% "procadcp.m" ... a routine for processing adcp files AFTER they have been
%  run through "bb2mat" routine
%
%  Five phases: 1) loading data & run "toabs.m"
%               2) subsample & plot ubar & vbar time series
%               3) Merging with ship's lat&lon & save yadcp
%               4) detide (subroutine tideadcp.m if needed)
%               5) contouring cross-section

% Jim Manning Summer '97

% specify a year
yr=97;

% specify a phase of processing
ph=menu('Which Phase of Processing?',...
 ' 1) loading data (combining bb2mat output files) and run "toabs.m"',...
 ' 2) subsample & plot ubar & vbat time series',...
 ' 3) Merging with ships lat&lon and prep for versec.m',...
 ' 4) detide (subroutine tideadcp.m if needed)',...
 ' 5) contouring cross-section');


if ph==1;
%  get the Inputfile prefix from user
%  inf=input('Inputfile prefix (including directory? ','s');
%  inf='/data5/jmanning/adcp/oc303/oc303_6';
%   inf='/data5/jmanning/adcp/al9805/000';
  inf='/data5/jmanning/adcp/al9707/f27';

% loop through set of files until no more are found (mf=0)
  exf=1;
  ii=0;
  while exf~=0;
    ii=ii+1;
    disp(['file #',num2str(ii)]);
    eval(['load ',inf,num2str(ii)]);

%   run the conversion to scientific units on this file
    toabs;
 
%   save only variables of interest from this file
    if ii==1;
      uat=ua;vat=va;wat=wa;ubart=ubar;vbart=vbar;julian_dayt=julian_day;
      bottom_deptht=bottom_depth;mean_latitudet=mean_latitude;mean_longitudet=mean_longitude;
      tempt=trans_temp;
    else;
      uat=[uat ua];vat=[vat va];wat=[wat wa];ubart=[ubart ubar];vbart=[vbart vbar];
      julian_dayt=[julian_dayt julian_day];
      bottom_deptht=[bottom_deptht bottom_depth];
      mean_latitudet=[mean_latitudet mean_latitude];
      mean_longitudet=[mean_longitudet mean_longitude];
      tempt=[tempt trans_temp];
    end;
    eval(['exf=exist(''',inf,num2str(ii+1),'.mat'');']);
  end;

% Here's where we might add another .mat file if one is available
  an=menu('Add another .mat file?','y','n');
  if an==1;
%   but first we save the variables just generated
    disp('Saving the present arrays ...');
    julian_dayt2=julian_dayt; ubart2=ubart; vbart2=vbart; uat2=uat; vat2=vat;bottom_deptht2=bottom_deptht;
    prev=input('Enter previous mat file name: ','s');
    eval(['load ',prev]);
%   now put them together & assume we just did the first
    disp('Builing the combined arrays ...')
    julian_dayt=[julian_dayt julian_dayt2];
    bottom_deptht=[bottom_deptht bottom_deptht2];
    clear julian_dayt2 
    ubart=[ubart ubart2];
    clear ubart2 
    vbart=[vbart vbart2];
    clear vbart2 
    uat=[uat uat2];
    clear uat2
    vat=[vat vat2];
    clear vat2
    wat=[wat wat2];
    clear wat2
  end;

  disp('Saving to unsubsampled inf file....')
  eval(['save ',inf,'  julian_dayt ubart vbart uat vat wat depth bottom_deptht mean_latitudet mean_longitudet tempt']);

end;

%%%%%%%%%%%%%%%%%%%%  PHASE 2 subsample & plot u & v time series %%%%%%

if ph==2;

% load the data from phase 1 if it isn't already in workspace
  exj=exist('julian_dayt');
  if exj==0;
    disp('Loading unsubsampled adcp data...');
%    inf='/data5/jmanning/adcp/oc303/oc303_6';
     inf='/data5/jmanning/adcp/al9805/phase1';
    eval(['load ',inf]);
  end;

% subsample vector to make basic plots
% ssr=input('Enter subsample rate: ');
  ssr=1;
  disp(['Subsample rate hardcoded to ',num2str(ssr)]);
%  stop=length(ubart);
  stop=length(wat);
  ua=uat(:,1:ssr:stop);
  va=vat(:,1:ssr:stop);
  wa=wat(:,1:ssr:stop);
  ubar=ubart(1:ssr:stop);
  vbar=vbart(1:ssr:stop);
  bottom_depth=bottom_deptht(1:ssr:stop);
  lat=mean_latitudet(1:ssr:stop);
  lon=mean_longitudet(1:ssr:stop);
  disp('Hardcoded to not remove one day from the adcp yearday');
  jd=julian_dayt(1:ssr:stop);
  temp=tempt(1:ssr:stop);

  % add drifter (or mooring) time series to compare with adcp
  addd=menu('Add other measure of velocity?',...
       'drifter','model','drift & model','drift/mod/candela/moody');
  if addd>2;load /data5/jmanning/model/bank150/may98/r71/ts.dat;end;
  if addd>=4;
    % add Candela tide
    leap=floor(abs(92-yr)/4);
    if yr>=98;sb1=0;disp('Candela tide NOT corrected for >=98');else;sb1=0;end; % corrects phase problem
    [utide,vtide]=detide(lon,lat,jd*24+(yr-94)*(365*24)+leap*24-sb1);
  end;

% plot raw vertically averaged current
%  subplot(2,1,1);
    plot(jd,ubar,'r');
    hold on;
    if addd~=2;sumplt;end;
    if addd>2;plot(ts(:,6),ts(:,10),'g','linewidth',2);
    id=find(ts(:,12)>-17&ts(:,12)<14);plot(ts(id,6),ts(id,13),'g--');end;
    if addd>=4;plot(jd,utide,'b','linewidth',2);moody(1);end;

%    title(inf);
    title('May 1998 Observed & Modeled Flow')
    ylabel('E/W vertically average flow (m/s)');
%    set(gca,'xticklabels',' ');
%    axis([min(jd) max(jd) -1 1]);
    axis([125 126.8 -.8 .8]);
    legend('ADCP','Drifter@15m','Drifter@33m','QUODDY','QUODDY@16m','Candela','Moody@15m',-1);
    yd2calab(98)
  subplot(2,1,2);

    hold on;
    plot(jd,vbar,'r');
    if addd~=2;sumplt;end;
    if addd>2;plot(ts(:,6),ts(:,11),'g');end;
    if addd>=4;plot(jd,vtide,'m');moody(2);end;
%    yd2calab(yr)
    ylabel('N/S vertically average flow (m/s)');
%    axis([min(jd) max(jd) -1 1]);
    axis([125.5 127 -.7 .7]); 
    
   legend('A','D','Q','C','M',-1);

%  subplot(3,1,3);
%    georbat;
%    plot(lon,lat)
%    %  axis([min(lon) max(lon) min(lat) max(lat)]);
%    axis([-68 -67.5 40.5 41.25])
%    lllabel;
%    pause(2);


% save all this in a .mat file
%  disp('Saving to subsampled inf_t file....')
  eval(['save ',inf,'t  julian_dayt tempt ubart vbart uat vat wat bottom_deptht jd ubar vbar ua va wa bottom_depth depth lat lon mean_latitudet mean_longitudet']);

end;

%%%%%%%%%%%%%%%%    PHASE 3 Plotting ship track and (if necessary)
%%%%%%%%%%%%%%%%    GETTING LAT/LON FROM SHIP %%%%%%%

if ph==3;

% load the data from phase 2 if it isn't already in workspace
  exj=exist('jd');
  if exj==0;
    disp('Loading processed adcp data..');
%    load /data5/jmanning/adcp/oc303/oc303_6t
%    load /data5/jmanning/adcp/al9805/phase1t
    load /data5/jmanning/adcp/al9707/f27t

  end;

  merge=menu('Do you need to get lat/lon from ship data?','y','n');
  if merge==1;

    % merge with ship lat/lon
    disp('Loading ship data...');
    %  wop=menu('which?','processed ship','unprocessed');
    wop=1;

    % load ship data
    if wop==1;
      load /data5/jmanning/ship/al9805/ship.dat
      d=ship;
      los=d(:,3);
      las=d(:,2);
      disp('');
      disp('WARNING: I AM ADDING 4 HOURS TO SHIP TIME TO GET GMT');
      disp('         AND SUBTRACTING ONE DAY')
      disp('');
      yds=d(:,1)+4/24-1;
    else;
      load /data5/jmanning/ship/al9707/gps1_178.mtl
      d=gps1_178;
      los=d(:,7);
      las=d(:,6);
      disp('');
      disp('WARNING: I AM ADDING 4 HOURS TO SHIP TIME TO GET GMT');
      disp('         AND SUBTRACTING ONE DAY')
      disp('');
      yds=d(:,3)+d(:,4)+4/24-1;
    end;

    % find min& max of shiptime
    mist=min(yds);mast=max(yds);

    % find the adcp times that fit within this
    % extra .001 added to insure time fits
    ida=find(jd>=mist+.001&jd<=mast-.001);
    jd=jd(ida);ubart=ubart(ida);vbart=vbart(ida);uat=uat(:,ida);vat=vat(:,ida);
    bottom_deptht=bottom_deptht(ida);

    % make sure there isn't any repeat times by not taking ALL ship records
    % because there are often roundoff errors causing negative time difference
    yds=yds(1:5:length(yds));
    los=los(1:5:length(los));
    las=las(1:5:length(las));

    % interpolate lat/lon onto the adcp's julian_day vector
    [la]=interp1(yds,las,jd);
    [lo]=interp1(yds,los,jd);

  else;
    la=mean_latitudet;
    lo=mean_longitudet;
  end;

% plot adcp ship track 
 clf;
 georbat;
 lllabel(10,10,40.5,41.5,-68,-67.5)
 axis([-68 -67.5 40.5 41.5])
  jd1l=julian_dayt-4/24+1;
  jd0g=julian_dayt;
  mint=input('Enter minumum yrday0_gmt? ');
  maxt=input('Enter maxumum yrday0_gmt? ');
  id=find(jd0g>=mint&jd0g<=maxt);
  plot(lo(id),la(id),'r','linewidth',3);
% post values along ship track
  dn=[jd1l;jd1l;jd1l;jd0g;lo;la]';
  [date1]=calday(yr,min(dn(:,1)));
  [date2]=calday(yr,max(dn(:,1)));
  disp('note: dn=[jd1l;jd1l;jd1l;jd0g; lo; la]');
  wval=input('Which one of these do you want? Enter 0 for none of these: ');
  if wval~=0;
    hmv=input('How many pts do you want labeled?');
    sgd=input('How many signif digits do you want labeled?');
    for i=1:floor(length(id)/hmv):length(id);
%     post only sgd sig.digit
      idp=find(num2str(dn(id(i),wval))=='.');
      if idp~=[];
        post=num2str(dn(id(i),wval),7);
        if sgd==0;sgd=sgd-1;end;
        [nr,nc]=size(post);
        if nc-idp<sgd;sgdn=nc-idp;else;sgdn=sgd;end;
        text(dn(id(i),5),dn(id(i),6),post(1:idp+sgdn),'color','b',...
           'horizontalalignment','center','fontsize',6);
      else;
        text(dn(id(i),5),dn(id(i),6),num2str(dn(id(i),wval)),'color','b',...
           'horizontalalignment','center','fontsize',6);
      end;
    end;
  end;
  pause(1);

    disp('plotting arrows from subsampled data');
    hmv=input('How many pts do you want ubar arrows?');
    skip=floor(length(id)/hmv);
    rod=menu('What?','raw current','detided w/Candela');
    if rod==1;
      vecplot2(lo(id(1:skip:length(id))),la(id(1:skip:length(id))),...
         ubart(id(1:skip:length(id))),vbart(id(1:skip:length(id))),...
            1,'m/s',-67.9,41,'g')
    else;
      leap=floor(abs(92-yr)/4);       %leap years
      [utide,vtide]=detide(lo(id(1:skip:length(id))),...
           la(id(1:skip:length(id))),...
           jd0g(id(1:skip:length(id)))*24+(yr-94)*365*24+(leap*24)-leap);
           % note the final "leap" is added to correct Julio's model by 1hr
      uf=ubart(id(1:skip:length(id)))-utide';
      vf=vbart(id(1:skip:length(id)))-vtide';
      vecplot2(lo(id(1:skip:length(id))),la(id(1:skip:length(id))),...
          uf,vf,...
          .2,'m/s');
%          1,'m/s',-67.9,41,'g')
    end;



% save all this in a .mat file that can be read by "getctd" or versec
  disp('Saving to subsampled lat/lon merged  file....')
  disp('For specific times of sections');

% first, make all variables one long vector
  lon=reshape([lo'*ones(1,length(depth))]',length(depth)*length(lo'),1);
  disp('lon reshaped');
  lan=reshape([la'*ones(1,length(depth))]',length(depth)*length(la'),1);
  disp('lat reshaped');
  bdn=reshape([bottom_deptht'*ones(1,length(depth))]',length(depth)*length(lo),1);
  uan=reshape(uat,length(depth)*length(lo),1);
  van=reshape(vat,length(depth)*length(lo),1);
  wan=reshape(wat,length(depth)*length(lo),1);
  disp('ua, va, and wa  reshaped');
  jdn=reshape([jd0g'*ones(1,length(depth))]',length(depth)*length(lo),1);
  depthn=reshape(depth'*ones(1,length(jd0g)),length(depth)*length(lo),1);

%  tt=[124.3 125.0; 125.1 125.15; 125.36 125.6; 125.89 126.166; 126.43 126.542; 126.76 126.854 ];
  k=0;
  tt=[177 179];
%  for i=[1 3 4 5 6 7];
  for i=1;
    k=k+1;
    disp(['Finding deployment #',num2str(i)]);
    id=find(jdn>=tt(k,1)&jdn<=tt(k,2));
    disp(['Saving deployment #',num2str(i)]);
     yadcp=[lon(id) lan(id)   depthn(id) bdn(id) uan(id) van(id) uan(id) van(id) jdn(id)  wan(id)];
    eval(['save /data5/jmanning/adcp/al9707/adcp',num2str(i),'tll yadcp']);
    disp('done');
  end;    

end;

%%%%%%%%%%%%%%%%%%%%%   PHASE 4 DETIDE & Vector plots %%%%%%%%%%

if ph==4;
  tideadcp;
end;


%%%%%%%%%%%%%%%%%  PHASE 5 CROSS _SECTION
if ph==5;
  clf
% assumes phase 1 & 2 has been run and jd1l is in yrday1_local
  jd1l=julian_dayt-4/24+1;
  miny=input('min yrday1_local?');
  maxy=input('max yrday1_local?');
  id=find(jd1l>=miny&jd1l<=maxy);
  uan=uat(:,id)*100;van=vat(:,id)*100;
  
  x=distcalc(mean_latitudet(id),mean_longitudet(id));
  [nru,ncu]=size(uan);
  for m=1:length(x);
    xx(:,m)=x(m)*ones(length(depth),1);
  end;
  for m=1:length(depth);
    yy(m,:)=depth(m)*ones(1,length(x));
  end;
  wvar=menu('Which?','u','v');
%  wvar=2;
  if wvar==1;
    zz=uan;lab='e/w velocity (cm/s)';
  else;
    zz=van;lab='n/s velocity (cm/s)';
  end;
  [rang,minz,maxz]=rangejm(zz);
  v=[floor(minz):floor(rang/10):ceil(maxz)];

% contour direct or grid & contour?
  dog=menu('How to contour?','direct','grid and then contour');
  if dog==1;  
    econtour(xx,yy,zz,v,'fill');
  elseif dog==2;

%   rid of NaN's
    xx1=reshape(xx,nru*ncu,1);
    yy1=reshape(yy,nru*ncu,1);
    zz1=reshape(zz,nru*ncu,1);    
    idnn=find(isnan(zz)~=1);
    [zi,xi,yi]=gridopts(xx1(idnn),yy1(idnn),zz1(idnn),4,.5,[.05 10]);
    econtour(xi,yi,zi,v,'fill');
%     econtour(xi,yi,zi,v,'label');
  end;
  cbar(gca,'v',floor(rang/10),floor(rang/10)*2);
  my=get(gca,'ylim');
  text(min(x),my(1)+abs(my(1))/10,lab,...
    'fontsize',18,'horizontalalignment','left')
  ylabel('Depth (m)');
  xlabel('Alongtrack Distance (km)')
%  title('May 1997 Bongo Transect 3 w/tide','fontsize',14)
%  title('June 1997 Transect 5')

% add blanking of bottom_depth in yadcp(:,4)
  
  axis([0 max(xx1) my(1) my(2)]);
  

end;


