function datpar2plot(filename,new_figure);
% datpar2plot(filename,make_new_figure)
%
% Reads filename.par, filename.rep and filename.dat and creates a 
% plot of the population trajectory, with one color per cohort. In
% Addition, the survey data divided by q is plotted with the same color as the 
% corresponding cohort, and with the marker for the survey data
% being the same for data from the same survey (but potentially 
% from different cohorts).
%
% Written by Lennart Frimannslund, April 2010

% Get N matrix from filename.rep
extractpar;
% Scale back N to match scaling in cod.tpl
N0scale = 1000;
N = N/N0scale;

% Get stds from filename.std
%clear std_*;
std_filename = [filename '.std'];
%who('std_*');
extractstd;

%whos('std_*')
std_coeff = 1.96;

% Get first year from dat file. 
[filename '.dat']
[m,startyears,endyearsy,C,num_surveys,survey] = read_dat([filename '.dat']);
firstyear = startyears(1);

% Helper arrays for plotting
colors = {'r' 'g' 'b' 'c' 'm' 'y' 'k'};
styles = {'-' '--' '-.' ':'};
dotstyles = {'+' 'o' '*' 'x'};

width = size(N,2);

num_cohorts = size(N,1);

if (exist('new_figure') == 1 && new_figure == 0),
 % No new figure
 fprintf(1,'%s debug: No new figure requested.\n',mfilename);
else,
 figure;
end;
hold on;

% Plot trajectories
fprintf(1,'%s debug: num_cohorts = %i\n',mfilename,num_cohorts);

for (i=1:num_cohorts),
  plot(firstyear+(i-1):firstyear+(i-1) + (width-1) - (i-1),N(i,1: width - (i-1)),colors{i});
  try,
    %errorbar(firstyear+(i-1),N(i,1),std_coeff * std_N0(i,2),colors{i});
    %
    % ERRORBAR doesn't quite work in octave, so make our own bar
    plot([ firstyear+(i-1) ; firstyear+(i-1)],[N(i,1) - std_coeff * std_N0(i,2); N(i,1) + std_coeff * std_N0(i,2)], [colors{i} styles{1}]);

  catch,
    %fprintf(2,'Errorbar plot failed for cohort born in %4i -- wrong variable name?\n',firstyear+(i-1));
    %disp(lasterr);
    errorbar(firstyear+(i-1),N(i,1),std_coeff * std_N0(i,2));
  end;
end;


% Plot observations divided by q
startyears = firstyear:firstyear+num_cohorts-1;

for (i=1:size(survey,1)),
  
  % Check if q_alpha is present, which indicates sigmoid function
  if (exist('q_alpha','var')),
    q_sigmoid = q0 * ((exp(q_alpha*survey(i,2)+q_beta))/(1+exp(q_alpha*survey(i,2)+q_beta)));
    
    plot(startyears(survey(i,1))+survey(i,2), ...
        (1/N0scale)*survey(i,4)/q_sigmoid, ...
        [colors{survey(i,1)} dotstyles{survey(i,3)}]);
  elseif (exist('q02','var')),
    if (survey(i,2) >= 3),
       plot(startyears(survey(i,1))+survey(i,2),(1/N0scale)*survey(i,4)/q3plus(survey(i,3)),...
        [colors{survey(i,1)} dotstyles{survey(i,3)}]);
    else,
       plot(startyears(survey(i,1))+survey(i,2),(1/N0scale)*survey(i,4)/q02(survey(i,3)),...
        [colors{survey(i,1)} dotstyles{survey(i,3)}]);
    end;

  else,
    plot(startyears(survey(i,1))+survey(i,2),(1/N0scale)*survey(i,4)/q(survey(i,3)),[colors{survey(i,1)} dotstyles{survey(i,3)}]); 
  end;

end;

title(sprintf('Cod population trajectory. File: %s(.dat,.rep)',string_convert_underscore(filename)));
xlabel('Year');
ylabel('Volume in billions');

%curr_axis = axis;
%try,
%  axis([curr_axis(1) curr_axis(2) 0 2500/N0scale curr_axis(5) curr_axis(6)]);
%catch,
%  axis([curr_axis(1) curr_axis(2) 0 2500/N0scale]);
%end;
hold off;

% If q is a sigmoid function, plot this function
if (exist('q_alpha','var')),

sf = figure;
xx = 0:1:20;
yy = q0 * (exp(q_alpha*xx+q_beta)./(1+exp(q_alpha*xx+q_beta)));
plot(xx,yy);
title('Sigmoid function q0 * exp(q_a x + q_b)/(1+exp(q_a x + q_b))');
xlabel('Age');
ylabel('Catchability');

end;


return;

% Convert _ to \_ so that latex interpretation doesn't garble the string
function [out] = string_convert_underscore(input);

 num_underscores = length(find(input == '_'));
 if (num_underscores == 0),
   out = input;
 else,

   out = char(32*ones(length(input)+num_underscores,1));
   counter = 1;
   for (i=1:length(input)),
     if (input(i) == '_'),
       out(counter:counter+1) = '\_';
       counter = counter + 2;
     else,
       out(counter) = input(i);
       counter = counter + 1;
     end;
   end;
 end;

return;