% Plot6.m  Plot TAPS-6 data
%
% 4/1/98   cfg: revised from p6t.m
% 7/31/98  cfg: revised to use rdtaps6 function on log files
% 8/5/98   cfg: revised to correct absorption using measured temperatures
% 8/12/98  cfg: added depth limit inputs for Sv vrs depth plot
% 1/26/00  cfg: revised depth limit code; improved plots
% 7/8/03   cfg: added decoding of yymmdd data

clear all

fprintf('\nPLOT6: program to bin and plot TAPS-6 CAST-MODE log files\n\n')

dsf=1.0;                               % scale depths here

f=[265 420 700 1100 1850 3000];        % TAPS-6 frequencies


[filename,pname] = uigetfile('*.*', 'Select TAPS-6 LOG file');
cd(pname);
cmd=['x=rdtaps6(' '''' filename '''' ');'];
fprintf('Reading log file ...')
eval(cmd)
fprintf(' done\n')

% -------------------------------------------------------------------------
% Enter corrections to the internal TAPS cals here
% -------------------------------------------------------------------------

cals=[0 0 0 0 0 0];          % if no delta's, use zeros
calflag=input('Input delta cals (y/n?):','s');
if calflag == 'y'
   for k=1:6
      cmd=['Delta-Cal(' int2str(k) ') = '];
      cals(k)=input(cmd);
   end
end

% -------------------------------------------------------------------------
% Correct absorption values?
% -------------------------------------------------------------------------

anom=[0 0 0 0 0 0];          			% if no absorption corrections, use zeros

abflag=input('Correct absorption for temperature (y/n?)','s');
if abflag=='y'
   tnom=input('    Enter temperature used for cals (15C?):');
   snom=input('    Enter salinity used for cals (34?):');
   snow=input('    Enter salinity to use for corrections: ');
   anom=absorp(f,tnom,snom);   			% compute vector of nominal absorptions
end

% -------------------------------------------------------------------------
% figure out how many casts we have; plot depth history vrs time of day
% -------------------------------------------------------------------------

date=x(:,1);
yr=fix(date/10000);
mon=fix((date-10000*yr)/100);
day=date-10000*yr-100*mon;
if min(yr)<90,yr=yr+2000;
else
	yr=yr+1900;
end
tod=x(:,2);
hr=fix(tod/10000);
minute=fix((tod-10000*hr)/100);
sec=tod-10000*hr-100*minute;
d=datenum(yr,mon,day,hr,minute,sec);    	% decimal time vector

for k=2:length(d);
	if d(k)<d(k-1),d(k)=d(k-1);end
end

for k=1:length(d)-1;
   dd(k)=abs(d(k+1)-d(k));              	% delta-time vector
end
del=10*mean(dd(1:10));
iv=find(dd>del);                     		% locate cast dividers
numcasts=length(iv)+1;

kc=1;istart=1;iend=length(x);
while kc<numcasts
   kc=kc+1;
   istart(kc)=iv(kc-1)+1;
   iend(kc)=length(x);
   iend(kc-1)=iv(kc-1);
end

h1=figure(1);clf
plot(d,-dsf*x(:,3));grid
datetick('x','HH:MM')
title(['Depth history for ' filename]);
ylabel('Depth - m')
xlabel('Time of day')

fprintf('There are %2.0f casts in this file\n',numcasts')
c2p=input('Which cast do you want to plot? :');

% =====================================================================
%                       MAIN LOOP
% ---------------------------------------------------------------------

while c2p>0 & c2p<=numcasts
   ctitle=['Cast #' int2str(c2p) ' of ' int2str(numcasts)];
   binsize=0;
   while binsize<.001
      binsize=input('Enter depth bin size in m:');
   end
   cflag=input('Do up (u), down (d), or both (b) casts:','s');
   nflag=input('Subtract minimums (y/n)? ','s');

   
   %                   Sort out up/down/both casts
   %                           Set up bins
   %                 NB bins from surface to max depth
   %                  Bin center is midpoint of bin
   
   if exist('atemp')
      clear atemp aSv1 aSv2 aSv3 aSv4 aSv5 aSv6 da;
      clear nd nt nd1 nd2 nd3 nd4 nd5 nd6;
      clear depth Sv1 Sv2 Sv3 Sv4 Sv5 Sv6;
      clear timecode daycode t inst1 inst2;
   end
   
   disp('Setting up depth bins')
   dep = -dsf*x(:,3);
   figure(2);clf;
   plot(dep(istart(c2p):iend(c2p)));grid
   mindepth=min(dep(istart(c2p):iend(c2p)));
   imd=find(dep(istart(c2p):iend(c2p))==mindepth);
   imd=imd(1)+istart(c2p)-1;
   nb=ceil(abs(mindepth/binsize-.5));
   adepth = -(1:nb)*binsize+binsize/2;
   
   atemp = zeros(nb,1);        % vector of binned temperatures
   aSv1 = eps*ones(nb,1);      % vector of binned Sv's
   aSv2 = aSv1;aSv3 = aSv1;
   aSv4 = aSv1;aSv5 = aSv1;
   aSv6 = aSv1;
   nd   = eps*ones(nb,1);nt=nd;
   nd1  = nd;nd2  = nd;nd3  = nd;
   nd4  = nd;nd5  = nd;nd6  = nd;
   
   fprintf('%4.0f bins setup\n',nb)
   fprintf('Maximum depth this cast = %5.2f m\n',mindepth);
   
%  Convert input data
   
   ks=istart(c2p);kend=iend(c2p);  % setup data indices
   if cflag=='d'; kend=imd;end     % downcast only
   if cflag=='u'; ks=imd;end       % upcast only
   
   depth=dep(ks:kend);             % depth
   daycode = x(ks:kend,1);         % not used in this code
   timecode = x(ks:kend,2);        % not used in this code
   t  =  x(ks:kend,4);             % water temperature
   inst1=x(ks:kend,5);             % sensor 1 output
   inst2=x(ks:kend,6);             % sensor 2 output
   Sv1 = x(ks:kend,7)  + cals(1);   % 265 KHz
   Sv2 = x(ks:kend,8)  + cals(2);   % 420 KHz
   Sv3 = x(ks:kend,9)  + cals(3);   % 700 KHz
   Sv4 = x(ks:kend,10) + cals(4);   % 1100 KHz
   Sv5 = x(ks:kend,11) + cals(5);   % 1850 KHz
   Sv6 = x(ks:kend,12) + cals(6);   % 3000 KHz
   ndat=length(Sv1);
   
%   Noise subtraction routine
   
   if nflag == 'y'
      msv1=min(Sv1);msv2=min(Sv2);msv3=min(Sv3);
      msv4=min(Sv4);msv5=min(Sv5);msv6=min(Sv6);
      i1=10^(msv1/10);i2=10^(msv2/10);i3=10^(msv3/10);
      i4=10^(msv4/10);i5=10^(msv5/10);i6=10^(msv6/10);
      for k=1:ndat
         if Sv1(k) ~= nan;Sv1(k)=10*log10(10^(Sv1(k)/10)-i1+eps);end
         if Sv2(k) ~= nan;Sv2(k)=10*log10(10^(Sv2(k)/10)-i2+eps);end
         if Sv3(k) ~= nan;Sv3(k)=10*log10(10^(Sv3(k)/10)-i3+eps);end
         if Sv4(k) ~= nan;Sv4(k)=10*log10(10^(Sv4(k)/10)-i4+eps);end
         if Sv5(k) ~= nan;Sv5(k)=10*log10(10^(Sv5(k)/10)-i5+eps);end
         if Sv6(k) ~= nan;Sv6(k)=10*log10(10^(Sv6(k)/10)-i6+eps);end
      end
   end
   
   %      loop over data set, sorting Sv values into depth bins
   
   fprintf('Binning data now ... ')
   for k=1:ndat;
      z=depth(k);     %if abs(z)<binsize;z=-binsize;end
      l = ceil((-z)/binsize-.5);  % index into bin array
      if l>0 & l<= nb
         nd(l) = nd(l) + 1;  % increment # samples in this bin
      
         atemp(l)=atemp(l)+t(k);
         nt(l)=nt(l)+1;

         z=10 .^(Sv1(k)/10);
         aSv1(l) = aSv1(l) + z;
         nd1(l) = nd1(l) +1;

         z=10 .^(Sv2(k)/10);
         aSv2(l) = aSv2(l) + z;
         nd2(l) = nd2(l) +1;

         z=10 .^(Sv3(k)/10);
         aSv3(l) = aSv3(l) + z;
         nd3(l) = nd3(l) +1;

         z=10 .^(Sv4(k)/10);
         aSv4(l) = aSv4(l) + z;
         nd4(l) = nd4(l) +1;

         z=10 .^(Sv5(k)/10);
         aSv5(l) = aSv5(l) + z;
         nd5(l) = nd5(l) +1;

         z=10 .^(Sv6(k)/10);
         aSv6(l) = aSv6(l) + z;
         nd6(l) = nd6(l) +1;
      end
   end
   atemp = atemp ./ nt;
   
   fprintf('  Done\n')
   
   %  adjust absorption for measured temperature, estimated salinity
   
   if abflag == 'y'
      for k=1:nb
         da(k,1:6)=3*(absorp(f,atemp(k),snow)-anom);
      end
   else
      da=zeros(nb,6);
   end
   
   aSv1 = aSv1 ./ nd1; aSv1 = 10 * log10(aSv1)-da(:,1);
   aSv2 = aSv2 ./ nd2; aSv2 = 10 * log10(aSv2)-da(:,2);
   aSv3 = aSv3 ./ nd3; aSv3 = 10 * log10(aSv3)-da(:,3);
   aSv4 = aSv4 ./ nd4; aSv4 = 10 * log10(aSv4)-da(:,4);
   aSv5 = aSv5 ./ nd5; aSv5 = 10 * log10(aSv5)-da(:,5);
   aSv6 = aSv6 ./ nd6; aSv6 = 10 * log10(aSv6)-da(:,6);
   
   iv=find(aSv1==0);     % code value for empty bins
   if length(iv)>0
      for k=1:length(iv) % replace zeroes with NaN's
         aSv1(iv(k))=NaN;aSv2(iv(k))=NaN;aSv3(iv(k))=NaN;
         aSv4(iv(k))=NaN;aSv5(iv(k))=NaN;aSv6(iv(k))=NaN;
      end
   end
   
   % Set depth scale for plots
   
   mnd=0;
   mxd=max(abs(adepth));
   dd=1;
   if mxd>24;dd=5;end
   if mxd>50;dd=10;end
   if mxd>100;dd=20;end
   if mxd>200;dd=50;end
   mxd = -ceil(mxd/dd+.1)*dd;

   % Set min and max temp for plots
   mnt = floor(min(t/2))*2;
   mxt = ceil(max(t/2))*2;
   
   % Set min and max Sv for plots
   iv=find(isfinite(aSv1));
   maxv=max([max(aSv1(iv)) max(aSv2(iv)) max(aSv3(iv)) max(aSv4(iv)) ...
         max(aSv5(iv)) max(aSv6(iv))]);
   maxv=fix(maxv/5)*5;
   minv = maxv-60;
   
   figure(3);clf
   bx=[aSv1 aSv2 aSv3 aSv4 aSv5 aSv6];  % proper format for TAPSINV program
   semilogx(f,bx,'o',f,bx);grid
   axis([100 10^4 minv maxv])
   xlabel('Frequency - KHz')
   ylabel('Sv - dB')
   s=['Volume scattering spectra'];
   title(s)
   
   %    Plot Sv and temperature vrs depth
   
   plotflag='y';
   while plotflag=='y'
      
      figure(4);clf
      
      subplot(121);
      dy=(mxd-mnd)/15;yt1=mnd+dy;yt2=mnd+2*dy;
      plot(t,depth,'-r');
      axis([mnt mxt mxd mnd]);
      xlabel('Temperature - C')
      ylabel('Depth in m')
      title('Temperature');
      grid;
      text(mnt+.25,yt1,filename)
      text(mnt+.25,yt2,ctitle)
      
      subplot(122);
      plot(aSv1,adepth,'k-',aSv2,adepth,'c-',aSv3,adepth,'m-',...
         aSv4,adepth,'y-',aSv5,adepth,'r-',aSv6,adepth,'g-');
      grid
      axis([minv maxv mxd mnd]);
      xlabel('Sv in dB')
      s=['Scattering - ',num2str(binsize),'M bins'];
      title(s)
      ans=input('Change plot axes (y/n?)','s');
      if ans=='y'
         fprintf('Depth range is %4.1f to %4.1f\n',mnd,mxd);
         ds=input('New upper depth = ','s');
         if length(ds)>0;mnd=str2num(ds);end
         ds=input('New lower depth = ','s');
         if length(ds)>0;mxd=str2num(ds);end
         plotflag='y';
		 if mxd>mnd
			mnd=0;
   			mxd=-max(abs(adepth));
		end
		fprintf('Temperature range is %4.1f to %4.1f\n',mnt,mxt);
		dt=input('New lower temp = ','s');
		if length(dt)>0;mnt=str2num(dt);end
		dt=input('New upper temp = ','s');
		if length(dt)>0;mxt=str2num(dt);end
		fprintf('Sv range is %4.1f to %4.1f\n',minv,maxv);
		dt=input('New lower Sv = ','s');
		if length(dt)>0;minv=str2num(dt);end
		dt=input('New upper Sv = ','s');
		if length(dt)>0;maxv=str2num(dt);end
	  else
         plotflag='n';
      end
      
   end
   
   outfile=input('File name for processed data: ','s');
   if length(outfile) > 1
      ss=['save ' outfile ' adepth bx depth t atemp'];
      eval(ss);
   end
   
   c2p=input('Which cast do you want to plot? :');
end
% ---------------------------------------------------------------------
%                          END
% =====================================================================
