clear all
close all

uctd1=70200059;
uctd2=70200054;

makeplots=2;                 % 0 for no plots, 1 for most recent, 2 for all
casts=[1:60];                % casts to be loaded 

% identify number of casts
prof=str2num(textread(['data/pne11_nav.asc'],'%3c%*[^\n]','headerlines',1)); 

% load parameters for uctd1 salinity correction (s. uctd_ctd_cmp.m)
load uctd_ctd_cmp.mat

% load navigation file
fid=fopen(['data/pne11_nav.asc'],'r');
n=0;
t=NaN.*ones(size(prof));
lat=t;lon=t;
while feof(fid)==0
      L=fgets(fid);
      n=n+1;
      ofs=0; 
      if n>1 & length(L)>2 
         IS=isspace(L);
         idx=find(IS==1); 
         idx=idx(1);
         m=str2num(L(idx+1+ofs));
         d=str2num(L(idx+3+ofs:idx+4+ofs));
         y=str2num(['20',L(idx+6+ofs:idx+7+ofs)]);
         hh=str2num(L(idx+10+ofs:idx+11+ofs));
         mm=str2num(L(idx+12+ofs:idx+13+ofs));
         t(n-1)=juliand([y m d hh mm 0]);
         ilat=findstr(L,'N');
         lat(n-1)=str2pos(L(ilat-10:ilat));
         ilon=findstr(L,'W');
         lon(n-1)=str2pos([L(ilat+2:ilat+4),' ',L(ilat+5:ilon)]);
         clear IS y m d hh mm ilat ilon
      end
      clear L ofs
end
fclose(fid);

q=NaN.*ones(size(casts));
for ii=1:length(casts)
    q(ii)=find(prof==casts(ii));
end
prof=prof(q);t=t(q);lat=lat(q);lon=lon(q);

% initialize matrices to contain all profiles
dp=1;p=[1:dp:1500]';
sal_raw=NaN.*ones(length(p),length(prof));
sal=sal_raw;temp=sal_raw;

% load data files
ic=0;ie=0;dummy=999;
sn=NaN.*ones(size(prof));cdelay=sn;
for castnum=prof'
	
	filename='pne11_uctd_';
	if (castnum<100) filename=[filename,'0']; end;
	if (castnum<10) filename=[filename,'0']; end;
	filename=[filename,num2str(castnum)];
    ic=ic+1;
	if exist(['data/',filename,'.asc'])  
	   header=char(71,100);data=NaN.*ones(100000,4);
	   fid=fopen(['data/',filename,'.asc'],'r');
	   n=0;nn=0;
	   while feof(fid)==0
             L=fgets(fid);
             if L(1)=='*' | length(L)==2
                n=n+1;
                header(n,1:length(L))=L;
             else
                nn=nn+1;
                data(nn,:)=str2num(L);
             end
             clear L
       end
	   fclose(fid);
   
       % S/N
       if castnum<3
          sn(ic)=str2num(deblank(header(7,15:end)));
       else
          sn(ic)=str2num(deblank(header(9,15:end)));
       end
    
       % determine downcast (sampling freq.: 1/16 1/s)
       if (castnum==43 | castnum==53 | castnum==57 | castnum==58) % "broken casts"
          P=dummy;T=P;C=P;S=P; 
       elseif (castnum==46 | castnum==47 | castnum==60) % calibration casts on CTD
          k=find(data(:,4)==max(data(:,4)));
          k=k(1);
          [kk,pp]=fitsoak(data(:,4),castnum,'verb');
          P=data(kk:k,4);
	      T=data(kk:k,3);
          C=data(kk:k,2);
          scan=data(kk:k,1);
          clear k kk
          % calculate conductivity delay:
          % best C delay is delay at which high frequency variability of 
          % T and C has highest correlation
          [b,a]=butter(3,0.01);
          Cff=filtfilt(b,a,C);
          Tff=filtfilt(b,a,T);
          Cd=C-Cff;
          Td=T-Tff;
          [cc,la]=xcorr(Cd,Td,100,'coeff');
          % example
          if castnum==60
             ie=ie+1; 
             cc_ex(ie,:)=cc;
             la_ex(ie,:)=la;
          end
          cmax=find(cc==max(cc));
          cdelay(ic)=la(cmax);
          clear b a Cff Tff Cd Td cc la
          % offset of 1 entry for best (least spiky on downcast) salinity
          ind1=1:length(P)-1;
          ind2=2:length(P);
	      S=gsw_SP_from_C(C(ind1).*10,T(ind2),P(ind2));
          S=[S(1);S];
          clear ind*
       else
          k=find(data(:,4)==max(data(:,4)));
          k=k(1);K=round(.5.*k);
          SC=data(1:16:K,1);
          p2=gradient(gradient(data(1:16:K,4)));
          kk=find(p2==max(p2));
	      P=data(SC(kk+1):k,4);
	      T=data(SC(kk+1):k,3);
          C=data(SC(kk+1):k,2);
          scan=data(SC(kk+1):k,1);
          clear SC p2 k kk K
          % calculate conductivity delay:
          % best C delay is delay at which high frequency variability of 
          % T and C has highest correlation
          [b,a]=butter(3,0.01);
          Cff=filtfilt(b,a,C);
          Tff=filtfilt(b,a,T);
          Cd=C-Cff;
          Td=T-Tff;
          [cc,la]=xcorr(Cd,Td,100,'coeff');
          % examples
          if (castnum==30 | castnum==35 | castnum==40 | castnum==44)
             ie=ie+1; 
             cc_ex(ie,:)=cc;
             la_ex(ie,:)=la;
          end
          cmax=find(cc==max(cc));
          cdelay(ic)=la(cmax);
          clear b a Cff Tff Cd Td cc la
          % offset of 1 entry for best (least spiky on downcast) salinity
          ind1=1:length(P)-1;
          ind2=2:length(P);
	      S=gsw_SP_from_C(C(ind1).*10,T(ind2),P(ind2));
          S=[S(1);S];
          clear ind*
       end
    
       % interpolate to grid: use downcast only
       if length(P)>1
          for ii=1:length(p)
              q=find(P>p(ii)-dp./2 & P<=p(ii)+dp./2);
              if ~isempty(q)
                 sal_raw(ii,ic)=nmedian(S(q));
                 temp(ii,ic)=nmedian(T(q));
                 clear q
              end
          end 
       end
       
       % salinity correction 
       if sn(castnum)==uctd1
          sal(:,ic)=sal_raw(:,ic)+(a1+b1.*temp(:,ic));
       else
          sal(:,ic)=sal_raw(:,ic); 
       end
       
       if ((castnum==length(lat) & makeplots==1) | makeplots==2)
        
          figure(1)
	      clf
	      orient portrait
		
	      subplot(121);
	      hold on;
	      plot(T,P,'b-','linewidth',1);
          if ceil(max(P)./1e2).*1e2<400
             dy=50;
          elseif ceil(max(P)./1e2).*1e2>=400 & ceil(max(P)./1e2).*1e2<800 
             dy=100;
          else
              dy=200;
          end
	      set(gca,'box','on','xlim',[floor(min(T))-1 ceil(max(T))+1], ...
		          'xtick',[floor(min(T))-1:2:ceil(max(T))+1],...
                  'ylim',[0 ceil(max(P)./1e2).*1e2],'ytick',[0:dy:ceil(max(P)./1e2).*1e2],'ydir','reverse');
	      grid on;
	      xlabel('T [{\circ}C]');
	      ylabel('P [dbar]');
	      ttl=['PNE11: uCTD cast ',filename(length(filename)-2:length(filename))];
	      title(ttl,'fontweight','bold');
		
	      subplot(122);
	      hold on;
	      plot(S,P,'b-','linewidth',1);
	      set(gca,'box','on','xlim',[floor(min(S.*10))./10-.1 ceil(max(S.*10))./10+.1], ...
		          'xtick',[floor(min(S.*10))./10-.1:.3:ceil(max(S.*10))./10+.1],...
                  'ylim',[0 ceil(max(P)./1e2).*1e2],'ytick',[0:dy:ceil(max(P)./1e2).*1e2],'ydir','reverse');
	      grid on;
	      xlabel('Salinity');
	      ylabel('P [dbar]');
	      title(ttl,'fontweight','bold');
        
          if (makeplots==2 & castnum<length(lat))
             disp('paused')
             pause(3);
          end;

       end
    end
    clear filename ttl header data P T C S
end;    % loop over castnum

save uctd_data.mat sn t lat lon p sal_raw sal temp -v6

figure(1)
clf
orient portrait

subplot(311);
hold on;
plot(prof(sn==uctd1),cdelay(sn==uctd1),'bx','linewidth',1);
plot(prof(sn==uctd2),cdelay(sn==uctd2),'rx','linewidth',1);
legend(['S/N ',num2str(uctd1)],['S/N ',num2str(uctd2)],'location','north');
plot(prof(prof==30),cdelay(prof==30),'bo','linewidth',1);
plot(prof(prof==35),cdelay(prof==35),'bo','linewidth',1);
plot(prof(prof==60),cdelay(prof==60),'bo','linewidth',1);
plot(prof(prof==40),cdelay(prof==40),'ro','linewidth',1);
plot(prof(prof==44),cdelay(prof==44),'ro','linewidth',1);
set(gca,'xlim',[prof(1)-1 prof(end)+1],'xtick',[prof(1):3:prof(end)],...
        'ylim',[-2 6],'ytick',[-2:2:6],'box','on');
grid on;    
xlabel('Cast');
ylabel('C-delay [1/16 s^{-1}]');
title('Cross Correlation','fontweight','bold');

subplot(312);
hold on;
plot(la_ex(1,:),cc_ex(1,:),'b-','linewidth',1);
plot(la_ex(2,:),cc_ex(2,:),'r-','linewidth',1);
plot(la_ex(5,:),cc_ex(5,:),'g-','linewidth',1);
legend('Cast 30','Cast 35','Cast 60','location','south');
plot([cdelay(prof==30) cdelay(prof==30)],[-1 1],'b--','linewidth',1);
plot([cdelay(prof==35) cdelay(prof==35)],[-1 1],'r--','linewidth',1);
plot([cdelay(prof==60) cdelay(prof==60)],[-1 1],'g--','linewidth',1);
set(gca,'xlim',[min(la_ex(:)) max(la_ex(:))],'xtick',[min(la_ex(:)):20:max(la_ex(:))],...
        'ylim',[-1 1],'ytick',[-1:1:1],'box','on');
grid on;
xlabel('Scan [1/16 s^{-1}]');
ylabel('Correlation');

subplot(313);
hold on;
plot(la_ex(3,:),cc_ex(3,:),'b-','linewidth',1);
plot(la_ex(4,:),cc_ex(4,:),'r-','linewidth',1);
legend('Cast 40','Cast 44','location','south');
plot([cdelay(prof==40) cdelay(prof==40)],[-1 1],'b--','linewidth',1);
plot([cdelay(prof==44) cdelay(prof==44)],[-1 1],'r--','linewidth',1);
set(gca,'xlim',[min(la_ex(:)) max(la_ex(:))],'xtick',[min(la_ex(:)):20:max(la_ex(:))],...
        'ylim',[-1 1],'ytick',[-1:1:1],'box','on');
grid on;
xlabel('Scan [1/16 s^{-1}]');
ylabel('Correlation');

print -depsc uctd_cdelay.eps
