% model used for all runs shown in the paper


clear all
close all
fclose all;

start_day    = 1;         % starting day
prey_size_min= 100;       % maximum prey size available (prey size in µm); smallest is alway 100
prey_size_max= 1500;      % maximum prey size available (prey size in µm); smallest is alway 100
class_size   = 100;
%class_size   = (prey_size_max-prey_size_min)/15; % prey class size (in µm)
T_start      = 2;         % Start temperature for temperature variation
T_end        = 20;        % End temperature for temperature variation
T_class      = 1;         % Temperature class size (step size between the variations) in °C
L_max        = 20;        % maximum length where simulation is stopped in mm
n_conc       = 15;        % amount of concentrations for the prey concentration variation (after eych step start conc is doubled)
day_max      = 365;       % end day when simulation is stopped
prey_conc_all= 0.0001;     % lowest concentration of prey in no/ml (sum of all prey sizes) for the concentrations variations
k_day        = 2;         % maximum k for active respiration (at the moment k = k_max)
k_night      = 1;
Q10GER       = 2.52;      % Q10 Gut evacuation
GER_incr     = 2.0;       % Faktor to multiply gut evacuation rate under constant feeding
CS_b         = 1.1;       % Capture success parameter (the lower the easier prey can be caught)
w            = 1.3;       % turbulent velocity [mm/s] from Dower et al 2002
prop_fac     = 0;         % if not equal 0 then this value determines the amount of food that is channeled into length growth 1=all in length growth
HT_a         = 0.264;     % 0.264;
HT_b         = 20;
PD_start     = 1.3;       % 1.3  pause duration mean from Myrons data > 9mm  < 13 mm
PF_start     = 0.35;      % 0.35 pause frequency 
max_p_a      = 2200;      % 1524.9
max_p_b      = 2;         % 1.63
L_after_y0   = 11;        % 9.75 Length after the yolk sac phase reduced by temperature
night_length = 10;
pi = 4*atan(1);
sc_day       = 0;        % slope and concentration julian day (0 for ignoring)

%Gauss_sd     = 100;
%Gauss_m      = 400;
%Gauss_sd     = 43;
%Gauss_m      = 280;
%Gauss_sd2    = 430;
%Gauss_m2     = 35;

% check box GUI for plots
% creating window with check boxes
cb1  = uicontrol('Style', 'checkbox', 'String', 'gain Rs R','Position', [20 400 400 20], 'Max', 1, 'Min',0);
cb2  = uicontrol('Style', 'checkbox', 'String', 'development time (contour)','Position', [20 380 400 20], 'Max', 1, 'Min',0);
cb3  = uicontrol('Style', 'checkbox', 'String', 'length over time','Position', [20 360 400 20], 'Max', 1, 'Min',0,'Value',0);
cb4  = uicontrol('Style', 'checkbox', 'String', 'length growth rate over T','Position', [20 340 400 20], 'Max', 1, 'Min',0,'Value',1);
cb5  = uicontrol('Style', 'checkbox', 'String', 'DW growth over length','Position', [20 320 400 20], 'Max', 1, 'Min',0);
cb6  = uicontrol('Style', 'checkbox', 'String', 'DW growth over food concentration (+ Kiörboe)','Position', [20 300 400 20], 'Max', 1, 'Min',0,'Value',1);
cb7  = uicontrol('Style', 'checkbox', 'String', 'GGE over length','Position', [20 280 400 20], 'Max', 1, 'Min',0,'Value',0);
cb8  = uicontrol('Style', 'checkbox', 'String', 'DW growth over ingestion (GGE Kiörboe)','Position', [20 260 400 20], 'Max', 1, 'Min',0,'Value',1);
cb9  = uicontrol('Style', 'checkbox', 'String', 'DW ingested over length','Position', [20 240 400 20], 'Max', 1, 'Min',0);
cb10 = uicontrol('Style', 'checkbox', 'String', 'DW ingested over fish DW','Position', [20 220 400 20], 'Max', 1, 'Min',0);
cb11 = uicontrol('Style', 'checkbox', 'String', 'Particles ingested over length','Position', [20 200 400 20], 'Max', 1, 'Min',0);
cb12 = uicontrol('Style', 'checkbox', 'String', 'Feeding stike frequ over length and conc (contour)','Position', [20 60 400 20], 'Max', 1, 'Min',0);
cb13 = uicontrol('Style', 'checkbox', 'String', 'preferred prey size over length','Position', [20 40 400 20], 'Max', 1, 'Min',0);
cb14 = uicontrol('Style', 'checkbox', 'String', 'DW over length highest conc','Position', [20 140 400 20], 'Max', 1, 'Min',0,'Value',0);
cb15 = uicontrol('Style', 'checkbox', 'String', 'temperature history','Position', [20 120 400 20], 'Max', 1, 'Min',0,'Value',0);
cb16 = uicontrol('Style', 'checkbox', 'String', 'DW over length all conc','Position', [20 100 400 20], 'Max', 1, 'Min',0,'Value',0);
cb17 = uicontrol('Style', 'checkbox', 'String', 'development time over temperature','Position', [20 80 400 20], 'Max', 1, 'Min',0,'Value',0);
cb18 = uicontrol('Style', 'checkbox', 'String', 'length growth over food concentration (C)','Position', [20 180 400 20], 'Max', 1, 'Min',0,'Value',1);
cb19 = uicontrol('Style', 'checkbox', 'String', 'dry weight over time','Position', [20 160 400 20], 'Max', 1, 'Min',0,'Value',0);
cb20 = uicontrol('Style', 'checkbox', 'String', 'PLOT ALL FIGURES','Position', [150 20 400 20], 'Max', 1, 'Min',0,'Value',0);
cb21 = uicontrol('Style', 'checkbox', 'String', 'starvation: dry weight over degree days','Position', [300 400 400 20], 'Max', 1, 'Min',0,'Value',0);
cb22 = uicontrol('Style', 'checkbox', 'String', 'Volume searched per day','Position', [300 380 400 20], 'Max', 1, 'Min',0,'Value',0);
cb23 = uicontrol('Style', 'checkbox', 'String', 'Particles in the gut','Position', [300 360 400 20], 'Max', 1, 'Min',0,'Value',0);
cb24 = uicontrol('Style', 'checkbox', 'String', 'Particles encountered per day','Position', [300 340 400 20], 'Max', 1, 'Min',0,'Value',0);
cb25 = uicontrol('Style', 'checkbox', 'String', 'Reactive distance [mm]','Position', [300 320 400 20], 'Max', 1, 'Min',0,'Value',0);
cb26 = uicontrol('Style', 'checkbox', 'String', 'Percent body weight eaten per day over length','Position', [300 300 400 20], 'Max', 1, 'Min',0,'Value',0);
cb27 = uicontrol('Style', 'checkbox', 'String', 'Length over time at 8°C and different conc. + Lit','Position', [300 280 400 20], 'Max', 1, 'Min',0,'Value',0);
cb28 = uicontrol('Style', 'checkbox', 'String', 'DW over time at 8°C and different conc. + Lit','Position', [300 260 400 20], 'Max', 1, 'Min',0,'Value',0);
cb29 = uicontrol('Style', 'checkbox', 'String', 'Ingestion [µg d] over food dens [µg l]','Position', [300 240 400 20], 'Max', 1, 'Min',0,'Value',0);
cb30 = uicontrol('Style', 'checkbox', 'String', 'Clearance rate (particles consumed/encounter rate)','Position', [300 240 400 20], 'Max', 1, 'Min',0,'Value',0);
cb31 = uicontrol('Style', 'checkbox', 'String', 'growth per day vs. ingestion per day','Position', [300 220 400 20], 'Max', 1, 'Min',0,'Value',1);
cb32 = uicontrol('Style', 'checkbox', 'String', 'length growth per day vs. age','Position', [300 200 400 20], 'Max', 1, 'Min',0,'Value',0);
cb33  = uicontrol('Style', 'checkbox', 'String', 'length over time with literature (winter)','Position', [300 180 400 20], 'Max', 1, 'Min',0,'Value',0);
cb34  = uicontrol('Style', 'checkbox', 'String', 'length over time with literature (summer)','Position', [300 160 400 20], 'Max', 1, 'Min',0,'Value',0);
cb35  = uicontrol('Style', 'checkbox', 'String', 'dry weight at length ','Position', [300 140 400 20], 'Max', 1, 'Min',0,'Value',0);

bstart = uicontrol('Style', 'pushbutton', 'String', 'start', 'Position', [20 20 100 20], 'Callback', 'uiresume');
 

% wait until all check boxes are set
uiwait
%get values from the boxes (only needed in the figure section)
val1 = get(cb1,'Value');
val2 = get(cb2,'Value');
val3 = get(cb3,'Value');
val4 = get(cb4,'Value');
val5 = get(cb5,'Value');
val6 = get(cb6,'Value');
val7 = get(cb7,'Value');
val8 = get(cb8,'Value');
val9 = get(cb9,'Value');
val10 = get(cb10,'Value');
val11 = get(cb11,'Value');
val12 = get(cb12,'Value');
val13 = get(cb13,'Value');
val14 = get(cb14,'Value');
val15 = get(cb15,'Value');
val16 = get(cb16,'Value');
val17 = get(cb17,'Value');
val18 = get(cb18,'Value');
val19 = get(cb19,'Value');
val20 = get(cb20,'Value');
val21 = get(cb21,'Value');
val22 = get(cb22,'Value');
val23 = get(cb23,'Value');
val24 = get(cb24,'Value');
val25 = get(cb25,'Value');
val26 = get(cb26,'Value');
val27 = get(cb27,'Value');
val28 = get(cb28,'Value');
val29 = get(cb29,'Value');
val30 = get(cb30,'Value');
val31 = get(cb31,'Value');
val32 = get(cb32,'Value');
val33 = get(cb33,'Value');
val34 = get(cb34,'Value');
val35 = get(cb35,'Value');

close all
% if all graphs shall be plotted at the end
if (val20)
val2 = 1;
val3 = 1;
val4 = 1;
val5 = 1;
val6 = 1;
val7 = 1;
val8 = 1;
val9 = 1;
val10 = 1;
val11 = 1;
val14 = 1;
val16 = 1;
val17 = 1;
val18 = 1;
val19 = 1;
val20 = 1;
val22 = 1;
val23 = 1;
val24 = 1;
val25 = 1;
val26 = 1;
val27 = 1;
val28 = 1;
val29 = 1;
val30 = 1;
val31 = 1;
val32 = 1;
val33 = 1;
val34 = 1;
val35 = 1;

end 
% initializing random with clock 
rand('state',sum(100*clock))
% loop for varying concentrations
for conc_count = 1:n_conc
   count_T = 0; % counter for writing the temperature variations  
   clear conc
   prey_classes = 0;
   for (n=prey_size_min:class_size:prey_size_max)
       prey_classes = prey_classes + 1;    % conting the of prey classes
   end
   % If slope and concentration should vary by day according to a sinus
   % function
   if (sc_day)
          slope = -2.5966-0.5527*sin(-0.017214*(sc_day+97.762));
          slope = -2;
          conc(1:prey_classes) = 100000.*((prey_size_min:class_size:prey_size_max).*100).^slope;
          prey_conc_all = (0.1034+0.02896*sin(-0.01721*(sc_day+46)));
          prey_conc_all = 0.1;
          prey_conc_class(1:prey_classes) = conc/sum(conc)*prey_conc_all;   
   else

%FUNCTION USED FOR PAPER RUNS 
      conc(1:prey_classes) = 695.73*exp(-0.0083.*(prey_size_min:class_size:prey_size_max)); % concentration in prey class determined by exponetial decrease (Daewel & Peck)
%%%%%

%   conc(1:prey_classes) = 100000.*((prey_size_min:class_size:prey_size_max).*100).^(-2);  % constant slope of 2
%   conc(1:prey_classes) = 13.7*exp(-0.0045.*(100:class_size:prey_size_max));              % from Lenis diploma thesis

% concentration in prey class as Gauss function (used for paper scenario runs
%   conc(1:prey_classes) =   (((1/(Gauss_sd*(2*3.14)^0.5))^2).*exp(-(((prey_size_min:class_size:prey_size_max)-Gauss_m).^2)./(2*Gauss_sd^2))); 
%   conc(1:prey_classes) = conc(1:prey_classes)+(((1/(Gauss_sd2*(2*3.14)^0.5))^2).*exp(-(((prey_size_min:class_size:prey_size_max)-Gauss_m2).^2)./(2*Gauss_sd2^2))); % concentration in prey class Gauss function

      prey_conc_class(1:prey_classes) = conc/sum(conc)*prey_conc_all;   % transferring exp curve in precantage values and multiplying with maximum prey concentration

   end    
   % loop for temperature variations
   %Temperature=[10.5685	10.5685	10.5178	10.5178	10.5685	10.5685	10.7716	10.7716	10.264	10.264	10.6701	10.6701	10.1624	10.2132	10.2132	9.7563	9.7563	8.9949	8.9949	9.7056	10.0609	10.0609	10.0609	10.264	10.264	11.8883	11.8883	13.2589	14.1218	14.1218	14.1218	13.8173	13.8173	13.7157	13.7157	14.0711	14.0711	14.0711	14.0711	14.0711	14.0711	14.2741	14.5787	14.5787	14.5787	15.4416	15.9492 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16];
   for temp_T = T_start:T_class:T_end 

      L_after_y = L_after_y0 -0.09 *temp_T;       % start length (after yolk phase) according to temperature  
      disp(sprintf('temperature: %d  \t concentration: %f  \tconc %d of %d',temp_T,prey_conc_all,conc_count,n_conc))
      clear length; 
      clear DW;
      clear GC;
      clear T_hist
      clear food
      % initializing values 
      time  = 1;  % counter for time in [h]
      day   = 0;  % counter for time in [d]
      count_T = count_T + 1;  % counter for temperature
      dead = 0;
      % initializing variables
      food_d(count_T,conc_count,1:day_max)   = 0;
      x(count_T,conc_count,1:day_max)        = 0;
      length_d(count_T,conc_count,1:day_max) = 0;
      DW_d(count_T,conc_count,1:day_max)     = 0;      
      Rs_d(count_T,conc_count,1:day_max)     = 0;      
      R_d(count_T,conc_count,1:day_max)      = 0;      
      gain_d(count_T,conc_count,1:day_max)   = 0;
      part_d(count_T,conc_count,1:day_max)   = 0;
      pos_cons_d(count_T,conc_count,1:day_max)= 0;
      pos_cons(1:1:day_max*24)               = 0;
      part_in_gut(1:1:day_max*24)            = 0;
      vol_searched(1:1:day_max*24)           = 0;
      vol_searched_d(count_T,conc_count,1:day_max) = 0;
      length(1:day_max*24) = 0;
      DW(1:day_max*24)     = 0;
      GC(1:day_max*24)     = 0;   
      food(1:1:day_max*24) = 0;
      
% including egg development (time of peak hatch from Myron)
      fract = 1; % if 1 than no egg development included (else 0)
      day         = 1+floor(time/24);
      T_hist(day) = temp_T;
      while (fract < 1)
         fract       = fract+(1/(4461.9*temp_T^(-1.232)));
         day         = 1+floor(time/24);
         T_hist(day) = temp_T;
         time        = time +1;
      end
      length(time) = -0.0023*temp_T^3+0.084*temp_T^2-1.035*temp_T+10.77;  % length at hatch [mm] Myrons data
      DW(time)     =  0.018521  * L_after_y^(3.613981*exp(mean(T_hist)*0.006266))*1.3;   % dry weight according to length and temperature good condition assumed    
      length(1:time) = length(time); 
      length_d(count_T,conc_count,1:day) = length(time);    
      hatchd(count_T,conc_count) = day;
      hatchl(count_T,conc_count) = length(time);
      hatchdw(count_T,conc_count)= 90; %hatch dry weight set to 90 (used only for plotting)
      
% including yolk sac phase      
% start of feeding according to RNA/DNA ratio after 59 degree days
      fract = 0;
      while (fract < 1)
         fract = fract + 1/((59*24)/temp_T);
         time  = time +1;
         day   = 1+floor(time/24);
         T_hist(day) = temp_T;
%         length (time) = L_after_y; % length at yolk end 
         length(time) = length(time-1);
         length_d(count_T,conc_count,day) = length(time);      
%         DW(time)  = 0.018521  * length(time)^(3.613981*exp(temp_T*0.006266));     % dry weight 
         DW(time) = DW(time-1);
         DW_d(count_T,conc_count,day)= DW(time);
      end
      day   = 1+floor(time/24);      
      length_d(count_T,conc_count,day) = length(time);      
      DW_d(count_T,conc_count,day)= DW(time);
%      DW(1:time)=DW(time);
      length(time) = L_after_y;
      DW(time)     =  0.018521  * length(time)^(3.613981*exp(mean(T_hist)*0.006266));   % dry weight in µg
      yolkd(count_T,conc_count) = day+1;         % day when yolk phase ends
      yolkl(count_T,conc_count) = length(time);  % length at the end of the yolk phase
      yolkdw(count_T,conc_count) = DW(time);     % dry weight at the end of the yolk phase
      x(count_T,conc_count,1:day+1) = 1:day+1;      
      if (day > day_max)
          disp ('error: day_max too short')
      end

      % growth loop: Larvae growth until Lmax or day max is reached or larvae dies
      deathday(count_T,conc_count)=day_max; 
      T_hist(1:day) = temp_T;
      while ((length(time) < L_max) & (dead==0) & day < day_max)
         time  = time + 1;  % time in [h]
         count = floor(time/24)+1;
         day = start_day + floor(time/24); % time in [d]
         T_hist(day) = temp_T;         
         if ((time - (count-1) * 24) < night_length)           % dividing into day and night 
            night = 1;
         else 
            night = 0;
         end
         beta   = 0.60 *(1-0.3*exp(-0.003*(DW(time-1)- yolkdw(count_T,conc_count))));    % assimilation efficiency (dimensionless)
%         beta   = beta * exp(-0.2*log(prey_conc_all+1));
         Rs = 0.00463 * 227.0 * (1.3921 *(DW(time-1)/1000)^0.7722*exp(0.05369*temp_T));  % Respiration in µg DW per fish and h
         Rs_mulO2 = (1.3921 *(DW(time-1)/1000)^0.7722*exp(0.05369*temp_T));  % Respiration in µl O2/h and fish
         Card_out = 0.6337*Rs_mulO2;
         Card_Vol = 0.00001326 *  DW(time) - 0.002751;
         DWref =  0.018521  * length(time-1)^(3.613981*exp(mean(T_hist)*0.006266)); % reference dry weight [µg] needed for gut size
         GCmax = 10^(1.72+1.02*log10(DWref/1000)) ; %µg  Pepin and Penny
%         GCmax =(10^(1.29+0.66*log10(DWref/1000)))+GC_fact*DWref;           % Pepin and Penny (wrong values from the paper) plus percentage
%         GCmax = 0.064 * DWref^0.995 ; %µg

         GER    = (1.792 * length(time-1)^(-0.828)*Q10GER^((temp_T-12)/10))*100;     % % of whole Volume per h generealized model for different species from Myron
%         GER    = (2.0755 * length(time-1)^(-0.5680)*Q10GER^((temp_T-12)/10))*100;  % % of whole gut Volume per h (herring data)
%         GER = 0.5457*length(time-1)^(-0.67347)*exp(0.069314*temp_T)*100;%          % of whole gut Volume per h (too high)
%         GER = (1.3582 *length(time-1)^(-0.8281)*Q10GER^((temp_T-8)/10))*100;       % % general function reestimated with T = 8°C
%         GER    = 100-100*exp(-(1.792 * length(time-1)^(-0.828)*Q10GER^((temp_T-12)/10)));     % % of content per h generealized model for different species from Myron

         if (food(time-1)>0)
             GER = GER*GER_incr*(1+(1/(1+exp(2-5*prey_conc_all)))); % increasing GER during constant feeding according to prey concentration
%            GER = GER*GER_incr; % GER increased during constant feeding by the factor GER_incr
%            GER =  GER*(prey_conc_all * 0.05 +GER_incr);
%            GER =  GER*GER_incr;
         end
         if (night)         % at night only standard metabolism and no consumption
            k = k_night;
            consum  = 0;
         else
            % coupling activity to prey concentration and starvation
            k = k_day * 0.812 * prey_conc_all^-0.1922;  % according to prey concnetration 
            k = k*(-2*(GC(time-1)/GCmax)^2 + 2*(GC(time-1)/GCmax)+k_night/k) ; % further according to gut fullness
            speed        = 127.21/(  1+exp(-(length(time-1)-30.566)/3.67) )+6.7; %larval swimming speed [mm/s]
            PS_max       = max_p_a / (1+(length(time-1)/14)^(-max_p_b)); %maximum prey size in [µm]
            alpha        = 0.0167*exp(9.14-2.4*log(length(time-1))+0.229*(log(length(time-1)))^2); % angle of visual acuity
            alpha        = alpha * pi/180;
            PF = PF_start; %1/s
            PD = PD_start; %s
%            PD = max(PD_start,8.4715*exp(-0.1608*length(time)));

            % calculations needed for prey size ranking
            for c_cl = 1: prey_classes
               prey_size          = c_cl * class_size - (class_size/2);       %prey size in the middle of the class µm
               prey_weight(c_cl)  = (10^(2.7722*log10(prey_size)-7.47607))*2; % 2* µg Carbon content of prey (Myrons literature research)    
               prey_speed         = prey_size/1000*3;                         % mm/s 3 times body length
               RD(c_cl)           = (prey_size/1000)/(2*tan(alpha/2));        % reactive distance [mm]
               CS(c_cl)           = 1.1-(CS_b*prey_size/PS_max);              % capture success µm / µm
%               CS(c_cl)          = 1.1-13.6*((prey_size/1000)/length(time-1)); %Munk 1992 ohne max prey
               if (CS(c_cl)<0) CS(c_cl)=0; end
               HT(c_cl)           = exp(HT_a*10^(HT_b*((prey_size/1000)/length(time-1))));   % handling time [s]
%               HT(c_cl)           = exp(HT_a*10^(7.0151*((prey_size/1000)/length(time-1))));   % handling time [s]
               rk(c_cl)           = (prey_weight(c_cl)*CS(c_cl))/HT(c_cl);                      % prey type ranking  µg/s
               V(c_cl)            = sqrt(prey_speed^2+2*w^2+speed^2);            % Velocity component of contact rate mm/s
               ER_s(c_cl)         = k*((2/3)*pi*RD(c_cl)^3*(prey_conc_class(c_cl)/1000)*PF)+(pi*RD(c_cl)^2*(prey_conc_class(c_cl)/1000)*PF*PD*V(c_cl));  % Encounter rate no/sec              
%               ER_s(c_cl)         = 0.5*pi*RD(c_cl)^2*(prey_conc_class(c_cl)/1000)*(1-PF*PD)*V(c_cl);  % Encounter rate no/sec cruise predators
            end

            % sort classes from lowest to highest priority
            clear rkind;
            clear rkind2;
            for c_cl = 1: prey_classes 
            rk_count = 1;
               for c_cl2 = 1: prey_classes
                   if (rk (c_cl)<rk (c_cl2))
                       rk_count = rk_count +1;
                   end
               end
               rkind2(c_cl) = rk_count;
            end
            for c_cl = 1: prey_classes
               rkind(c_cl) = 999;                
            end
            for c_cl = 1: prey_classes
               if (rk(c_cl)>0)
                  rkind(rkind2(c_cl)) = c_cl;                
               end
            end
            best_prey = rkind(1);   

            % Calculating sums from lowest to highest rank
            summe1    = 0;
            summe2    = 0;
            summe_alt = 0;
            endloop = 0;
            for c_cl = 1: prey_classes
               if (endloop == 0)
                  if (rkind(c_cl)~=999)
                     summe1 = summe1 + (ER_s(rkind(c_cl)) * prey_weight(rkind(c_cl))*CS(rkind(c_cl))); % n/s * µg
                     
                     % % % Here is an error I presume (Anna Akimova, 20180409) - there should not be 1 in the sum
                     summe2 = summe2 + (1+ER_s(rkind(c_cl))*HT(rkind(c_cl))); % n/s * s = n
                     summe  = summe1/(1+summe2); %µg / s
                     if (summe<summe_alt & endloop == 0)
                        endloop = 1;
                     end
                     summe_alt = summe;
                  end 
               end
            end
            consum       = summe_alt*60*60;          % consumption in µg DW per hour
            pos_cons(time)     = consum/prey_weight(best_prey);   % maximum possible amount of particles ingested per hour
            vol_searched(time) = (((2/3)*pi*RD(best_prey)^3*PF)+(pi*RD(best_prey)^2*PF*PD*V(best_prey)))*(60*60)/(1000*1000); % searched volume per l/h
            part_in_gut(time)  = GC(time-1)/prey_weight(best_prey); % amount of particles actually in the gut n
         end
%         GC(time) = GC(time-1) - (GCmax * GER/100);  %if stomach is full than clear with evacuation rate and then add consumption
         GC(time) = GC(time-1) - (GC(time-1) * GER/100);  %if stomach is full than clear with evacuation rate and then add consumption
         if (GC(time) < 0)
             GC(time) = 0;
         end
         GC(time) = GC(time) + consum;  %µg + µg/h * h
         if (GC(time) > GCmax)
            GC(time)  = GCmax;                       % if stomach is full no more food can be added
         end
         food(time)  = GC(time-1)* GER/100;    % exponential evacuation (dependend on content); the amount evacuated is available for energy gain  µg
%         food(time)  = GCmax * GER/100;    % linear evacuation (dependent on gut size); the amount evacuated is available for energy gain  µg
         if (food(time) > GC(time))
             food(time) = GC(time);
         end
         DW_T   = 0.018521  * length(time-1)^(3.613981*exp(mean(T_hist)*0.006266));       % reference dry weight at length µg
         CI     = DW(time-1)/DW_T;
         R = k*Rs;                                             % active metabolism µg DW per h 
         SDA = 0.10;
%         SDA      = 0.11+ 4.91*10^-6 * DW(time-1);            % specific dynamic action     
         dw_growth = food(time) * beta *(1-SDA)-R;             % amount of dry weight usable for growth in µg
         gain = food(time)* beta *(1-SDA);                    % weight gain without respiration losses (brutto)
         part_h = food(time)/prey_weight(best_prey);           % ~ Particles ingested per hour (best on the best prey size)
         if (dw_growth > 0)
%            prop_length_growth  = 0.5+0.36*atan(0.25*(CI*360-360)/(2*pi)); % proportion of the available DW that is used for length growth at upDW 100 % at lowDW 0%
            prop_length_growth = (CI-0.65)*(1/(1.35-0.65)); % proportion of the available DW that is used for length growth at upDW 100 % at lowDW 0%
%            prop_length_growth = 2.01*CI^3-8.08*CI^2+11.29*CI-4.48; % proportion of the available DW that is used for length growth at upDW 100 % at lowDW 0%
%            prop_length_growth  = 0.59+0.41*atan(0.25*(CI*360-420)/(2*pi)); % proportion of the available DW that is used for length growth at upDW 100 % at lowDW 0%
            if (prop_fac)
               prop_length_growth = prop_fac; % can be defined in the beginning as a fixed value
            end
            % dry weight and length changes according to growth 
            DW(time) = DW(time-1)+dw_growth*prop_length_growth;
            L_incr = (DW(time)/CI/0.018521)^(1/(3.613981*exp(mean(T_hist)*0.006266)))-length(time-1);
            length(time) = length(time-1)+L_incr;
            % no shrinking possible
            if (length(time) < length(time-1))
               length(time) = length(time-1);
            end
            DW(time)= DW(time) + dw_growth * (1-prop_length_growth);    % addition of weight 
         else
            % dry weight loss due to starvation no length loss
            DW(time)= DW(time-1) + dw_growth;
            length(time) = length(time-1);
         end
         if (DW(time)< 0.65*(0.018521*length(time)^(3.613981*exp(mean(T_hist)*0.006266))))
            dead = 1; % animal dies
            length(time) = length(time-1);               % animal dies length set to the previous value     
            deathday(count_T,conc_count)=day;            % saving death day    
         end
         length_d(count_T,conc_count,count) = length(time);              % length per day for plotting
         DW_d(count_T,conc_count,count)     = DW(time);                  % DW per day for plotting
         pos_cons_d (count_T,conc_count,count)= pos_cons_d(count_T,conc_count,count) + pos_cons(time); % max amount that can be consumed
         vol_searched_d(count_T,conc_count,count) = vol_searched_d(count_T,conc_count,count)+ vol_searched(time); % searched volume
         food_d(count_T,conc_count,count)   = food_d(count_T,conc_count,count)+food(time); % daily food ingestion
         part_d(count_T,conc_count,count)   = part_d(count_T,conc_count,count)+part_h; % daily particle ingestion
         Rs_d(count_T,conc_count,count)     = Rs_d(count_T,conc_count,count)+Rs; % daily respiration (routine)
         R_d(count_T,conc_count,count)      = R_d(count_T,conc_count,count)+R;   % daily respiration (active)
         gain_d(count_T,conc_count,count)   = gain_d(count_T,conc_count,count)+gain; % energy gain without respiration 
         x(count_T,conc_count,count)        = day;  
      end
      % for contour plot of development days per conc and T
      Erg(conc_count,count_T) = day; 
      x_t(count_T)            = mean(T_hist);
      y_conc(conc_count)      = prey_conc_all;
%      x_t(count_T)            = temp_T;
%      y_conc(conc_count)      = max_p_b;

      % if calculation was stopped before day_max (animal dies or reaches Lmax)all remaining fields are
      % set to the previous value
      if (day<day_max)
         length_d(count_T,conc_count,count:day_max) = length(time-1);              % length per day for plotting
         DW_d(count_T,conc_count,count:day_max)     = DW(time-1);                  % length per day for plotting
         food_d(count_T,conc_count,count:day_max)   = food_d(count_T,conc_count,count);
         part_d(count_T,conc_count,count:day_max)   = part_d(count_T,conc_count,count);
         x(count_T,conc_count,count:day_max)        = day;  
         deathday(count_T,conc_count)               = day;
      end          
   end
   prey_conc_all = prey_conc_all *2;  % increasing prey concentration for the next step 
%   prey_conc_all = prey_conc_all + 0.005;  
%   max_p_b = max_p_b + 0.1; 
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% FIGURES %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% dw for growth Rs and R according to length and mean over temperature
if(val1)
count_T = 1;  
conc_count = n_conc;
clear datay1;
clear datay2;
clear datay3;
for temp_T = T_start:T_class:T_end
   clear tempy1;
   clear tempy2;
   clear tempy3;
   clear datax;
   for n=1:size(gain_d,3)
      if ((Rs_d(count_T,conc_count,n)~=0) & (length_d(count_T,conc_count,n)<15))
         datax(n)  = length_d(count_T,conc_count,n);
         tempy1(n) = gain_d(count_T,conc_count,n); 
         tempy2(n) = Rs_d(count_T,conc_count,n); 
         tempy3(n) = R_d(count_T,conc_count,n); 
      end
   end
   datay1 (count_T)= mean(tempy1);
   datay2 (count_T)= mean(tempy2);
   datay3 (count_T)= mean(tempy3);
   count_T = count_T+1;
   figure 
   hold on
   plot (datax,tempy1,'k')
   plot (datax,tempy2,'r')
   plot (datax,tempy3,'b')
   legend('gain','Rs','R')   
   title (num2str(temp_T))
end
figure 
hold on
plot (x_t,datay1,'k');
plot (x_t,datay2,'r');
plot (x_t,datay3,'b');
xlabel('temperature');
ylabel('gain, Rs and R');
legend('gain','Rs','R')    
title ('mean of whole model runtime')
end

% figure contour plot for duration time in days
if(val2)
figure
contourf(x_t,y_conc,Erg,30);
colorbar;
colormap jet
xlabel('temperature');
ylabel('prey concentration');
title('development time until reaching 35 mm');
end

if(val17)
figure
hold on
for cx = 1:size(y_conc,2)
   plot (x_t, Erg(cx,:))
   text (x_t, Erg(cx,:),num2str(y_conc(cx)) )
end
xlabel('temperature');
ylabel('time needed to reach maximum length');
end

% figure length growth rate against temperature
if(val4)
figure
hold on
clear mean_l_growth
for conc_count = 1:n_conc
   count_T = 1;    
   for temp_T = T_start:T_class:T_end
      clear length_growth; 
      n = 1;
      for c =  (yolkd(count_T,conc_count))+2:deathday (count_T,conc_count)-1
         length_growth(n) = (length_d(count_T,conc_count,c)-length_d(count_T,conc_count,c-1));
         n = n+1;
      end
      if (deathday (count_T,conc_count)-1<=yolkd(count_T,conc_count)+2)
          length_growth(n)=0;
      end
      mean_l_growth(count_T)  = mean(length_growth);
      count_T = count_T+1;
   end
   plot (x_t,mean_l_growth)
   text (x_t,mean_l_growth,num2str(y_conc(conc_count)))
   xlabel('temperature [°C]')
   ylabel('growth rate (mm/d)')
   title (strcat ('max prey size: ',num2str(prey_size)))
   load T_GR_lit.txt
   plot (T_GR_lit(:,1),T_GR_lit(:,2),'x')
   load T_GR_lit2.txt
   plot (T_GR_lit2(:,1),T_GR_lit2(:,2),'+')
   plot((4:20),[0.16 0.18 0.24 0.25 0.35 0.32 0.35 0.38 0.44 0.50 0.55 0.58 0.65 0.66 0.56 0.57 0.66],'--rs');
end
end

% figure length growth rate against temperature contour plot
if(val4)

clear mean_l_growth
for conc_count = 1:n_conc
   count_T = 1;
   for temp_T = T_start:T_class:T_end
      clear length_growth; 
      n = 1;
      for c =  (yolkd(count_T,conc_count))+2:deathday (count_T,conc_count)-1
         length_growth(n) = (length_d(count_T,conc_count,c)-length_d(count_T,conc_count,c-1));
         n = n+1;
      end
      if (deathday (count_T,conc_count)-1<=yolkd(count_T,conc_count)+2)
          length_growth(n)=0;
      end
      mean_l_growth(count_T,conc_count,1)  = mean(length_growth);
      count_T = count_T+1;
   end
end
figure
mesh (x_t,y_conc,mean_l_growth')
hold on
xlabel('temperature [°C]')
zlabel('growth rate (mm/d)')
ylabel('concentration')
title (strcat ('max prey size: ',num2str(prey_size)))
load T_GR_lit.txt
T_GR_lit(1:size(T_GR_lit,1),3)=y_conc(conc_count);
plot3 (T_GR_lit(:,1),T_GR_lit(:,3),T_GR_lit(:,2),'o')
load T_GR_lit2.txt
T_GR_lit2(1:size(T_GR_lit2,1),3)=y_conc(conc_count);
plot3 (T_GR_lit2(:,1),T_GR_lit2(:,3),T_GR_lit2(:,2),'+')
plot3((4:20),(4:20)./(4:20)*y_conc(conc_count),[0.16 0.18 0.24 0.25 0.35 0.32 0.35 0.38 0.44 0.50 0.55 0.58 0.65 0.66 0.56 0.57 0.66],'--rs');
save concTgrowth.dat  mean_l_growth -ascii
end

% figure length against time for different temperatures and concentrations
if(val3)
conc_count = n_conc;
figure
count_T = 1;    
for temp_T = T_start:T_class:T_end
   clear datax;
   clear datay;
   for c = 1: day_max
      datax(c) = x(count_T,conc_count,c);
      datay(c) = length_d(count_T,conc_count,c);
   end
   hold on   
   plot (datax,datay)
   text(datax(c),datay(c),num2str(temp_T))
   larvday  = hatchd(count_T,conc_count);
   larvlen  = hatchl(count_T,conc_count);
   yolkday  = yolkd(count_T,conc_count);
   yolklen  = yolkl(count_T,conc_count);
   plot (larvday,larvlen,'or')
   plot (yolkday,yolklen,'ok')
   count_T = count_T +1;   
end
title(num2str(y_conc(conc_count)))
xlabel('time [d]')
ylabel('length [mm]')
title (strcat (num2str(y_conc(conc_count)),' no/ml',' max prey size: ',num2str(prey_size)))
end

% figure dry weight against time for different temperatures and concentrations
if(val19)
conc_count = n_conc;
figure
count_T = 1;    
for temp_T = T_start:T_class:T_end
   clear datax;
   clear datay;
   for c = 1: day_max
      datax(c) = x(count_T,conc_count,c);
      datay(c) = DW_d(count_T,conc_count,c);
   end
   hold on   
   plot (datax,datay)
   text(datax(c),datay(c),num2str(temp_T))
   larvday  = hatchd(count_T,conc_count);
   larvlen  = hatchdw(count_T,conc_count);
   yolkday  = yolkd(count_T,conc_count);
   yolklen  = yolkdw(count_T,conc_count);
   plot (larvday,larvlen,'or')
   plot (yolkday,yolklen,'ok')
   count_T = count_T +1;   
end
title(num2str(y_conc(conc_count)))
xlabel('time [d]')
ylabel('dry weight [µg]')
title (strcat (num2str(y_conc(conc_count)),' no/ml',' max prey size: ',num2str(prey_size)))
end


% figure DW growth rate against fish length
if(val5)
figure
hold on
conc_count = n_conc;
count_T = 0;    
for temp_T = T_start:T_class:T_end
   count_T = count_T+1;
   clear datax;
   clear dw_growth
   n=1;
   for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)-1
      if (DW_d(count_T,conc_count,c-1)>0)
         dw_growth(n) = (DW_d(count_T,conc_count,c)/DW_d(count_T,conc_count,c-1))*100-100;
      else
         dw_growth(n) = 0;
      end
      datax(n) = length_d(count_T,conc_count,c);
      n = n+1;
   end
   if (deathday (count_T,conc_count)-1<= yolkd(count_T,conc_count)+2)
       datax(n) = 0;
       dw_growth(n) = 0;
   end
   plot (datax,dw_growth) 
   text (datax(size(datax,2)),dw_growth(size(dw_growth,2)),num2str(temp_T))
end
xlabel('fish length [mm]')
ylabel('dry weight growth [%/d]')
title (strcat (num2str(y_conc(conc_count)),' no/ml',' prey size: ',num2str(prey_size)))
end

% figure DW growth rate against concentration
if(val6)
figure
hold on
count_T = 1;    
clear datax;
clear datay;
for temp_T = T_start:T_class:T_end
   for conc_count = 1:n_conc
      clear dw_growth;       
      n=1;
      for c =  yolkd(count_T,conc_count):deathday (count_T,conc_count)
         if (DW_d(count_T,conc_count,c-1)>0)
            dw_growth(n) = (DW_d(count_T,conc_count,c)/DW_d(count_T,conc_count,c-1))*100-100;
         else
            dw_growth(n) = 0;
         end
         datay(conc_count)=mean(dw_growth);
         n = n+1;
      end
      prey_conc_class(1:prey_classes) = conc/sum(conc)*y_conc(conc_count);
      sum_dw_prey = 0;
      for c_cl = 1: prey_classes
         prey_size    = c_cl * class_size - (class_size/2); % µm
         prey_weight  = (10^(2.7722*log10(prey_size)-7.47607))*2; % 2* µg Carbon content of prey (Myrons literature research)    
         sum_dw_prey  = sum_dw_prey + prey_conc_class(c_cl)*prey_weight; %  n/ml * µg/n = µg/ml = mg/l
      end
      datax(conc_count) = sum_dw_prey*1000; % µg/l
   end
   plot (datax,datay)
   if (x_t(count_T)==8)
       plot(datax,datay,'k-o')
       savex=datax;
       savey=datay;
   end
   text (datax(n_conc),datay(n_conc),num2str(x_t(count_T)))
   xlabel('food concentration [µg DW/l]')
   ylabel('dry weight growth rate (%/d)')
   title (strcat (' max prey size: ',num2str(prey_size)))
   count_T = count_T+1;   
end
Kiorx = (1:200);
Kiory = -1.36+2.00.*log(1:200);
plot (Kiorx,Kiory,'r')
load kjoerb.txt;
plot (kjoerb(:,1), kjoerb(:,2),'x')
end

% figure DW growth rate against concentration in Carbon
if(val6)
figure
hold on
count_T = 1;    
clear datax;
clear datay;
for temp_T = T_start:T_class:T_end
   for conc_count = 1:n_conc
      clear dw_growth;       
      n=1;
      for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)
         if (DW_d(count_T,conc_count,c-1)>0)
            dw_growth(n) = (DW_d(count_T,conc_count,c)/DW_d(count_T,conc_count,c-1))*100-100;
         else
            dw_growth(n) = 0;
         end
         datay(conc_count)=mean(dw_growth);
         n = n+1;
      end
      prey_conc_class(1:prey_classes) = conc/sum(conc)*y_conc(conc_count);
      sum_dw_prey = 0;
      for c_cl = 1: prey_classes
         prey_size    = c_cl * class_size - (class_size/2); % µm
         prey_weight  = (10^(2.7722*log10(prey_size)-7.47607)); % µg Carbon content of prey (Myrons literature research)    
         sum_dw_prey  = sum_dw_prey + prey_conc_class(c_cl)*prey_weight; %no/ml*µg C
      end
      datax(conc_count) = sum_dw_prey; %µg C/ml
      datax(conc_count) = datax(conc_count) * 1000;
   end
   plot (datax,datay)
   text (datax(n_conc),datay(n_conc),num2str(x_t(count_T)))
   xlabel('food concentration [mg C /m³]')
   ylabel('dry weight growth rate (%/d)')
   title (strcat (' max prey size: ',num2str(prey_size)))
   count_T = count_T+1;   
end
end


% figure length growth rate against C concentration
if(val18)
figure
hold on
count_T = 1;    
clear datax;
clear datay;
clear ERGsx
clear ERGsy
for temp_T = T_start:T_class:T_end
    endloop = 0;
   for conc_count = 1:n_conc
      clear length_growth;       
      n=1;
      for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)
         length_growth(n) = (length_d(count_T,conc_count,c)-length_d(count_T,conc_count,c-1));
         n = n+1;
      end
      if   (deathday (count_T,conc_count)<yolkd(count_T,conc_count)+2)
          length_growth(n) = 0;
      end
      datay(conc_count)= mean(length_growth);
      prey_conc_class(1:prey_classes) = conc/sum(conc)*y_conc(conc_count); %n/l
      sum_dw_prey = 0;
      for c_cl = 1: prey_classes
         prey_size    = c_cl * class_size - (class_size/2); % µm
         prey_weight  = (10^(2.7722*log10(prey_size)-7.47607)); % µg Carbon content of prey (Myrons literature research)    
         sum_dw_prey  = sum_dw_prey + prey_conc_class(c_cl)*prey_weight; %n/ml* µg C/n = µgC / l
      end
      datax(conc_count) = sum_dw_prey; %µg C/ml = mg C/l
      datax(conc_count) = datax(conc_count) * 1000; %µgC / m³
      if (mean(length_growth)>0 & endloop == 0)
         datay2(count_T) = datax(conc_count);
         datax2(count_T) = x_t(count_T);
         endloop = 1;
      end
   end
   semilogx (datax,datay)
   ERGsx(count_T,:) = datax;
   ERGsy(count_T,:) = datay;
   count_T = count_T+1;   
end
xlabel('food concentration [mg C/m³]')
ylabel('length growth rate (mm/d)')
title (strcat (' prey size: ',num2str(prey_size)));
figure
semilogx(ERGsx',ERGsy')
text (ERGsx(:,size(ERGsx,2)),ERGsy(:,size(ERGsy,2)),num2str(x_t(:)));
clear crit
crit(:,1)  = ERGsx(1,:);
crit(:,2:size(ERGsy',2)+1) = ERGsy';

save crit_conc.dat crit -ascii

figure
semilogy(datax2,datay2,'ko')
xlabel('Temperature')
ylabel('critical prey concentration [mg C/m³]')

end

% figure GGE against length
if(val7)
figure
hold on
conc_count = n_conc;
count_T = 0;
for temp_T = T_start:T_class:T_end
   count_T = count_T+1;
   clear datax;
   clear GGE;
   n=1;
   for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)-1
      if  (food_d(count_T,conc_count,c-1) ~= 0)
         GGE(n)     = ((DW_d(count_T,conc_count,c)-DW_d(count_T,conc_count,c-1))/food_d(count_T,conc_count,c-1))*100;
      else
         GGE(n) = 0
      end
      n=n+1;
   end
   n=1;
   for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)-1
      datax(n) = length_d(count_T,conc_count,c);
      n = n+1;
   end
   if (deathday (count_T,conc_count)-1<= yolkd(count_T,conc_count)+2)
       datax(n) = 0;
       GGE(n) = 0;
   end
   plot (datax,GGE) 
   text (datax(size(datax,2)),GGE(size(GGE,2)),num2str(temp_T))
end
xlabel('fish length [mm]')
ylabel('GGE = dry weight growth / consumed food *100 [%]')
title (strcat (num2str(y_conc(conc_count)),' no/ml',' prey size: ',num2str(prey_size)))
end

% figure ingestion and dry weight growth
if(val8)
figure
hold on
conc_count = n_conc;
clear datax;
clear datay;
n=1;
%for conc_count = 1: n_conc
count_T = 0;   
    
for temp_T = T_start:T_class:T_end
   count_T = count_T+1;
   for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)-1
      datay(n) = DW_d(count_T,conc_count,c)-DW_d(count_T,conc_count,c-1);
      datax(n) = food_d(count_T,conc_count,c-1);
      n=n+1;
   end
   if (deathday (count_T,conc_count)<= yolkd(count_T,conc_count))
       datax(n:n+1) = 0;
       datay(n:n+1) = 0;
   end
   plot (datax,datay,'+') 
end
%end

xlabel('ingested food µg/d')
ylabel('dry weight growth µg/d')
title (strcat (num2str(y_conc(conc_count)),' no/ml',' prey size: ',num2str(prey_size)))
Kiorboex = (10:100);
Kiorboey1 = (0.17.*(10:100)-2.5); % paper Environmental biology of fishes 1986
Kiorboey2 = (0.30.*(10:100)-5.85);% paper Environmental biology of fishes 1986
Kiorboey3 = (0.40.*(10:100)-7.26);% paper Environmental biology of fishes 1986
plot(Kiorboex,Kiorboey1)
plot(Kiorboex,Kiorboey2)
plot(Kiorboex,Kiorboey3)
axis ([0 150 0 150])
end

% figure ingested DW per day against fish length
if(val9)
figure
hold on
conc_count = n_conc;
count_T = 0;    
for temp_T = T_start:T_class:T_end
   count_T = count_T+1;
   c=2;
   clear datax;
   clear datay;
   n=1;
   for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)-1
      datax(n) = length_d(count_T,conc_count,c);
      datay(n) = food_d(count_T,conc_count,c);
      n=n+1;
   end
   if (deathday (count_T,conc_count)-1<= yolkd(count_T,conc_count)+2)
       datax(n:n+1) = 0;
       datay(n:n+1) = 0;
   end
   plot (datax,datay) 
   text (datax(size(datax,2)),datay(size(datay,2)),num2str(temp_T))
end
xlabel('fish length [mm]')
ylabel('ingested dry weight[µg/d]')
title (strcat (num2str(y_conc(conc_count)),' no/ml',' prey size: ',num2str(prey_size)))
end

% figure ingested DW per day against fish DW
if(val10)
figure
hold on
conc_count = n_conc;
count_T = 0;    
for temp_T = T_start:T_class:T_end
   count_T = count_T+1;
   c=2;
   clear datax;
   clear datay;
   n=1;
   for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)-1
      datax(n) = DW_d(count_T,conc_count,c);
      datay(n) = food_d(count_T,conc_count,c);
      n=n+1;
   end
   if (deathday (count_T,conc_count)-1<= yolkd(count_T,conc_count)+2)
       datax(n:n+1) = 0;
       datay(n:n+1) = 0;
   end
   plot (datax,datay) 
   text (datax(size(datax,2)),datay(size(datay,2)),num2str(temp_T))
end
xlabel('fish dry weight[µg]')
ylabel('ingested dry weight[µg/d]')
title (strcat (num2str(y_conc(conc_count)),' no/ml',' prey size: ',num2str(prey_size)))
end

% figure ingested particles per day against fish length
if(val11)
figure
hold on
conc_count = n_conc;
count_T = 0;    
for temp_T = T_start:T_class:T_end
   count_T = count_T+1;
   n=1;
   clear datax;
   clear datay;
   for c = yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)-1
      datax(n) = length_d(count_T,conc_count,c);
      datay(n)= part_d(count_T,conc_count,c);
      n=n+1;
   end
   if (deathday (count_T,conc_count)-1<= yolkd(count_T,conc_count)+2)
       datax(n:n+1) = 0;
       datay(n:n+1) = 0;
   end
   plot (datax,datay) 
   text (datax(size(datax,2)),datay(size(datay,2)),num2str(temp_T))
end
xlabel('fish length [mm]')
ylabel('ingested particles [n/d]')
title (strcat (num2str(y_conc(conc_count)),'no/ml',' prey size: ',num2str(prey_size)))
end

% plot feeding strike frequency according to length and prey concentration
if(val12)
clear datax
clear datay
clear FSF
c_L=1;
for L = 10:35
   c_conc = 1; 
   FSF_L = (0.0237*log(L)-0.0484);
   if (FSF_L<0)FSF_L  = 0.0002;end
   datax(c_L) = L;
   for conc = 0.01:0.01:0.2
      FSF_conc           = (0.0216*(conc)^0.3296);
      FSF(c_L,c_conc)    = FSF_L * FSF_conc*1000;
      datay(c_conc) = conc;
      c_conc = c_conc +1;
   end
   c_L = c_L+1;
end
figure
contourf(datay,datax,FSF)
title ('feeding strike frequency according to length and prey concentration')
xlabel ('prey concentration [n/ml]')
ylabel ('larval length [mm]')
colorbar
end

% plot ranking function CS HT and prey weight
if(val13)
figure 
hold on
l=1;
clear datax
clear datay
clear rk
for cL = 10:35 
   for c_cl = 1: prey_classes
      PS_max = max_p_a / (1+(cL/14)^(-max_p_b));
      prey_size2(c_cl)   = c_cl * class_size - (class_size/2); %prey size in the middle of the class
      prey_weight(c_cl)  = (10^(2.7722*log10(prey_size2(c_cl))-7.47607))*2; % 2* µg Carbon content of prey (Myrons literature research)
      CS(c_cl)           = 1.1-(CS_b*prey_size2(c_cl)/PS_max);              % capture success 
      if (CS(c_cl)<0) CS(c_cl)=0; end
      HT(c_cl)     = exp(HT_a*10^(HT_b*((prey_size2(c_cl)/1000)/cL)));   % handling time [s]
      rk(c_cl)     = (prey_weight(c_cl)*CS(c_cl))/HT(c_cl);                      % prey type ranking
   end
   % sort classes from lowest to highest priority
   clear rkind;
   clear rkind2;
   for c_cl = 1: prey_classes 
      rk_count = 1;
      for c_cl2 = 1: prey_classes
         if (rk (c_cl)<rk (c_cl2))
            rk_count = rk_count +1;
         end
      end
      rkind2(c_cl) = rk_count;
   end
   for c_cl = 1: prey_classes
      rkind(c_cl) = 999;                
   end
   for c_cl = 1: prey_classes
      if (rk(c_cl)>0)
         rkind(rkind2(c_cl)) = c_cl;                
      end
   end
   datay1(l)=prey_size2(rkind(1));
   datay2(l)=prey_size2(rkind(2));
   datay3(l)=prey_size2(rkind(3));
   if (rkind(4) ~= 999)
   datay4(l)=prey_size2(rkind(4));
   end
   if (rkind(5) ~= 999)
   datay5(l)=prey_size2(rkind(5));
   end
   if (rkind(6) ~= 999)
   datay6(l)=prey_size2(rkind(6));
   end
   if (rkind(7) ~= 999)
   datay7(l)=prey_size2(rkind(7));
   end
   if (rkind(8) ~= 999)
   datay8(l)=prey_size2(rkind(8));
   end
   if (rkind(9) ~= 999)
   datay9(l)=prey_size2(rkind(9));
   end
   datax(l)=cL;
   l=l+1;
end
plot (datax,datay1,'r')    
text (datax(size(datax,2)),datay1(size(datay1,2)),'rank 1')
plot (datax,datay2,'b')    
text (datax(size(datax,2)),datay2(size(datay2,2)),'rank 2')
plot (datax,datay3,'k')    
text (datax(size(datax,2)),datay3(size(datay3,2)),'rank 3')
plot (datax,datay4,'g')    
text (datax(size(datax,2)),datay4(size(datay4,2)),'rank 4')
if (rkind(5) ~= 999)
plot (datax,datay5,'k')    
text (datax(size(datax,2)),datay5(size(datay5,2)),'rank 5')
end
if (rkind(6) ~= 999)
plot (datax,datay6,'c')    
text (datax(size(datax,2)),datay6(size(datay6,2)),'rank 6')
end
if (rkind(7) ~= 999)
plot (datax,datay7,'y')    
text (datax(size(datax,2)),datay7(size(datay7,2)),'rank 7')
end
if (rkind(8) ~= 999)
plot (datax,datay8,'y')    
text (datax(size(datax,2)),datay8(size(datay8,2)),'rank 8')
end
if (rkind(9) ~= 999)
plot (datax,datay9,'y')    
text (datax(size(datax,2)),datay9(size(datay9,2)),'rank 9')
end
title ('Best prey size according to ranking')
xlabel ('larval length [mm]')
ylabel ('prey size [µm]')
x_myron = [12 14 15 16 17 20];
y_myron_mean = [286 613 685 690 722 663];
plot(x_myron,y_myron_mean,'ko')
y_myron_max = [325 768 820 877 879 889];
plot(x_myron,y_myron_max,'k+')
load PL_HL.txt
x_Susan = PL_HL(:,1);
y_Susan = PL_HL(:,2)*1000;
plot(x_Susan,y_Susan,'ro')

%plot((10:35), (max_p_a ./ (1+((10:35)./14).^(-max_p_b))))

%plot((10:35),(10:35).*2.384/100*1000) 
%%%%%%%%%%%%%%%%%%%%%
figure
hold on
plot (datax,datay1,'r')    
text (datax(size(datax,2)),datay1(size(datay1,2)),'rank 1')
title ('Best prey size according to ranking')
xlabel ('larval length [mm]')
ylabel ('prey size [µm]')
x_myron = [12 14 15 16 17 20];
y_myron_mean = [286 613 685 690 722 663];
plot(x_myron,y_myron_mean,'ko')
y_myron_max = [325 768 820 877 879 889];
plot(x_myron,y_myron_max,'k+')
load PL_HL.txt
x_Susan = PL_HL(:,1);
y_Susan = PL_HL(:,2)*1000;
plot(x_Susan,y_Susan,'ro')
plot((10:35), (max_p_a ./ (1+((10:35)./14).^(-max_p_b))))
plot((10:35),(10:35).*2.384/100*1000) 

end

% figure DW at length
if(val14)
figure
hold on
conc_count = n_conc;
count_T = 0;    
for temp_T = T_start:T_class:T_end
   count_T = count_T+1;
   clear datax;
   clear datay;
   n=1;
   for c = yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)-1
      datax(n) = length_d(count_T,conc_count,c);
      datay(n) = DW_d(count_T,conc_count,c);
      n=n+1;
   end
   if (deathday (count_T,conc_count)-1<= yolkd(count_T,conc_count)+2)
       datax(n:n+1) = 0;
       datay(n:n+1) = 0;
   end
   plot (datax,datay) 
   text (datax(size(datax,2)),datay(size(datay,2)),num2str(temp_T))
end
xlabel('fish length [mm]')
ylabel('fish dry weight[µg]')
plot(10:20, 0.018521.*(10:20).^(3.613981*exp(7*0.006266)),'r')
text (20, 0.018521.*(20)^(3.613981*exp(7*0.006266)),'7°C laboratory')
plot(10:20, 0.018521.*(10:20).^(3.613981*exp(10*0.006266)),'r')
text (20, 0.018521.*(20).^(3.613981*exp(10*0.006266)),'10°C laboratory')
plot(10:20, 0.018521.*(10:20).^(3.613981*exp(13*0.006266)),'r')
text (20, 0.018521.*(20).^(3.613981*exp(13*0.006266)),'13°C laboratory')
end

%temperature history
if(val15)
figure
hold on
for c=1:size(T_hist,1)
   plot(T_hist(c,:))
   xlabel('time')
   ylabel('temperature')
end
end

if(val16)
figure
hold on
conc_count = n_conc;
count_T = 0;    
for temp_T = T_start:T_class:T_end
    count_T = count_T+1;
    for conc_count = 1:n_conc
      clear datax;
      clear datay;
      for c = 1: size(x,3)
         datax(c) = length_d(count_T,conc_count,c);
         datay(c) = DW_d(count_T,conc_count,c);
      end
      if (temp_T == 10)
         plot (datax,datay,'o') 
         text (datax(size(datax,2)),datay(size(datay,2)),strcat(num2str(y_conc(conc_count)),' ',num2str(temp_T)))
      end
   end
end
xlabel('fish length [mm]')
ylabel('fish dry weight[µg]')
plot(10:20, 0.018521.*(10:20).^(3.613981*exp(10*0.006266)),'r')
text (20, 0.018521.*(20).^(3.613981*exp(10*0.006266)),'10°C laboratory')
end

% figure dry weight against degree days for starvation plot
if(val21)
clear datax
clear datay
conc_count = n_conc;
figure
count_T = 1;    
for temp_T = T_start:T_class:T_end
   clear datax;
   clear datay;
   for c = (floor(59/temp_T)+1): day_max
      datax(c) = x(count_T,conc_count,c)*temp_T;
      datay(c) = DW_d(count_T,conc_count,c);
   end
   hold on   
   plot (datax,datay,'o')
   count_T = count_T +1;   
end
title(num2str(y_conc(conc_count)))
xlabel('degree-days-post-hatch [°d]')
ylabel('dry weight [µg]')
title ('starvation')
plot (0:160,89.1-0.33.*(0:160))
end

% Volume searched per day according to length 
if(val22)
figure
hold on
count_T = 1;    
for temp_T = T_start:T_class:T_end
   for conc_count = 1:n_conc
      clear datax;
      clear datay;
      n = 1;
      for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)-1
         datax(n) = length_d(count_T,conc_count,c);
         datay(n) = vol_searched_d(count_T,conc_count,c);
         n = n+1;
      end
      plot (datax, datay);
      if (temp_T == T_start)
          text (datax(:),datay(:),num2str(y_conc(conc_count))); 
      end
   end  
end
xlabel('larval length [mm]')
ylabel('Volume searched per day (l/d)')
end

% particles in the gut
if(val23)
figure
%plot (length,part_in_gut)
%semilogy (length,part_in_gut,'o')
plot(length,part_in_gut,'o')
axis ([5 25 1 100])
xlabel('larval length [mm]')
ylabel('particles in the gut ')
load part_gut.txt
hold on
%semilogy (part_gut(:,1),part_gut(:,2),'x') 
plot (part_gut(:,1),part_gut(:,2),'x') 
end

% Encounters per day  
if(val24)
figure
hold on
count_T = 1;    
for temp_T = T_start:T_class:T_end
   for conc_count = 1:n_conc
      clear datax;
      clear datay;
      n = 1;
      for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)
         datax(n) = length_d(count_T,conc_count,c);
         datay(n) = pos_cons_d(count_T,conc_count,c);
         n = n+1;
      end
      plot (datax, datay);
      text (datax(:,size(datax,2)),datay(:,size(datay,2)),num2str(temp_T));
   end
end
xlabel('larval length [mm]')
ylabel('number of particles encountered per day')
end

% reactive distance
if(val25)
figure
hold on
for c_L = 10:25
   clear datax
   clear datay
   n=1;
   for c_cl = 1: prey_classes
      alpha        = 0.0167*exp(9.14-2.4*log(c_L)+0.229*(log(c_L))^2); % angle of visual acuity
      alpha        = alpha * pi/180;
      prey_size    = c_cl * class_size - (class_size/2);       %prey size in the middle of the class µm
      RD(c_cl)     = (prey_size/1000)/(2*tan(alpha/2));        % reactive distance [mm]
      datax (n) = prey_size;
      datay (n) =  RD(c_cl);
      n=n+1;
   end   
   plot (datax, datay)
   text (prey_size,RD(c_cl),num2str(c_L))
end
xlabel('prey size [µm]')
ylabel('Reactive distance [mm]')   
end
   
      
% percentage body weight ingested per day 
if(val26)
figure
hold on
conc_count = n_conc;
count_T = 0;    
for temp_T = T_start:T_class:T_end
   count_T = count_T+1;
   c=2;
   clear datax;
   clear datay;
   n=1;
   for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)-1
      datax(n) = length_d(count_T,conc_count,c);
      datay(n) = food_d(count_T,conc_count,c)/DW_d(count_T,conc_count,c)*100;
      n=n+1;
   end
   if (deathday (count_T,conc_count)-1<= yolkd(count_T,conc_count)+2)
       datax(n:n+1) = 0;
       datay(n:n+1) = 0;
   end
   plot (datax,datay) 
   text (datax(size(datax,2)),datay(size(datay,2)),num2str(temp_T))
end
xlabel('fish length [mm]')
ylabel('ingested dry weight as % body weight [%/d]')
end

% figure length against time for different temperatures and concentrations
% together with literature data
if(val27)
conc_count = n_conc;
figure
hold on
count_T = 1;
for temp_T = T_start:T_class:T_end
    if temp_T == 11
        count_T2=count_T;
    end
    count_T = count_T +1;;    
end
count_T = count_T2;
for conc_count = 1:n_conc
   clear datax;
   clear datay;
   for c = 1: day_max
      datax(c) = x(count_T,conc_count,c);
      datay(c) = length_d(count_T,conc_count,c);
   end   
   plot (datax,datay)
   text(datax(c),datay(c),num2str(y_conc(conc_count)))
   larvday  = hatchd(count_T,conc_count);
   larvlen  = hatchl(count_T,conc_count);
   yolkday  = yolkd(count_T,conc_count);
   yolklen  = yolkl(count_T,conc_count);
   plot (larvday,larvlen,'or')
   plot (yolkday,yolklen,'ok') 
   datax'
   datay'
end
title(num2str(x_t(count_T)))
xlabel('time [d]')
ylabel('length [mm]')

load conc_gr.txt
temp(1,1) = 0;
temp(2,1) = yolkday;
for c = 3: size(conc_gr,1)
   temp(c,1) = temp(c-1,1)+ conc_gr(c,1)-conc_gr(c-1,1);
end
conc_gr(:,1)=temp(:,1);
for c = 2: size(conc_gr,2)
   plot(conc_gr(2:size(conc_gr,1),1),conc_gr(2:size(conc_gr,1),c),'r')
   text(conc_gr(size(conc_gr,1),1),conc_gr(size(conc_gr,1),c),num2str(conc_gr(1,c)))
end
%Geffen
plot([7 14 22],[8.49 11.29 11.82],'go')
plot([7 14 22],[10.06 11.13 13.17],'g+')

%Pedersen
plot([7	9.0109	13.5574	18.5847	23.5246	28.6831	33.6667	38.6066],[7.6552	8.1049	9.3683	10.5032	12.1306	13.3298	13.9293	14.4433],'bo')
plot([7	9.0109	13.5574	18.5847	23.5246	28.6831	33.6667	38.6066],[7.6552	8.1049	9.1328	9.8822	11.6381	12.5375	13.0086	13.1585],'bx')
end


% figure length against time for different temperatures and concentrations
% together with literature data
if(val28)
conc_count = n_conc;
figure
hold on
count_T = 1;
for temp_T = T_start:T_class:T_end
    if temp_T == 8
        count_T2=count_T;
    end
    count_T = count_T +1;;    
end
count_T = count_T2;
for conc_count = 1:n_conc
   clear datax;
   clear datay;
   for c = 1: day_max
      datax(c) = x(count_T,conc_count,c);
      datay(c) = DW_d(count_T,conc_count,c);
   end   
   plot (datax,datay)
   text(datax(c),datay(c),num2str(y_conc(conc_count)))
   larvday  = hatchd(count_T,conc_count);
   larvlen  = hatchdw(count_T,conc_count);
   yolkday  = yolkd(count_T,conc_count);
   yolklen  = yolkdw(count_T,conc_count);
   plot (larvday,larvlen,'or')
   plot (yolkday,yolklen,'ok') 
end
title(num2str(x_t(count_T)))
xlabel('time [d]')
ylabel('dry weight [µg]')

clear temp
load conc_gr2.txt
temp(1,1) = 0;
temp(2,1) = yolkday;
for c = 3: size(conc_gr2,1)
   temp(c,1) = temp(c-1,1)+ conc_gr2(c,1)-conc_gr2(c-1,1);
end

conc_gr2(:,1)=temp(:,1);
for c = 2: size(conc_gr2,2)
   plot(conc_gr2(2:size(conc_gr2,1),1),conc_gr2(2:size(conc_gr2,1),c),'r')
   text(conc_gr2(size(conc_gr2,1),1),conc_gr2(size(conc_gr2,1),c),num2str(conc_gr2(1,c)))
end
end


% comparison of food ingested vs. food concentration with Mac Kenzie
% Literature review
if(val29)
figure
n=1;
clear datax;
clear datay;
for conc_count = 1:n_conc
   prey_conc_class(1:prey_classes) = conc/sum(conc)*y_conc(conc_count);
   sum_dw_prey = 0;
   for c_cl = 1: prey_classes
      prey_size    = c_cl * class_size - (class_size/2); % µm
      prey_weight  = (10^(2.7722*log10(prey_size)-7.47607))*2; % 2* µg Carbon content of prey (Myrons literature research)    
      sum_dw_prey  = sum_dw_prey + prey_conc_class(c_cl)*prey_weight; %  n/ml * µg/n = µg/ml = mg/l
   end
   datax_m(conc_count)= sum_dw_prey*1000; % µg/l;
   count_T = 0; 
   for temp_T = T_start:T_class:T_end
      count_T = count_T+1;
      for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)-1
         datax(n) = sum_dw_prey*1000;
         datay(n) = food_d(count_T,conc_count,c)/DW_d(count_T,conc_count,c);
%         datay(n) = food_d(count_T,conc_count,c);
         n=n+1;
      end
      if (deathday (count_T,conc_count)-1<= yolkd(count_T,conc_count)+2)
         datay(n:n+1) = 0;
      end
   end
end
loglog (datax,datay,'o') 
xlabel('food concentration [µg/l]')
ylabel('ingested dry weight[µg/d]')
hold on
load MacKenz.txt
loglog (MacKenz(:,1),MacKenz(:,2),'x') 

end

% clearance rate particles ingested/encounters
if(val30)
figure
hold on
count_T = 1;    
datax(1) = 1;
datay(1) = 1;

for temp_T = T_start:T_class:T_end
   for conc_count = 1:n_conc
      clear datax;
      clear datay;
      n = 1;
      for c =  yolkd(count_T,conc_count)+2:deathday (count_T,conc_count)
         if ( pos_cons_d(count_T,conc_count,c)>0)
            datax(n) = length_d(count_T,conc_count,c);
            datay(n) = part_d(count_T,conc_count,c)/pos_cons_d(count_T,conc_count,c);
            n = n+1;
         end
      end
      plot (datax, datay);
   end
   count_T = count_T+1;
end
xlabel('larval length [mm]')
ylabel('n ingested / n encountered per day')
end


% growth per day vs. ingestion per day 

if(val31)
figure
hold on
count_T = 1;   
clear datax;
clear datay;  
datax(1) = 1;
datay(1) = 1;
n =1;
for conc_count = 1:n_conc
   for c =  yolkd(count_T,conc_count):deathday (count_T,conc_count)
      if (DW_d(count_T,conc_count,c)>600)
         datax(n) = food_d(count_T,conc_count,c)/DW_d(count_T,conc_count,c-1);
         datay(n) = (DW_d(count_T,conc_count,c)-DW_d(count_T,conc_count,c-1))/DW_d(count_T,conc_count,c-1);
         n=n+1;
      end
   end
   if (n>1)
      plot ((datax),(datay),'+') 
   end
end
plot((0.01:0.01:0.5), 0.12.*(1-exp(-7.6.*(0.01:0.01:0.5)))-0.029) % Kiorb (Checkley)

end

if(val32)
for conc_count = 1:n_conc
   figure
   hold on
   title (strcat (num2str(y_conc(conc_count)),' no/ml'))
   xlabel('age')
   ylabel('growth mm/d')
   count_T = 1;    
   for temp_T = T_start:T_class:T_end
      clear datax; 
      clear datay; 
      n = 1;
      for c =  (yolkd(count_T,conc_count))+2:deathday (count_T,conc_count)-1
         datax(n) = c - yolkd(count_T,conc_count)+2;
         datay(n) = (length_d(count_T,conc_count,c)-length_d(count_T,conc_count,c-1));
         n = n+1;
      end
      if (deathday (count_T,conc_count)-1<=yolkd(count_T,conc_count)+2)
         datax(n) = 0;
         datay(n) = 0;
      end
      plot (datax, datay)
      text (datax(size(datax,2)),datay(size(datay,2)),num2str(temp_T))
      count_T = count_T+1;
   end
end
end

% figure length against time for different temperatures and concentrations
if(val33)
Bolz_and_Burns_age=[31.1	52	131.1	148.4	27.3	44.5	63.2	31.5	50.5	74.8	30.9	48.8	60.7	89.1	27.1	53.9	65.5	81.2	31	59	80.5	31.7	46.9	68.6];
Bolz_and_Burns_length=[16.9	21.5	30.3	31.8	14.4	21.8	25.4	17.3	23.4	28	15.9	22.2	24.3	29.7	13.5	24.8	25.8	28.6	14.1	24.9	27.7	17	22.1	27.7];
Bolz_and_Burns_age = Bolz_and_Burns_age;
figure
for conc_count = 1:n_conc
   %figure
   hold on
   %title (strcat (num2str(y_conc(conc_count)),' no/ml'))
   xlabel('age')
   ylabel('length [mm]')
   count_T = 1;    
   for temp_T = T_start:T_class:T_end
      clear datax; 
      clear datay; 
      n = 1;
      for c =  (yolkd(count_T,conc_count))+2:deathday (count_T,conc_count)-1
         datax(n) = x(count_T,conc_count,c);
         datay(n) = length_d(count_T,conc_count,c);
         n = n+1;
      end
      if (deathday (count_T,conc_count)-1<=yolkd(count_T,conc_count)+2)
         datax(n) = 0;
         datay(n) = 0;
      end
      plot (datax, datay)
      plot(Bolz_and_Burns_age,Bolz_and_Burns_length,'k+')
      text (datax(size(datax,2)),datay(size(datay,2)),num2str(y_conc(conc_count)))
      count_T = count_T+1;
      
   end
end
end


% figure length against time for different temperatures and concentrations
if(val34)
Henderson_age1 = [0	1	6	7	8	13	14	15	17	18	19	20	23	24	26	28	35	42];
Henderson_age2 = [0	1	3	4	5	6	8	11	12	14	15	16	21	22	30	35	48	49];
Henderson_L1   = [7.0708	7.4788	8.6006	9.0595	9.0595	10.2323	9.8754	10.2323	10.5382	11.3031	9.9773	11.2011	11.2011	10.6402	11.864	14.7195	16.6062	19.0538];
Henderson_L2   = [6.5099	6.9688	7.4788	7.9887	8.1416	7.8357	8.9065	8.4476	7.9377	8.4986	8.4476	9.3144	9.7224	10.3343	12.5269	14.4136	21.5014	23.0822];

for conc_count = 1:n_conc
   figure
   hold on
   title (strcat (num2str(y_conc(conc_count)),' no/ml'))
   xlabel('age')
   ylabel('length [mm]')
   count_T = 1;    
   for temp_T = T_start:T_class:T_end
      clear datax; 
      clear datay; 
      n = 1;
      for c =  (yolkd(count_T,conc_count))+2:deathday (count_T,conc_count)-1
         datax(n) = x(count_T,conc_count,c);
         datay(n) = length_d(count_T,conc_count,c);
         n = n+1;
      end
      if (deathday (count_T,conc_count)-1<=yolkd(count_T,conc_count)+2)
         datax(n) = 0;
         datay(n) = 0;
      end
      plot (datax, datay)
      plot(Henderson_age1,Henderson_L1,'k+')
%      plot(Henderson_age2,Henderson_L2,'k+')
      text (datax(size(datax,2)),datay(size(datay,2)),num2str(temp_T))
      count_T = count_T+1;
      
   end
end
end


% Dry weight condition index (dry weight at length)
if(val35)
figure
hold on
for conc_count = 1:n_conc
   xlabel('length [mm]')
   ylabel('dryweight / length')
   count_T = 1;    
   for temp_T = T_start:T_class:T_end
      clear datax; 
      clear datay; 
      n = 1;
      for c =  (yolkd(count_T,conc_count))+2:deathday (count_T,conc_count)-1
         datax(n) = length_d(count_T,conc_count,c);
         datay(n) = (DW_d(count_T,conc_count,c)/1000)/length_d(count_T,conc_count,c)^3*100;
         n = n+1;
      end
      if (deathday (count_T,conc_count)-1<=yolkd(count_T,conc_count)+2)
         datax(n) = 0;
         datay(n) = 0;
      end
      plot (datax, datay,'.k')
%      text (datax(size(datax,2)),datay(size(datay,2)),num2str(temp_T))
      count_T = count_T+1;
   end
end
plot (10:0.01:15, 0.0074, '-r')  % Checkley 1984
plot (10:0.01:15, 0.0165, '-r')  % Checkley 1984
plot (10:0.01:20, 0.0005.*(10:0.01:20)+0.0098,'-g') % Westernhagen & Rosenthal 1979
axis ([10 20 0 0.03])
end






