%  routine to read RV OCEANUS shipboard data 
%
%  OLD DATA LOAD PHASE BECAME OBSOLETE in 1998
%  DUE TO LAURA"S WEB SERVED ASCII DATA
%
%  note: this routine calls other routines in this directory

%  note: I have only modified phase 1 (ph=1: loading raw data).
%  Phase 2 (rough plots) and phase 3 (cleaning data)
%  still need work but those functions may be part of "ATSEE.M"
%  routines in the future.

 
%ph=menu('What?','load data','plot dat');
ph=2;

if ph==1;
   
   %  wmth=menu('Which Method?','old: reading raw data','new: reading web served data');
   wmth=2;
   if wmth==1;

  day1=input('Enter day of the month or 99 for range of days:  ');

% generate a new number "dayn" given the day of the month

  if day1==99;
    day1=input('First day?');
    day2=input('Last day?');
  else;
    day2=day1;
  end;

% make this a yearday
  dayn1=day1;
  dayn2=day2;

  disp('');
%  inc=input('How often do you want data in minutes?');
inc=1;

% note that pcnfs says mounted the h drive as 
% "net use h: april:/home/april/pcnfs/athena"
% However, since we are working directly on ships computer
% or assuming the data is on the lab computer we hardcode:
  tree='/data5/jmanning/ship/oc303/oc9705';
  disp('CURRENTLY HARDCODED to FILL HEADING MISSING WITH 999.9');

% loop through all the days
  for j=dayn1:dayn2;

    if j<10;zer='0';else;zer='';end;

% read all the files one at a time
    filena=[tree,zer,num2str(j),'.001'];
    fid=fopen(filena,'r');

% allow a second chance the file might have a .002 extension
    if fid==-1;
      filena=[tree,zer,num2str(j),'.002'];
      fid=fopen(filena,'r');
    end;

    if fid==-1;
      disp(['Trouble finding file: ',filena]);
      return;
    else;   
      disp(' ')
      disp(['Loading ',filena, ' every ',num2str(inc),' minutes...']);
      disp(' ')
    end;

% skip by header
    for kl=1:16;fgetl(fid);end;

% keep reading until we pass the blank line

% read each line and extra info according to the leading three digit code
    k=0;
    alw=0;
    lin=fgetl(fid);
    while lin ~=''
      lin=fgetl(fid);
    end;
    lin=fgetl(fid);
    m=0;

    while lin~=-1;
      if lin~='';
      
%       look for day,mth,yr
        if lin(1:3)=='002';
          k=k+1;
          m=m+1;
          mth=str2num(lin(6:7));
          day(k)=str2num(lin(9:10));
          yr=str2num(lin(12:13));
        end;

%       read time (GMT)
        if lin(1:3)=='003';
          m=m+1;
          hr(k)=str2num(lin(6:7));      
          mn(k)=str2num(lin(9:10));
          sc(k)=str2num(lin(12:13));
          yd(k)=julian(mth,day(k),yr,hr(k)*100+mn(k));
%         write out every hour
          if mn(k)==0;
            disp(['Hour ',num2str(hr(k))]);
          end;
        end;

%       read ships heading from gyro
        if lin(1:3)=='011';
          if lin(6)~='-';
            m=m+1;
            h(k)=str2num(lin(6:length(lin)));
          else;
            h(k)=999;
          end;
        end;
      

%       read IMET data (w/wind relative to ship)      
        if lin(1:3)=='013'
          m=m+1;
          pswind(k)=-1*str2num(lin(6:11));
          idn=find(lin(11:length(lin))==' ');
          fawind(k)=-1*str2num(lin(11+idn(1)+1:11+idn(2)-1));

          id=find(lin==';');
          if length(id)>=2;
%           get baro pressure
            barp(k)=str2num(lin(id(1)+1:id(2)-1));
          else;
            barp(k)=-999;
          end;
          if length(id)>=3;
            ids=find(lin(id(2):id(3))==' ');
            airt(k)=str2num(lin(id(2)+ids(length(ids)):id(3)-1));
          else;
            airt(k)=999;
          end;
          if length(id)>=4;
            swr(k)=str2num(lin(id(3)+1:id(4)-1));
          else;
            swr(k)=-999;
          end;
        end;

%       read lat& lon      
        if lin(1:3)=='029'
          m=m+1;
%         check to make sure we have a fix. If not, substitute with previous
          if length(lin)<34;
            if k>1;
              lat(k)=lat(k-1);
              lon(k)=lon(k-1);
            else;
%             use previous days
              [nrr,ncc]=size(d);
              lat(k)=d(length(d),6);
              lon(k)=d(length(d),5);
            end;
          else;
            lat(k)=str2num(lin(13:21));
            lon(k)=str2num(lin(26:34));
          end;
        end;

%       read seasurface conductivity
        if lin(1:3)=='024';
          m=m+1;
          ll=length(lin);
          cond(k)=str2num(lin(6:ll));
        end;

%       read seasurface salt
        if lin(1:3)=='064';
          if length(lin)>7;
            m=m+1;
            ll=length(lin);
            sal(k)=str2num(lin(6:ll));
          else;
            sal(k)=999;
          end;
        end;

%       read sea temperature
        if lin(1:3)=='025';
          if length(lin)>8;
            m=m+1;
            ll=length(lin);
            sst(k)=str2num(lin(6:ll));
          else;
            sst(k)=999;
          end;
        end;

%       loop testing allowing for blank line
      elseif lin=='';

%       if not all the variables were there redo this point
%       or fill the data as missing value

        if m~=8;
          hh=[length(day) length(hr) length(h) length(fawind) ...
            length(lat) length(cond) length(sal) length(sst)];
          if alw==1;
%           disp('CURRENTLY HARDCODED to FILL HEADING MISSING WITH 999.9');
%            if find(hh==min(hh))==3;
              h(k)=999.9;
%            else;
%              disp('fill the missing point with some value');
%              keyboard
%            end;
              
          else;
            disp('There should be 8 sets of variables');
            disp('day hr head imet lat cond salt sst');
            disp(['Missing the ',num2str(find(hh==min(hh))),'th data']);
%            wtd=menu('What to do?','skip','fill w/999.9 always','fill this time');
            wtd=2;
            if wtd==1;
              k=k-1;
            elseif wtd==2;
              alw=1;
              sal(k)=999.9;
            elseif wtd==3;
              disp('fill the missing point with some value');
              keyboard
            end;
          end;
        end;
        m=0;

      end;


% now read another line and go back to beginning of the "while" loop
      lin=fgetl(fid);

    end;

% convert lat lon to decimal degrees
     [la,lo]=ddmm2dec(lat,lon);

% make longtitude negative
%    lo=lo*(-1);

% TRY correcting  wind for ships motion
%    twind=menu('try correcting wind?','y','n');
    twind=1;
    if twind==1;
     disp('Hardcoded to use ships bearing as heading');
     for i=2:length(pswind);

%     calculate ships speed & heading
      [distkm(i),h(i)]=dist(la(i-1),lo(i-1),la(i),lo(i));
      sspd(i)=distkm(i)*1000/((yd(i)-yd(i-1))*86400);

%     subtract ship speed from fore & aft wind
%     remove ship motion (note: that we add instead of subtract
%     to convert ships wind to oceanographic convention)
      fawind2(i)=fawind(i)+sspd(i);

%     calculate measured e/w wind from port/strd wind & heading
%      [wu(i),wv(i)]=rot(-1*h(i),pswind(i),fawind2(i));
       [ss(i),dd(i)]=uv2sd(pswind(i),fawind2(i));
       dd(i)=dd(i)+h(i);
       if dd(i)>=360;dd(i)=dd(i)-360;end;
       [wu(i),wv(i)]=sd2uv(ss(i),dd(i));
     end;
     wu(1)=wu(2);wv(1)=wv(2);fawind(1)=fawind(2);pswind(1)=pswind(2);
% estimate stress assumes m/s wind and anemomitor at 20m on mast
     disp('calculating stress... please wait');
     fac=1.0;hgt=20;
     for i=1:length(wu);
      [taux(i),tauy(i)]=wstr(wu(i)*fac,wv(i)*fac,hgt);
      rslt(i)=sqrt(taux(i)^2+tauy(i)^2);
     end;
    else;
     wu=ones(length(la),1)*999;
     wv=ones(length(la),1)*999;
     taux=wu;
     tauy=wv;
     rslt=wu;
    end;

% this day is done so add these to final output matrix
    if j==dayn1;
      d=[day; hr; mn; yd; lo; la; sst; sal; cond; airt; wu; wv; barp; swr; taux; tauy]';
      if dayn1~=dayn2; 
        clear day  hr  mn  yd  lon lat lo  la ewind nwind  wu  wv  barp ...
            h pswind fawind sst  cond taux tauy airt sal swr
      end;
    else;
      d1=[day; hr; mn; yd; lo; la; sst; sal; cond; airt; wu; wv; barp; swr; taux; tauy];
      d=[d; d1'];
      clear day  hr  mn  yd  lon lat lo  la ewind nwind  wu  wv  barp ...
            h pswind fawind sst  cond taux tauy airt sal swr
    end;

% get another day of data
  end;

% save this big "d" matrix to a .mat file so that you don't need
% to load read it again

  elseif wmth==2;
     % reading laura's web served data
     load c:\gbs99\ship\ship.web
     s=ship;
     
     %correct year day to yrday0_gmt
     s(:,1)=s(:,1)-1;
     
     %calculate hr mn day and wind stress
     disp('Hold on ... calculating wind stress...');
     [uu,vv]=sd2uv(s(:,13),s(:,14));
     for i=1:length(s);
        dd(i,:)=calday(99,s(i,1));
        [taux(i),tauy(i)]=wstr(uu(i),vv(i),17);
     end;   
     d=[dd(:,3) dd(:,4) dd(:,5) s(:,1) s(:,3) s(:,2) s(:,5) s(:,6) s(:,6) s(:,9) uu' vv' s(:,7) s(:,8)];    
  end;   
  
  save c:\gbs99\ship\ship.dat d -ascii;
%  saved ship d;

% we are done the 1st phase of loading data
end;

  
%%%%%%%%%%%%%%%%%%%%%%%%  Rough PLOT Phase  %%%%%%%%%%%%%%%%%%%%%  

if ph==1|ph==2;

% load the data matrix in "ship.mat" file
%  exi=exist('d');
%  if exi==0;
%    load ship;
%  end;

%  nsp=input('How many subplots for displaying results?');
  nsp=1;
  for ksp=1:nsp;

% call a function that sets subplot command based on how many
  if nsp>1;
    hmsp(nsp,ksp);
  end;

% user needs to edit lines below to specify reasonable ranges
% quick clean
  id=find(d(:,6)>39&d(:,6)<42);
  id2=find(d(id,5)>-69&d(id,5)<-65);
  miyd=input('Enter min yearday: (enter 999 for whole record) ');
%  miyd=999;
  if miyd~=999;
    mayd=input('Enter max yearday: ');
    id3=find(d(id(id2),4)>=miyd&d(id(id2),4)<=mayd);
    dn=d(id(id2(id3)),:);
  else;
    dn=d(id(id2),:);
  end;

wp=menu('What do yoyu want?','Ship Track',...
         'Airtemp & SST time series',...
         'salt time series',...
         'Barom & SWR time series',...
         'SST contour',...
         'wind',...
         'all of the above');

if wp==1|wp==7; 
% plotting ship track
  figure(1);

% put down a base map with bathymetry
% with candela shades and no 80m and no labels
  georbat(1,0,0,2);

% plot the ship track
%  plot(dn(:,5),dn(:,6),'r.','markersize',2);
  plot(dn(:,5),dn(:,6),'r')
% fix up the axis
%  axis([floor(min(dn(:,5))) max(dn(:,5)) min(dn(:,6)) max(dn(:,6))]);
  axis([floor(min(dn(:,5))) max(dn(:,5)) min(dn(:,6)) max(dn(:,6))]);
  [date1]=calday(99,min(dn(:,4)));
  [date2]=calday(99,max(dn(:,4)));
  lllabel(10,30,40,42,-71,-66);

% post values along ship track
  disp('note: dn=[day; hr; mn; yd; lo; la; wu; wv; barp; sst; cond; taux; tauy]');
  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(dn)/hmv):length(dn);
%     post only sgd sig.digit
      idp=find(num2str(dn(i,wval))=='.');
      if idp~=[];
        post=num2str(dn(i,wval));
        if sgd==0;sgd=sgd-1;end;
          text(dn(i,5),dn(i,6),post(1:idp+sgd),'color','b',...
           'horizontalalignment','center');
      else;
        text(dn(i,5),dn(i,6),num2str(dn(i,wval)),'color','b',...
           'horizontalalignment','center');
      end;
    end;
  end;

% title with date
  title(['From ',num2str(date1(:,2)),'/',num2str(date1(:,3)),' to ',...
              num2str(date2(:,2)),'/',num2str(date2(:,3))])  

  
% put this figure in upper right portion of the screen 
  set(gcf,'units','normal','position',[.5 .5 .4 .4])    


end;

if wp==2|wp==7;
% example of a time series plot (seasurface & air temp vs yearday)
  figure(2);
  id3=find(dn(:,10)>0&dn(:,10)<30);
  id4=find(dn(:,7)>0&dn(:,7)<30);
  plot(dn(id3,4),dn(id3,10),'r','linewidth',2);
  hold on;
  plot(dn(id4,4),dn(id4,7),'g');
  title('Raw Hull Mounted Thermistor Records');
  legend('Air','Seasurface',2);
  ylabel('degC')
  yd2calab(99);
  set(gcf,'units','normal','position',[.5 .1 .4 .4])    
  figure(7);
  [ax,pl,ylab]=plotyy(dn(id3,4),dn(id3,10),dn(id4,4),dn(id4,7),'seasurface temp',...
    [0 0 0 min(dn(:,4)) max(dn(:,4)) 0 20 2 7]);
  ylabel('air temp');
  set(pl(1),'color','r');
  set(ax(1),'ycolor','r');
  set(pl(2),'color','g');
  set(ax(2),'ycolor','g');
  set(ylab,'color','g');
  set(gca,'xtick',[76 77 78 79 80 81]);
  yd2calab(99)
  title('Raw Hull Mounted Thermistor Records');
%  print
end;

if wp==3|wp==7;
% example of a time series plot (yd vs salt)
  figure(3);
  id3=find(dn(:,8)>30&dn(:,8)<37);
  plot(dn(id3,4),dn(id3,8),'r.');
  title('Raw SAIL System Salinity Record');
  ylabel('PSU')
  set(gca,'position',[.1 .1 .8 .8])
  axis([min(dn(:,4)) max(dn(:,4)) 31 37.5])
  %set(gca,'xtick',[127 130 133 136]);
  yd2calab(99);
  set(gcf,'units','normal','position',[.1 .1 .4 .4])    
%  print
end;

if wp==4|wp==7;
% example of a time series plot (barom & swr vs yearday)
  figure(4);
  subplot(2,1,1);
  id3=find(dn(:,13)>800&dn(:,13)<1200);
  id4=find(dn(:,14)>0&dn(:,14)<1500);
  plot(dn(id3,4),dn(id3,13),'r','linewidth',2);
  ylabel('Barometric Pressure');
  title('Raw Hull Mounted Sensors');
  set(gca,'xticklabels',' ');
  subplot(2,1,2);
  plot(dn(id4,4),dn(id4,14),'g');
%  legend('Barometric Pressure','Short Wave Radiation');
  ylabel('Short Wave Radiation')
  %set(gca,'xtick',[127 130 133 136 139 142]);
  yd2calab(99);
  set(gcf,'units','normal','position',[.1 .5 .4 .4])    
%  print
end;

if wp==5|wp==7;
% try contouring the sst
  figure(5);
  [ro,milo,malo]=range(dn(:,5));
  [ra,mila,mala]=range(dn(:,6));
  xi=milo:ro/25:malo;
  yi=mila:ra/25:mala;
  z=dn(1:10:length(dn),10);
  disp('Hardcoded to 25 grid pts');
%  zi=griddata(dn(1:10:length(dn),5),...
%     dn(1:10:length(dn),6),...
%     dn(1:10:length(dn),10),...
%     xi,yi');
%  
%  v=[2:1:10];
%  disp('Hardcoded contour intervals: ');
%  v
%  %  econtour(xi,yi,zi,v,'fill','label');
%  cd c:\gbs99\atsee\saga
%  extcontour(xi,yi,zi,v,'label');
%  cd c:\gbs99\atsee
  title('OC    Seasurface Temperture');
  lllabel;
  hold on;
  georbat;

  axis([-71 -65 40 43]);

% now overlay the ship track on this 
  plot(dn(:,5),dn(:,6),'r.','markersize',10)

% now blank out areas not sampled
  disp('Click 5 triangles to remove unwanted contours');
  for kk=1:5;
    [xb,yb]=ginput(3);
    fill(xb,yb,'k');
  end;

  set(gcf,'units','normal','position',[.1 .1 .4 .4])    
%  print
end;

% wind
if wp==6|wp==7;
   figure(6);
   wu=plfilt(dn(:,11));
   wv=plfilt(dn(:,12));
   stkcombo(99,dn(:,4),[wu wv wu wv],1,3,4,...
      ['eastward ';'northward';'         '],15,10)
   %plot(dn(:,4),dn(:,11))
   %plot(dn(:,4),wu,'r','linewidth',1);
   %hold on;
   %plot(dn(:,4),wv,'g','linewidth',2);
   %legend('eastward','northward');
  ylabel('m/s');
end;
%
% fix axis on subplots
  if rem(ksp,2)==0;
    ylabel(' ');
  end;
  if ksp<5;
    xlabel(' ');
  end;
  end;

% end of rough plot phase of this routine
end;
