function [scan,pp] = fitsoak(p,prof,verb)
% function [scan,res] = fitsoak(p,prof,verb)
%
% fit a surface soak to the pressure record and
% determine the scan from which the downcast starts
%
%  input :      p               - pressure record in dbar
%               prof      []    - profile number
%               verb      ''    - set to 'verbose' to get more info
%
%  output :     scan            - index from which the downcast starts
%               pp              - scan model parameters
%                               [deck-pressure,soak-pressure,start-pressure,...
%                                to-water-scan,soak-start-scan,soak-end-scan,...
%                                profile-start-scan];
%
% version 0.9   last change 04.05.2009

% G.Krahmann, IFM-GEOMAR Kiel, January 2007
% changed display, removed bad values         GK, 30.01.2007  0.1->0.2
% optimized start values                      GK, 13.03.2007  0.2-->0.3
% parameter output                            GK, 14.04.2007  0.3-->0.4
% first p can be NaN, interpolate p's         GK, 09.06.2007  0.4-->0.5
% handle Yoyo casts                           GK, 21.06.2007  0.5-->0.6
% added 'verbose' setting                     GK, 03.07.2008  0.6-->0.7
% figure(1), profile number                   GK, 06.01.2009  0.7-->0.8 
% reduce number of points to 1/10             GK, 04.05.2009  0.8-->0.9

% surface soak model has seven parameters which have to be
% determined via non-linear fit
%
% ---\                              - pressure-on-deck  p_deck
%     \             /\              - pressure-start    p_start
%      \           /  \
%       \_________/    \            - pressure-soak     p_soak
%                       \
%    |            |      \
%  s_water    s_soak_end  \
%       |            |     \
%   s_soak-start   s_start  \

global plotnew handle_model_graph profnum
figure(1)
plotnew = 0;
profnum = prof;
handle_model_graph = nan;

%
% make sure we have a column vector and get rid of NaN's
%
p = p(:);
p = nans(p,nan,-9.990e-29,'==');    
good = find(~isnan(p));
p = interp1(good,p(good),[1:length(p)],'nearest');
ind = [1:length(p)];
if good(1)>1
  bad = find(ind<good(1));
  p(ind(bad)) = p(good(1));
end
p = p(:);

%
% check whether we are deep only, e.g. in a yoyo cast
%
if isempty(find(p<30))
  scan = 1;
  pp = [nmean(p(1:10)),0,0,0,0,0,0,0];
  return
end

%
% first delimit the pressure record to data only above 50 dbar
%
ind = find(p>30);
if ~isempty(ind)
  p = p(1:ind(1));
  if strcmp(verb,'verbose')
    disp(['cut record to first ',int2str(ind(1)),...
	' scans (pressure less than 30dbar)'])
  end  
end
p_max = p(end);

%
% special case of a last Yoyo cast
%
if length(p)<10
  scan = 1;
  pp = [p(1),0,0,0,0,0,0,0];
  return
end

%
% set start parameters for the model
% these could be better optimized
%
first_off_deck = find( p-nmean(p(1:20)) > 5 );
first_off_deck = first_off_deck(1);
s_water = first_off_deck-100;
s_soak_start = first_off_deck+100;
s_soak_end = s_soak_start + round( (length(p)-s_soak_start)/2 );
s_start = s_soak_end + round( (length(p)-s_soak_end)/2 );
params = [min(p),min(p)+12,min(p)+5,s_water,s_soak_start,s_soak_end,s_start];
if min(p(1:100))<5
  params = [min(p(1:100)),min(p(1:100))+12,...
            min(p)+5,s_water,s_soak_start,...
            s_soak_end,s_start];
else
  s_water = 1;
  s_soak_start = 100;
  [dummy,s_start] = nmin(p); 
  s_soak_end = round( (s_soak_start+s_start)/2 );
  params = [min(p(1:100)),nmean(p(s_soak_start:s_soak_end)),...
            min(p),s_water,s_soak_start,s_soak_end,s_start];
end

%
% scale them to have the same order of magnitude
%
params(4:7) = params(4:7)/1000;

%
% set some stuff for the optimisation function
%
pp = params;
op = optimset('display','iter','maxfunevals',1000);
op = optimset('maxfunevals',200);
op = optimset('display','off','maxfunevals',1000);
res = 100;
count = 1;
if strcmp(verb,'verbose')
  disp(pp)
end
while res>0.5 & count<=4
    pp = fminsearch('soak_model',pp,op,p);
    res = soak_model(pp,p);
    count = count+1;
    if res<=0.5
      if strcmp(verb,'verbose')
        disp('cutting fit short since fit better than 0.5m')
      end
    end
end
disp(['final soak model misfit : ',num2str(res),' m'])
if strcmp(verb,'verbose')
  disp(pp)
end
disp(' ')

%
% determine highest point between end of soak and 30dbar
%
p = nans(p,nan,0,'==');
[dummy,ind] = min(p(round(pp(6)*1000:end)));
scan = round(pp(6)*1000)+ind-1;

%
% special case that the CTD wasn't raised again
%
if abs(pp(2)-pp(3))<1
  if strcmp(verb,'verbose')
    disp('special case: no raising to surface')
  end
  scan = round(pp(7)*1000);
end

%
% plot result
%
soak_model(pp,p);
hold on
ax = axis;
plot([scan,scan],ax(3:4),'g')
pause(2)
