function misfit = soak_model(params,p)
% function misfit = soak_model(params,p)
%
% misfit function for the determination of the beginning of the downcast
% this function is called by fitsoak.m
%
% input  :  params    - parameters determining the soak model
%           p         - pressure record in dbar
%
% output :  misfit    - rms misfit in m
%
% version 0.5	last change 06.01.2009

% G.Krahmann, IFM-GEOMAR, February 2007

% added penalty for params too close to each other	
%                                           GK, Mar 2007    0.1-->0.2
% modified display for speed                GK, 14.03.2007  0.2-->0.3
% added penalty for out of bounds params    GK, 21.06.2007  0.3-->0.4
% profile number as title                   Gk, 06.01.2009  0.4-->0.5

global plotnew handle_model_graph profnum

misfit = 0;

%
% parse parameters
%
p_deck = params(1);			% pressure on deck
p_soak = params(2);			% pressure during soak
p_start = params(3);			% pressure at the start of the cast (highest point after soak)
params(4:7) = params(4:7)*1000;		% conversion to have all params of the same order of magnitude
p_max = p(end);				% maximum pressure
s_water = round(params(4));		% time (in scans) when the CTD first records higher pressures
if s_water<1				% this can not be before or after the record
  misfit = misfit + abs(s_water);
  s_water = 1;
elseif s_water>length(p)
  misfit = misfit + abs(s_water);
  s_water = length(p);
end

%
% determine the misfit
%
fac = 1;

s_soak_start = round(params(5));	% time (in scans) when the soak period starts
if s_soak_start<1			% this can not be before or after the record
    misfit = misfit + (1-s_soak_start);
    s_soak_start = 1;
elseif s_soak_start>length(p)
    misfit = misfit + (s_soak_start-length(p));
    s_soak_start = length(p);
end

s_soak_end = round(params(6));		% time (in scans) when the soak period ends
if s_soak_end<1				% this can not be before or after the record
    misfit = misfit + (1-s_soak_end);
    s_soak_end = 1;
elseif s_soak_end>length(p)
    misfit = misfit + (s_soak_end-length(p));
    s_soak_end = length(p);
end

s_start = round(params(7));		% time (in scans) when the downcast starts
if s_start<1				% this can not be before or after the record
    misfit = misfit + (1-s_start);
    s_start = 1;
elseif s_start>length(p)
    misfit = misfit + (s_start-length(p));
    s_start = length(p);
end

%
% assemble the soak model (see fitsoak.m for a 'drawing')
%
mp = repmat(0,[length(p),1]);
mp(1:s_water-1) = p_deck;
mp(s_water:s_soak_start) = linspace(p_deck,p_soak,s_soak_start-s_water+1);
mp(s_soak_start:s_soak_end) = p_soak;
mp(s_soak_end:s_start) = linspace(p_soak,p_start,s_start-s_soak_end+1);
mp(s_start:end) = linspace(p_start,p_max,length(mp)-s_start+1);

%
% plot the record and the model (update only every 10th cycle)
%
if floor(plotnew/10)*10==plotnew
  if isnan(handle_model_graph)
    clf
    handle_model_graph = plot(mp);
    hold on
    plot(p,'r')
    title(int2str(profnum))
  else
    set(handle_model_graph,'ydata',mp)
  end
  drawnow
end
plotnew = plotnew+1;

%
% penalized are:
% timing params too close to each other
% soak shorter than 40 seconds (actually likes soak longer than 1500 scans)
% shallow soak (3m gets penalty)
% deck pressure far from first pressure value
%
fac2 = 100;
pen = 0;
pen = pen + exp( -(s_soak_start - s_water)/fac2 );
pen = pen + exp( -(s_soak_end - s_soak_start -700)/fac2 );
pen = pen + exp( -(s_start - s_soak_end)/fac2 );
if p_soak-p_deck < 3
  pen = pen + (3-p_soak+p_deck)*10;
end
pen = pen + (exp( abs(p(1) - p_deck) )-1);
misfit = misfit + pen;

%
% this is the real misfit
%
misfit = misfit + rms(p-mp);
