% 
%  routine to read shipboard data and plot it
%
ph=menu('what?','load data',...
                'clean data',...
                'plot individual variables',...
                'stack plot all',...
                'buoy winds');

if ph==1;

  day=input('enter day of the month or 99 for range of days:');
  tree='/data5/jmanning/ship/sj9505/';
    let=['a';'b';'c'];
  if day==99;
    day1=input('First day?');
    dayn=input('Last day?');
    disp(' ');
    disp('loading data.. please wait')
    disp(' ');

% make a single file of appended data
    for im=1:3;
      eval(['!cp ',tree,num2str(day1),'_apr_',let(im),'.lis  ',tree,'99_apr_',let(im),'.lis']);
    end;
    for kd=day1+1:dayn;
      for im=1:3;
        eval(['!cat ',tree,num2str(kd),'_apr_',let(im),'.lis >> ',tree,'99_apr_',let(im),'.lis']);
      end;
    end;
  end;

% find out if raw wind is wanted
%   rw=menu('Do you want raw wind calulated?','y','n');      
    rw=2;

% read all the files one at a time
    for im=1:3;
      filena=[tree,num2str(day),'_apr_',let(im),'.lis'];
      fid=fopen(filena,'r');

% find out how many records exist   
      if im==1;
%       read to end of file
        [aa,count]=fscanf(fid,'%s');
        fseek(fid,0,'bof');
      end;




% do a loop and decode every line
      i=0;
      for k=1:count;
        this=fscanf(fid,'%s',1);
        if this~=[]&length(this)~=4;
          if str2num(this)==95;
            i=i+1;
            if im==1;
              yd=fscanf(fid,'%f',1);
              time=fscanf(fid,'%s',1);
              ta(i)=yd+(str2num(time(1:2))+(str2num(time(4:5))/60))/24.;
              if i>1;
                if ta(i)==ta(i-1);
                  disp('Case of time repeat on yd= ')
                  disp(num2str(ta(i)));     
                  ta(i)=ta(i)+.0001;
                end;
              end;
              latd=fscanf(fid,'%f',1);
              latm=fscanf(fid,'%f',1);
              dummy=fscanf(fid,'%s',1);
              lond=fscanf(fid,'%f',1);
              lonm=fscanf(fid,'%f',1);
              [lat(i),lon(i)]=ddmm2dec(latd*100+latm,lond*100+lonm);

              % read raw wind and ship speed to calculate
              if rw==2;
                dummy=fscanf(fid,'%s',7);
              else;
                dummy=fscanf(fid,'%s',1);
                shps=fscanf(fid,'%f',1);
                shpd=fscanf(fid,'%f',1);
            
%               skip over to wind spd & dir 
                suc=fseek(fid,46,'cof');
                if suc<0; disp('fseek failed');end;
                wnds=fscanf(fid,'%f',1);
                wndd=fscanf(fid,'%f',1);

% calculate corrected wind
                if shps~=[]&wnds~=[]&shpd~=[]&wndd~=[];
                  [us,vs]=sd2uv(shps,shpd);
                  [um,vm]=sd2uv(wnds,wndd);
% you would intuitively want to subtract the ship speed BUT, since
% the ship speed is given 180 different "oceanographic" convention, you add
% The caculation gives wind in meteorologic form which is later converted
% prior to stick plotting.

                  uw(i)=um+us;
                  vw(i)=vm+vs;
                else;
                  i=i-1;
                end;
%             skip over to end of line 
                suc=fseek(fid,37,'cof');
                if suc<0; disp('fseek failed');end;
%               dep(i)=fscanf(fid,'%f',1);

              end;
% now read the met data file
            elseif im==2;
              yd=fscanf(fid,'%f',1);
              time=fscanf(fid,'%s',1);
              tb(i)=yd+(str2num(time(1:2))+(str2num(time(4:5))/60))/24.;
              if i>1;
                if tb(i)==tb(i-1);
                  disp(['Case of time repeat on yd= ',num2str(tb(i))]);     
                  tb(i)=tb(i)+.0001;
                end;
              end;
              dummy=fscanf(fid,'%s',6);
              wspd(i)=fscanf(fid,'%f',1);
              wdir(i)=fscanf(fid,'%f',1);
              barop(i)=fscanf(fid,'%f',1);
              byt1=ftell(fid);
              dummy=fscanf(fid,'%f',1);
              airt(i)=fscanf(fid,'%f',1);
              relh(i)=fscanf(fid,'%f',1);
              sun(i)=fscanf(fid,'%f',1);
% fix the negative values to zero
%              if sun(i)<0;
%                 sun(i)=0;
%              end;
              byt2=ftell(fid);
              if byt2-byt1>25|wspd(i)>900;
                 disp(['drop out at day ',num2str(yd),' at ',time])
                 i=i-1;
              end;

% now read the hull mounted sail sensor data
            elseif im==3;
              yd=fscanf(fid,'%f',1);
              time=fscanf(fid,'%s',1);
              tc(i)=yd+(str2num(time(1:2))+(str2num(time(4:5))/60))/24.;
              if i>1;
                if tc(i)==tc(i-1);
                  disp(['Case of time repeat on yd= ',num2str(tc(i))]);     
                  tc(i)=tc(i)+.0001;
                end;
              end;
              dummy=fscanf(fid,'%s',5);
              sst(i)=fscanf(fid,'%f',1);
              dummy=fscanf(fid,'%f',1);
              salt(i)=fscanf(fid,'%f',1);
              nextf=fscanf(fid,'%s',1);
              if length(nextf)==4;
                if nextf=='****';
                  fluoro(i)=NaN;
                else;
                  fluoro(i)=str2num(nextf);
                end;
              else;
                  fluoro(i)=str2num(nextf);
              end;
           end;
        end;
      end;
    end;
  fclose(fid)
  end;

end;
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%  PHASE 2  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
if ph==2;
%
%                        CLEAN THE DATA
%
% first, set the range and deltas for each variable
% where the first two elements are the range and the third is the delta.
% Note that we check the "time" of each variable as well as the
% the variable itself. Note all times have the same range and delta ("trd").
% We use a function called "rngdelchk.m" to do most of the work.

%
  airtrd=[-50 50 1];
  sstrd=[2 25 2];
  saltrd=[30 36 1];
  fluorord=[0.01 50 1];
  wspdrd=[0 100 1];
  wdirrd=[0 360 180];
  uwrd=[-40 40 20];
  vwrd=[-40 40 20];
  latrd=[39 43 .1];
  lonrd=[-71 -64 .1];
  relhrd=[50 100 10];
  baroprd=[950 1200 20];
  sunrd=[-1 1000 1000];
%
% keep track of the minimum # of pts
%
  minpts=1000000;
%
    disp(' ');
    disp('cleaning air temp')
    ctb=rngdelchk(tb,[0 366 mean(abs(diff(tb)))+2*std(diff(tb))],0);
    cairt=rngdelchk(airt,[airtrd(1:2) mean(abs(diff(airt)))+2*std(diff(airt))],0);
    if length(ctb)<minpts; minpts=length(ctb); end;
%
    disp(' ');
    disp('cleaning sea temp')
    ctc=rngdelchk(tc,[0 366 mean(abs(diff(tc)))+2*std(diff(tc))],0);
    csst=rngdelchk(sst,[sstrd(1:2) mean(abs(diff(sst)))+2*std(diff(sst))],0);
    if length(ctc)<minpts; minpts=length(ctc); end;
%
    disp(' ');
    disp('cleaning salinity')
    csalt=rngdelchk(salt,[saltrd(1:2) mean(abs(diff(salt)))+2*std(diff(salt))],0);
%
    disp('cleaning fluorescene')
    cfluoro=rngdelchk(fluoro,[fluorord(1:2) mean(abs(diff(fluoro)))+2*std(diff(fluoro))],0);
%
    disp(' ');
    disp('cleaning wind')
    cwspd=rngdelchk(wspd,[wspdrd(1:2) mean(abs(diff(wspd)))+2*std(diff(wspd))],0);
    cwdir=rngdelchk(wdir,[wdirrd(1:2) mean(abs(diff(wdir)))+2*std(diff(wdir))],0);
%   convert wind to oceanographic convention u & v
    [uw,vw]=sd2uv(wspd,wdir);uw=-1*uw;vw=-1*vw;
    cuw=rngdelchk(uw,[uwrd(1:2) mean(abs(diff(uw)))+2*std(diff(uw))],0);
    cvw=rngdelchk(vw,[vwrd(1:2) mean(abs(diff(vw)))+2*std(diff(vw))],0);
%
    disp('cleaning GPS')
    cta=rngdelchk(ta,[0 366 mean(abs(diff(ta)))+2*std(diff(ta))],0);
    clat=rngdelchk(lat,[latrd(1:2) mean(abs(diff(lat)))+2*std(diff(lat))],0);
    clon=rngdelchk(lon,[lonrd(1:2) mean(abs(diff(lon)))+2*std(diff(lon))],0);
    if length(cta)<minpts; minpts=length(cta); end;
%  end;
%
    disp(' ');
    disp('cleaning relative humidity')
    crelh=rngdelchk(relh,[relhrd(1:2) mean(abs(diff(relh)))+2*std(diff(relh))],0);
%
    disp(' ');
    disp('cleaning barop')
    cbarop=rngdelchk(barop,[baroprd(1:2) mean(abs(diff(barop)))+2*std(diff(barop))],0);
%
    disp(' ');
    disp('cleaning sun')
    csun=rngdelchk(sun,[sunrd(1:2) mean(abs(diff(sun)))+2*std(diff(sun))],0);
%
end; 
%


if ph==3;

% plot
  ch=menu('What?','Ship Track','depth','wind',...
                  'barometer','air temp','rel. humidity','radiation',...
                  'sst','surfsalt','fluoro');
  if ch==1;
    plot(lon,lat,'c');
  elseif ch==3;
    for kk=1:length(cwspd);
      [u(kk),v(kk)]=sd2uv(cwspd(kk),cwdir(kk));
    end;
    plot(fracday,cwspd);
  elseif ch==4;
    plot(fracday,cbarop);
 elseif ch==5;
    plot(fracday,cairt);
  elseif ch==6;
    plot(fracday,crelh);
  elseif ch==7;
    plot(fracday,csun);
  elseif ch==8;
    plot(fracday,csst);
  elseif ch==9;
    plot(fracday,csalt);
  elseif ch==10;
    plot(fracday,cfluoro);
  end;
end;


if ph==4;

% here's where we make a stack plot 

  clf


% choose to limit the time until after the 4th of apr
%  idt=find(ctc>93);
%  t=ctc(idt)';
%  y=[cfluoro(idt); csalt(idt); csst(idt)]';
   t=ctc(1:length(ctc)-1)';
   y=[cfluoro; csalt; csst]';

% in cases where there is one more time element fix it
   if length(y)-length(t)==1;
     y=y(1:length(y)-1,:);
   end;
 
  stkcombo(95,t,y,0,3,2,['fluoro (volt)';'  salt (psu) ';'   sst (degC)'])
  keyboard
  clf

% choose to limit the time until after the 4th of apr

  idt=find(ctb>93);
  t=ctb(idt)';

% First, find wind components

    for kk=1:length(idt);
      [u(idt(kk)),v(idt(kk))]=sd2uv(cwspd(idt(kk)),cwdir(idt(kk)));
    end;
%  y=[csun(idt); crelh(idt); cairt(idt); cbarop(idt); u(idt);  v(idt)]';
  y=[csun(idt); cairt(idt); cbarop(idt); cwspd(idt)]';
  stkcombo(95,t,y,0,4,2,...
 ['   sun (rad) ';'  airt (degC)';' barop (mbar)';'  wind (knts)'])

% save these variables?
  sv=menu('save?','y','n');
  if sv==1;
    save /data5/jmanning/ship/sj9505/sj9505 ...
      clat clon cfluoro csalt csst csun cairt cbarop cwspd cwdir cta ctb ctc;
  end; 
end;

%
% PHASE 5  Plot more wind
if ph==5;
%  rcnv=menu('convert file?','y','n');
  rcnv=2;
  if rcnv==1;
    !/home3/ocn/jmanning/gbs/wind/str2mtl.exe
  end;
  subplot(3,1,1);
  load /data4/jmanning/wind/z44011.mtl
  b=z44011;
  for i=1:length(b);
    byd(i)=julian(b(i,3),b(i,4),b(i,2),b(i,5)*100+b(i,6));
  end;
  ib=find(byd>=130&byd<138);
  by=[b(ib,10) b(ib,11)];
  stkcombo(95,byd(ib),by,1,1,2,['Buoy 44011']);
  set(gca,'position',[.1 .7 .8 .3])

  subplot(3,1,2);
  load /data4/jmanning/wind/z44008.mtl
  b2=z44008;
  for i=1:length(b2);
    byd2(i)=julian(b2(i,3),b2(i,4),b2(i,2),b2(i,5)*100+b2(i,6));
  end;
  ib2=find(byd2>=130&byd<138);
  by2=[b2(ib2,10) b2(ib2,11)];
  % spline fit to same time as 1st data set
  sb2u=spline(byd2(ib2),by2(:,1),byd(ib));
  sb2v=spline(byd2(ib2),by2(:,2),byd(ib));
  stkcombo(95,byd(ib),[sb2u; sb2v]',1,1,2,['Buoy 44008']);
  ylabel('meters/second');
  set(gca,'position',[.1 .4 .8 .3])

  subplot(3,1,3);
    for kk=1:length(cwspd);
      [u(kk),v(kk)]=sd2uv(cwspd(kk),cwdir(kk));
    end;
  by3=[u; v]';
  % spline fit to same time as 1st data set and multiply to get m/s

  sb3u=spline(ctb(find(isnan(by3(:,1))~=1)),by3((find(isnan(by3(:,1))~=1)),1),...
    byd(ib))*(0.51);

  sb3v=spline(ctb(find(isnan(by3(:,2))~=1)),by3((find(isnan(by3(:,2))~=1)),2),...
    byd(ib))*(0.51);
  stkcombo(95,byd(ib),[sb3u; sb3v]',1,1,2,['R/V SEWARD JOHNSON']);
  ylabel(' Wind Speed (meters/second)');
  set(gca,'position',[.1 .1 .8 .3])

% now put them together on single plot
  subplot(1,1,1)
  t=byd(ib);
  y=[by sb2u sb2v];
  stkcombo(95,t,[sb3u; sb3v]'*(0.51),1,1,2,['R/V SEWARD JOHNSON']);
  ylabel(' Wind Speed (meters/second)');
%  set(gca,'position',[.1 .1 .8 .3])

end;
