function [] = get_align(crsid,verb)
% function [] = get_align(crsid,verb)
%
% determine the delays of conductivity and oxygen sensors
%
% input  :  crsid        - cruise identifier to compose file names
%           verb     ''  - set to 'verbose' to get more info
%
% version 0.7   last change 06.05.2009

% G.Krahmann, IFM-GEOMAR, Mar 2007
% A.Funk, IFM-GEOMAR, March 2008

% bug in handling dual sensors          GK, Oct 2007    0.1-->0.2
% changed setup of input/output paths   AF, Mar 2008    0.2-->0.3
% moved settings & saving into file     GK, 25.06.2008  0.3-->0.4
% verbose and other changes             GK, 03.07.2008  0.4-->0.5
% new settings scheme                   GK, 03.01.2009  0.5-->0.6
% bug in cases of more than 1 setting   GK, 06.05.2009  0.6-->0.7


%
% load settings
%
eval([crsid,'_settings'])


%
% load start scans of down casts
%
load(['mat/',crsid,'_startscan'])


%
% depthrange      - range over which the data shall be correlated
%
depthrange = allsettings.aligndepthrange;


%
% signs           - correlation or anticorrelation
%                   e.g. [1,-1]  correlation for t-c and anticorrelation for t-o
signs = allsettings.alignsigns;


%
% load previously calculated sensor delays
%
if exist(['mat/',crsid,'_delays.mat'])
  load(['mat/',crsid,'_delays.mat'])
else
  delays.pri_c_ind = nan*ones(1,999);
  delays.sec_c_ind = nan*ones(1,999);
  delays.pri_o_ind = nan*ones(1,999);
  delays.sec_o_ind = nan*ones(1,999);
  delays.pri.c = nan*ones(1,999);
  delays.sec.c = nan*ones(1,999);
  delays.pri.o = nan*ones(1,999);
  delays.sec.o = nan*ones(1,999);
end


%
% loop over primary and seconday sensors
%
for prisec = 1:2

  if prisec==1
    settings = pri_settings;
    str = 'primary';
  else
    settings = sec_settings;
    str = 'secondary';
  end


  %
  % loop over the settings
  %
  for se=1:length(settings)


    disp(' ')
    disp(['deriving delays for ',str,' setting #',int2str(se)])
  

    %
    % get all station numbers for that setup
    %
    prof = settings(se).prof;


    %
    % set range of delays to test
    %
    cond_interp_res = 10;
    shift_range = 200;
    cdelay = [-0.1:0.005:0.1]/86400;
    odelay = [-10:0.5:10]/86400;


    %
    % loop over profiles within this setting
    %
    for n=prof

      %
      % figure out whether this has already been calculated
      %
      do_calc = 1;
      if prisec==1
        if ~isnan(delays.pri.c(n))
          do_calc = 0;
        end
      else
        if ~isnan(delays.sec.c(n))
          do_calc = 0;
        end
      end


      %
      % if this is a new one, do the calculation
      %
      if do_calc==1


        % 
        % load and prepare CTD data
        %  
        [hdr,data] = read_sbe_cnv2mat(['mat/',crsid,'_',int2str0(n,3),...
           allsettings.addname,'.mat'],start.scan(n),'d');
        fn = fieldnames(data);
        p = data.p;
        good = find( (p>depthrange(1) & p<depthrange(2)) );
  

        %
        % do the calculation only if there is sufficient data in the
        % depth range
        %
        if length(good)<1000

          if strcmp(verb,'verbose')
            disp(['no delay calculation for prof #',int2str(n),' because of too few points'])
          else
            fprintf(1,'x',[])
          end

        else

          if strcmp(verb,'verbose') 
            disp(['analyzing profile ',int2str(n),' ',str,' sensors']) 
          else
            fprintf(1,'.',[])
          end

   
          good = [good(1):good(end)];
          p = p(good);
          if prisec==1
            t = data.t_pri(good);
            c = data.c_pri(good);
            o = data.o_pri(good);
          else
            t = data.t_sec(good);
            c = data.c_sec(good);
            o = data.o_sec(good);
          end
          tim = data.timej(good);
  
          %
          % for obscure reasons there occur time jumps
          %
          dtim = diff(tim);
          if any(dtim>0.101e-5)
            disp('time jump !!!')
          end
          tim = linspace(tim(1),tim(end),length(tim));
          tim2 = linspace(tim(1),tim(end),length(tim)*cond_interp_res);
     
          %
          % interpolate NaNs as they prevent the calculation
          % 
          vars = {'c','t','o','p'};
          for m=1:length(vars)
            if exist(vars{m})
              eval(['bad = find(isnan(',vars{m},'));'])
              eval(['good = find(~isnan(',vars{m},'));'])
              if length(good)>2
                eval([vars{m},'(bad) = interp1(tim(good),',vars{m},...
                  '(good),tim(bad),''linear'',0);'])
              end
            end
          end

  
          %
          % calculate the conductivity delay
          %
          % the old and slow approach was
          % we define the best C delay as the delay at which the salinity 
          % curve which is derived from T and C differs least from a filtered
          % version of the salinity. That is so because any misalignment
          % causes additional spikes in salinity.
          %
          % the new one is just like for oxygen. The best C delay is the delay at
          % which the high frequency variability of T and C has the highest
          % correlation. I.e. quick changes of T and C are correlated.
          % Since these delays are rather small, T and C need to be interpolated
          % to 10 times their original resolution (tim2)
          % 
          [b,a] = butter(3,0.01);
  
          ca = interp1(tim,c,tim2,'spline')';
          ta = interp1(tim,t,tim2,'spline')';

          cam = filtfilt(b,a,ca);
          tam = filtfilt(b,a,ta);
          cad = ca-cam;
          tad = ta-tam;
          [cres,la] = xcorr(cad,tad*signs(1),shift_range);
          cdelay = la/86400/24/cond_interp_res;

          [dummy,ind] = nmax(cres);

          if prisec==1
            delays.pri_c_ind(n) = ind;
            delays.pri_prof(n) = n;
          else
            delays.sec_c_ind(n) = ind;
            delays.sec_prof(n) = n;
          end

          figure(1)
          clf
          subplot(2,1,1)
          plot(cdelay*86400,cres,'x')
          title(['Conductivity-Temp HF lagged crosscorrelation ',str])

          %
          % calculate the oxygen delay
          %
          % we define the best oxygen delay as the delay at which the high
          % frequency variability of oxygen agrees best with the high
          % frequency delay of temperature
          %
          [b,a] = butter(3,0.003);
          pa = p(260:end-260);
    
          om = filtfilt(b,a,o);
          tm = filtfilt(b,a,t);
          od = o-om;
          td = t-tm;
          [ores,la] = xcorr(od,td*signs(2),shift_range);
          odelay = la/86400/24;
  
          [dummy,ind] = nmax(ores);
          if prisec==1
            delays.pri_o_ind(n) = ind;
          else
            delays.sec_o_ind(n) = ind;
          end
          subplot(2,1,2)
          plot(odelay*86400,ores,'x')
          title('Oxygen-Temp HF lagged crosscorrelation')
          pause(0.1)

        end

      end

    end  
    fprintf(1,'\n',[])

    figure(2)
    clf
    if prisec==1
      cresall = delays.pri_c_ind(prof);
    else
      cresall = delays.sec_c_ind(prof);
    end

    good = find(~isnan(cresall));
    [n,x] = hist(cdelay(cresall(good)),cdelay);
    bar(cdelay*86400,n);
    xlabel('time to advance conductivity [s]')
    ylabel('stations')
    title('Conductivity HF variability at time+delay is optimally aligned with temperature')
    ax = axis;
    good = find( cresall>1 & cresall~=length(cdelay) );
    mind = mean( cresall(good) );
    meind = median( cresall(good) );
    text(ax(1)+(ax(2)-ax(1))/2,ax(4)-(ax(4)-ax(3))/5,[str,' sensor'])
    disp([str,' sensor'])
    if prisec==1
      setdelays.pri(se).c = interp1([1:length(cdelay)],cdelay,meind);
    else
      setdelays.sec(se).c = interp1([1:length(cdelay)],cdelay,meind);
    end
    disp(['Average optimal conductivity delay is : ',...
         num2str(86400*interp1([1:length(cdelay)],cdelay,mind)),'s'])
    disp(['Median optimal conductivity delay is : ',...
         num2str(86400*interp1([1:length(cdelay)],cdelay,meind)),'s'])
    text(ax(1)+(ax(2)-ax(1))/2,ax(4)-(ax(4)-ax(3))/10,...
         ['Average optimal conductivity delay is : ',...
         num2str(86400*interp1([1:length(cdelay)],cdelay,mind)),'s'])
    text(ax(1)+(ax(2)-ax(1))/2,ax(4)-(ax(4)-ax(3))/7,...
         ['Median optimal conductivity delay is : ',...
         num2str(86400*interp1([1:length(cdelay)],cdelay,meind)),'s'])

    print('-dpsc',['plots/',crsid,'_c_delays_',int2str(se),'.ps'])

    figure(3)
    clf
    if prisec==1
      oresall = delays.pri_o_ind(prof);
    else
      oresall = delays.sec_o_ind(prof);
    end
    good = find(~isnan(oresall));
    [n,x] = hist(odelay(oresall(good)),odelay);
    bar(odelay*86400,n);
    xlabel('time to advance oxygen [s]')
    ylabel('stations')
    title('Oxygen HF variability at time+delay is optimally aligned with temperature')
    ax = axis;
    good = find( oresall>1 & oresall~=length(odelay) );
    mind = mean( oresall(good) );
    meind = median( oresall(good) );
    text(ax(1)+(ax(2)-ax(1))/2,ax(4)-(ax(4)-ax(3))/5,[str,' sensor'])
    disp([str,' sensor'])
    if prisec==1
      setdelays.pri(se).o = interp1([1:length(odelay)],odelay,meind);
    else
      setdelays.sec(se).o = interp1([1:length(odelay)],odelay,meind);
    end
    disp(['Average optimal oxygen delay is : ',...
         num2str(86400*interp1([1:length(odelay)],odelay,mind)),'s'])
    disp(['Median optimal oxygen delay is : ',...
         num2str(86400*interp1([1:length(odelay)],odelay,meind)),'s'])
    text(ax(1)+(ax(2)-ax(1))/2,ax(4)-(ax(4)-ax(3))/10,...
         ['Average optimal oxygen delay is : ',...
         num2str(86400*interp1([1:length(odelay)],odelay,mind)),'s'])
    text(ax(1)+(ax(2)-ax(1))/2,ax(4)-(ax(4)-ax(3))/7,...
         ['Median optimal oxygen delay is : ',...
         num2str(86400*interp1([1:length(odelay)],odelay,meind)),'s'])

    print('-dpsc',['plots/',crsid,'_o_delays_',int2str(se),'.ps'])

  end

end


%
% distribute the delay results for each setting so that we have one
% delay per profile
%
for prisec = 1:2
  if prisec==1
    settings = pri_settings;
  else
    settings = sec_settings;
  end
  for se=1:length(settings)
    for prof = settings(se).prof;
      if prisec==1
        delays.pri.c(prof) = setdelays.pri(se).c;
        delays.pri.o(prof) = setdelays.pri(se).o;
      else
        delays.sec.c(prof) = setdelays.sec(se).c;
        delays.sec.o(prof) = setdelays.sec(se).o;
      end
    end
  end
end


save6(['mat/',crsid,'_delays'],'delays')

