 function [Edata EW_bj EW_je EN_em] = predict_Venturia_canescens(par, data, tW_bj, tW_je, t_em)
% par: (27)-vector with parameters (see below)
% data:(r_d,1)-matrix with zero-variate data (not all elements are used)
%     notice that nmregr requires column 2 and 3 as well
% tW_bj,.,t_em: data matirces
% Edata:(r_d,1)-matrix with predicted values for data(:,2)
% EW_bj - EN_em: predicted values for tW_bj(:,2),., tN(:,2)

  global dwm % import values set in mydata_Venturia_canescens

%% unpack par
  T_ref  = par(1); % K, temp for which rate pars are given 
  T_A    = par(2); % K, Arrhenius temp
  f      = par(3); % -, scaled functional response
  z      = par(4); % -, zoom factor
  del_M  = par(5); % -, shape coefficient to convert vol-length to physical length
  F_m    = par(6); % l/d.cm^2, {F_m}, max spec searching rate
  kap_X  = par(7); % -, digestion efficiency of food to reserve
  v      = par(8); % cm/d, energy conductance
  kap    = par(9); % -, allocation fraction to soma = growth + somatic maintenance
  kap_R  = par(10);% -, reproduction efficiency
  p_M    = par(11);% J/d.cm^3, [p_M], vol-specific somatic maintenance
  p_T    = par(12);% J/d.cm^2, {p_T}, surface-specific som maintenance
  k_J    = par(13);% 1/d, maturity maint rate coefficient
  E_G    = par(14);% J/cm^3, [E_G], spec cost for structure
  E_Hx   = par(15);% J, E_H^b, maturity at pre-birth (when embryo starts feeding)
  E_Hb   = par(16);% J, E_H^b, maturity at birth
  y_EV_l = par(17);% mol/mol, yield of imago reserve on larval structure
  E_He   = par(18);% J, E_H^p, maturity at emergence
  E_H0   = par(19);% J, maturity at egg maturation
  h_a    = par(20);% 1/d^2, Weibull aging acceleration
  s_G    = par(21);% -, Gompertz stress coefficient
  s_j    = par(22);% -, reprod buffer/structure at pupation as fraction of max
  s_s    = par(23);% -, reprod buffer/structure at silk production, as fraction of max [ER]
  kap_SA = par(24);% -, weight coefficient for assimilation
  kap_SG = par(25);% -, weight coefficient for growth
  kap_SD = par(26);% -, weight coefficient for dissipation
  k_E_im = par(27);% 1/d, reserve turnover and reproduction buffer turnover in the imago
  
  d_Vx = dwm(2,4);  d_Ex = dwm(3,4);% g/cm^3, specific densities for structure and reserve for the embryo
  d_V = dwm(2,1);   d_E = dwm(3,1); % g/cm^3, specific densities for structure and reserve
  w_V = dwm(2,2);   w_E = dwm(3,2); % g/mol,  molecular weights for structure and reserve
  mu_V = dwm(2,3);  mu_E = dwm(3,3);% J/mol,  chemical potentials for structure and reserve

  % Selected copy-paste from parscomp & statistics
  p_Am = z * p_M/ kap;             % J/d.cm^2, {p_Am} max spec assimilation flux
  J_E_Am = p_Am/ mu_E;             % mol/d.cm^2, {J_EAm}, max surface-spec assimilation flux
  k_M = p_M/ E_G;                  % 1/d, somatic maintenance rate coefficient
  k = k_J/ k_M;                    % -, maintenance ratio
  % p_Xm = p_Am/ kap_X;            % J/d.cm^2, max spec feeding power

  M_V = d_V/ w_V;                  % mol/cm^3, [M_V], volume-specific mass of structure
  kap_G = M_V * mu_V/ E_G;         % -, growth efficiency
  y_V_E = mu_E * M_V/ E_G;         % mol/mol, yield of structure on reserve
  y_E_V = 1/ y_V_E;                % mol/mol, yield of reserve on structure

  E_m = p_Am/ v;                   % J/cm^3, [E_m] reserve capacity 
  m_Em = y_E_V * E_m/ E_G;         % mol/mol, reserve capacity 
  g = E_G/ kap/ E_m;               % -, energy investment ratio
  w = m_Em * w_E/ w_V;             % -, contribution of reserve to weight

  L_m = v/ k_M/ g;                 % cm, maximum structural length
  L_T = p_T/ p_M;                  % cm, heating length (also applies to osmotic work)
  l_T = L_T/ L_m;     
  
  
  %maturity at pre-birth
  M_Hx = E_Hx/ mu_E;               % mol, maturity at pre-birth  
  U_Hx = M_Hx/ J_E_Am;             % cm^2 d, scaled maturity at pre-birth
  u_Hx = U_Hx * g^2 * k_M^3/ v^2;  % -, scaled maturity at pre-birth
  V_Hx = U_Hx/ (1 - kap);          % cm^2 d, scaled maturity at pre-birth
  v_Hx = V_Hx * g^2 * k_M^3/ v^2;  % -, scaled maturity at pre-birth
  % maturity at birth
  M_Hb = E_Hb/ mu_E;               % mol, maturity at birth  
  U_Hb = M_Hb/ J_E_Am;             % cm^2 d, scaled maturity at birth
  u_Hb = U_Hb * g^2 * k_M^3/ v^2;  % -, scaled maturity at birth
  V_Hb = U_Hb/ (1 - kap);          % cm^2 d, scaled maturity at birth
  v_Hb = V_Hb * g^2 * k_M^3/ v^2;  % -, scaled maturity at birth
  % maturity at pupation
  M_He = E_He/ mu_E;               % mol, maturity at emergence
  U_He = M_He/ J_E_Am;             % cm^2 d, scaled maturity at emergence 
  u_He = U_He * g^2 * k_M^3/ v^2;  % -, scaled maturity at emergence 
  V_He = U_He/ (1 - kap);          % cm^2 d, scaled maturity at emergence
  v_He = V_He * g^2 * k_M^3/ v^2;  % -, scaled maturity at emergence
 
  %% zero-variate data
  TC = tempcorr(data(1,1), T_ref, T_A);   % -, Temperature Correction factor for a_b; ONLY WORK WITH T_A
  %kT_J = k_J * TC;       % 1/d, temp corrected maintenance rate coeff
  
  % initial
  pars_UE0 = [V_Hx; g; k_J; k_M; v];
  [U_E0, L_x, info] = initial_scaled_reserve(f, pars_UE0);
  if info ~= 1 % numerical procedure failed
    fprintf('warning: invalid parameter value combination for initial scaled reserve \n')
  end
  E_0 = mu_E * J_E_Am * U_E0;       % J, initial reserve
  Wd_0 = w_E * J_E_Am * U_E0;       % g, initial dry weight
 
  %pre-birth
  pars_tb = [g; k; l_T; v_Hx; v_Hb]; % compose parameter vector
  %pars_lx = [g; k; v_Hx];           % compose parameter vector
  [t_b, t_x, l_b, l_x, info] = get_tp(pars_tb, f); % -, scaled length at birth at f
  %[t_x, l_x, info] = get_tb(pars_lx, f, L_x/ L_m); % -, scaled and length at birth at f
  if info ~= 1 % numerical procedure failed
    fprintf('warning: invalid parameter value combination for get_tb \n')
  end
  L_x = L_m * l_x;                   % cm, structural length at pre-birth at f
  L_b = L_m* l_b;                    % cm, structural length at birth at f
 % Lw_x = L_x/ del_M;                % cm, physical length at pre-birth at f
  Wd_x = L_x^3 * d_Vx * (1 + f * w); % g, dry weight at pre-birth at f (remove d_V for wet weight)
  aT_x = t_x/ k_M/TC;                % d, age at pre-birth at f and T
  Wd_b = L_b^3 * d_V * (1 + f * w);  % g, dry weight at pre-birth at f (remove d_V for wet weight)
  aT_b = t_b/k_M/TC;                 % d, age at pre-birth at f and T
  k_E = v/ L_b;        % 1/d, reserve turnover and spec decay of larval structure on pupa
  p_Amb = p_Am/ L_b;                % J/d.cm^3, spec max assimilation at birth
 
     
   
  % from birth to pupation (from b to j,metamorphosis)
  p_J = k_J * E_Hb;                  % J/d, maturity maint rate during larval stage
  
  %r = (f * k_E - g * k_M)/ (f + g); % 1/d, specific growth rate
  E_R_ref = (1 - kap) * E_m * g * (k_E + k_M)/ (k_E - g * k_M); % J/cm^3, reference value for [E_R]
  E_Rj = s_j * E_R_ref; % J/cm^3, threshold reprod buffer density at pupation
  E_Rs = s_s * E_R_ref; % J/cm^3, threshold reprod buffer density at silk production
  aELRS_0 = [0; f * E_m * L_b^3; L_b; 0;  0]; % initial state vector [time ,E ,L ,E_R ,E_Silk]
  [ER aELRS] = ode45(@dget_aELRS, [0; E_Rj], aELRS_0, [], f, k_M, k_E , p_J, p_Amb, E_m, g, kap, kap_SA, kap_SG, kap_SD, p_M, E_Rs);

  t_j = aELRS(end,1); aT_j = aT_b + t_j/ TC;    % d, age at pupation 
  E_j = aELRS(end,2);                % J, reserve just before pupation
  L_j = aELRS(end,3);                % cm, structural length at pupation
  E_Rj = aELRS(end,4);                % J, reprod buffer just before pupation
  E_Sj = aELRS(end,5);                % J, silk as a product at pupation
  V_j = L_j^3;                      % cm3; structure of larva at the end larval period
  E_Vj=(V_j.*d_V*mu_V )/w_V;        % J, structure in joules at the end larval period
  
  % Ana: below I calculate the joules in the body of the larva at the end of the
  % larval period, it should be around 100-150 joules(Howell and Fisher 1977)
  JoulesLarv=E_j+E_Rj+E_Sj+E_Vj; % J, total J in larva at the end larval period
    
  
  Wd_Sj = E_Sj/24370;% g, silk, there are aproximately 24370 joules in 1 g of dry weight of cocoon (Howell and Fisher 1977)
  Wd_j = Wd_Sj+(L_j^3 + (E_j+E_Rj) * w_E/ mu_E/ d_E) * d_V; % g, dry weight at pupation (with silk)
  
  del_S= Wd_Sj / Wd_j; % - , proportion weight of coccoon by weight at pupation with cocoon.

  
  % from pupation to emergence (from j to e)
  v_j = v * L_j/ L_b; % cm/d, energy conductance of imago
  L_mj = L_m * L_j/ L_b; % cm, max length of imago
  aVELR0_j =[0; L_j^3; E_j; 1e-4; E_Rj; 0]; %initial state vector at pupation [time,V_larv,E, L ,E_R ,E_eggs]
  [E_H aVELR0] = ode45(@dget_aVELR0, [0; E_He], aVELR0_j, [], k_J, k_E, v_j, E_m, L_mj, kap, g, M_V, y_EV_l, mu_E, E_H0, kap_R);
  t_e = aVELR0(end,1);     % d, time since pupation at emergence
  tT_e = t_e/ TC;        % d, duration of pupal stage at T
  V_e = aVELR0(end,2);     % cm^3, larval structure at emergence (should be close to zero)
  E_e = aVELR0(end,3);     % J, reserve at emergence
  L_e = aVELR0(end,4);     % cm, length at emergence
  E_Re = aVELR0(end,5);   % J, reproduction buffer at emergence
  E_0e = aVELR0(end,6);   % J, reserve in eggs at emergence
  E_e = E_e + y_EV_l * V_e * M_V * mu_E; % J, convert last bit of larval structure to reserve instantaneously
  
  Lw_e = L_e/ del_M;            % cm, physical length at puberty at f
  aT_p= aT_j+ tT_e;      % d, age at wasp emergence
  Wd_e = (L_e^3 + (E_e+E_Re +E_0e ) * w_E/ mu_E/ d_E)*d_V;   % g, dry weight at emergence, le sumamos el rep buffer construido en la larva
  
  
   % reproduction
  N_e = E_0e/ E_0; % total number of eggs (at emergence)
  
  % life span
  pars_tm = [g; l_T; h_a/ k_M^2; s_G];     % compose parameter vector at T_ref
  t_m = get_tm_s(pars_tm, f, l_b);         % -, scaled mean life span at T_ref
  aT_m = t_m/ k_M/ TC;                     % d, mean life span at T


 
    
%% uni-variate data
  
  % correct rates for temp
  kT_J = k_J * TC; kT_M = k_M * TC; kT_E = k_E * TC;  vT = v_j * TC; 
  kT_E_l=vT/L_j; % spec decay rate of larval structure in pupa
  % Weight-larva data
  
  ELRS_0 = [f * E_m * L_b^3; L_b; 0; 0]; % initial state vector [E; L; E_R; E_S]
  [t ELRS] = ode45(@dget_ELRS, tW_bj(:,1), ELRS_0, [], f, kT_M, kT_E, p_J*TC, p_Amb*TC, E_m, g, kap, kap_SA, kap_SG, kap_SD, p_M*TC, E_Rs);
  E_bj=ELRS(:,1);
  V_bj=ELRS(:,2).^3;
  E_R_bj=ELRS(:,3);
  E_S_bj = ELRS(:,4); % silk in Joules through time
  
  EW_S_bj = E_S_bj/24370;% g, silk, there are aproximately 24370 joules in 1 g of dry weight of cocoon (Howell and Fisher 1977)
  EW_bj = EW_S_bj+(V_bj + (E_bj+E_R_bj) * w_E/ mu_E/ d_E) * d_V; % g, dry weight at pupation (with silk)
  
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
    
  % Ana: col 1: E, col 2: L, col 3: R buffer
    % Weight-pupa data
    % pupal data
  % t-Wj-data
   VHELR0_j = [L_j^3; 0; E_j; 1e-4;E_Rj; 0]; % state at start pupation [V; E_H; E;L_imago;E_R; E_eggs]
  [t VHELR0] = ode45(@dget_VHELR0, tW_je(:,1), VHELR0_j, [], kT_J, kT_E_l, vT, E_m, L_mj, kap, g, M_V, y_EV_l, mu_E,E_H0,kap_R);
  V_l_je = VHELR0(:,1); E_je = VHELR0(:,3); L_je = VHELR0(:,4); V_je = L_je.^3; e_je = E_je ./ V_je/ E_m;  % -, scaled reserve density
  E_R_je=VHELR0(:,5);E_0_je=VHELR0(:,6);
  EW_je = Wd_Sj + d_V*((V_l_je + V_je) + (E_je+E_R_je+E_0_je)* w_E/ mu_E/ d_E); % g, dry weight of pupa  with cocoon
  
  kT_E_im=k_E_im*TC; 
  
  % Reproduction, progeny since wasp emergence
  pT_M = TC * p_M * L_e^3;  % J/d, maint rate (som)
  pT_J=TC * E_He * k_J;    % J/d, maint rate (mat)
  %pT_A =0;
  pT_A =pT_M+pT_J; % assim equals maint
  kT_E = k_E * TC; % 1/d, reserve mobilisation rate of adult
  
  ER0_e = [E_e,E_Re, E_0e]; %total reserve at wasp emergence, reprod buffer (excluding energy in eggs) and egg buffer
  [t ER0] = ode45(@dget_ER0, t_em(:,1), ER0_e, [], pT_A,pT_J,pT_M, kT_E_im, L_e,E_m,kap_R);
  

  EN_em = ER0(:,3)/ E_0;% #, cumulative number of eggs

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
 
  
  %Finally I would like to know the number of eggs produced when the wasp
  %face starving conditions after wasp emergence, note that it lives aprox only 2.53 days. The
  %total number of eggs produced two days after emergence should be around 75 accoring to Harvey et al. 2001
  timeStepWasp=0.001; %selected time step for wasp
  t_wasp_Star=0:timeStepWasp:2; % two days after emergence
  
  pT_A=0;
  %pT_A =pT_M+pT_J; % assim equals maint
  [t ER0] = ode45(@dget_ER0, t_wasp_Star, ER0_e, [], pT_A,pT_J,pT_M, kT_E_im, L_e,E_m,kap_R);
   E_e_starv=ER0(:,1);  % joules 
   E_R_em_starv=ER0(:,2);% joules 
   E_0_em_starv=ER0(:,3);% joules 
   
    EN_em_starv=E_0_em_starv(end)/E_0;%number of eggs after 2.53 days when they are supposed to die
    
    
  
   % pack output for zero-variate data
  Edata = [aT_b; aT_j; aT_p; aT_m; Wd_x; Wd_b; Wd_j; Wd_e; del_S; N_e; % real data
        v; kap; kap_R; p_M; p_T; k_J; kap_G];        % pseudo data
 
  
 end

 
      
  
function daEL = dget_aEL(H, aEL,k_J,p_M, v, E_m, L_m,k, g)
        % embryo development before it starts assimilating
        a = aEL(1); % d, time since the start of life
        E = aEL(2); % reserve  sincethe start of life
        L = aEL(3); % structural lenght the start of life
        
        V=L^3;
        
        
        e = E/ L^3/ E_m;      % -, scaled reserve density
        r = v * (e/ L - 1/ L_m)/ (e + g);  % 1/d, specific growth rate
        
        p_C = E * (v/ L - r); % J/d mobilisation flux
        p_J = k_J * H;        % J/d mat maintenance flux
        p_S = p_M * V;        % J/d somatic maintenance flux
        dE =-p_C;             % J/d; change in reserve
        dH = (1-k)*p_C-p_J;   % J/d, change in maturity
        dL = r * L/ 3;        % cm/d, change in length
        
        %% pack derivatives
       daEL = [1; dE; dL]/dH;
    end

    function daEL1 = dget_aEL1(H, aEL,k_J,p_M, v, E_m, L_m,k, g, p_Am,f)
        % embryo development after it starts assimilating
        a = aEL(1); %d, time since pre-birth
        E = aEL(2); % reserve at pre-birth
        L = aEL(3); % lenght at pre-birth
        
        V=L^3;
        
        
        e = E/ L^3/ E_m;      % -, scaled reserve density
        r = v * (e/ L - 1/ L_m)/ (e + g);  % 1/d, specific growth rate
        
        p_A = f * p_Am * V^(2/3); % J/d assimilation flux
        p_C = E * (v/ L - r); % J/d mobilisation flux
        p_J = k_J * H;        % J/d mat maintenance flux
        p_S = p_M * V;        % J/d somatic maintenance flux
        dE =p_A-p_C;             % J/d; change in reserve
        dH = (1-k)*p_C-p_J;   % J/d, change in maturity
        dL = r * L/ 3;        % cm/d, change in length
        
        %% pack derivatives
       daEL1 = [1; dE; dL]/dH;
    end
        
                                   
    function dELRS = dget_ELRS(t, ELRS, f, k_M, k_E, p_J, p_Amb, E_m, g, kap, kap_SA, kap_SG, kap_SD, p_M, E_Rs)
      % larval development
      % t: time since birth
      
      % unpack variables
      E  = ELRS(1);  % J, reserve
      L  = ELRS(2);  % cm, structural length
      E_R= ELRS(3);  % J, reproduction buffer
      E_S  = ELRS(4);  % J, silk
                
      V = L^3;                          % cm^3, structural volume
      e = E/ V/ E_m;                    % -, scaled reserve density
      r = (e * k_E - g * k_M)/ (e + g); % 1/d, specific growth rate
      p_C = E * (k_E - r);              % J/d, mobilisation rate
      p_S = p_M * V;                    % J/d, somatic maint rate during larval stage
      p_A = f * p_Amb * V;               % J/d, assimilation flux of larvae
      p_G = kap * p_C - p_S;            % J/d, growth rate (actual growth and overheads)
      p_D = p_J + p_S;                  % J/d, dissipation flux (somatic and maturity maint) not reprod overhead because it does not build eggs yet
      
      
      dE = p_A  - p_C;          % J/d, change in reserve
      dL = r * L/ 3;                    % cm/d, change in length
      dE_R = (1 - kap) * p_C - p_J;     % J/d, change in reprod buffer
      if E_R/V < E_Rs;
      dE_S = 0; % %joules/d, silk flux
      else 
      dE_S = kap_SA*p_A + kap_SG*p_G + kap_SD*p_D; %joules/d, silk flux
      end
          
      dELRS = [dE; dL; dE_R; dE_S]; 
     
    end
    



    function daELRS = dget_aELRS(ER, aELRS, f, k_M, k_E, p_J, p_Amb, E_m, g, kap, kap_SA, kap_SG, kap_SD, p_M, E_Rs)
      % larval  - development
      % ER condition index: [E_R] = E_R/V in J/cm^3
      
      % unpack variables
      %a  = aELRS(1); % d, time since birth
      E  = aELRS(2); % J, reserve
      L  = aELRS(3); % cm, structural length
      E_R= aELRS(4); % J, reproduction buffer
      E_S  = aELRS(5); % J, silk as product
      
      
       
      
      V = L^3;                          % cm^3, structural volume
      e = E/ V/ E_m;                    % -, scaled reserve density
      r = (e * k_E - g * k_M)/ (e + g); % 1/d, specific growth rate
      p_C = E * (k_E - r);              % J/d, mobilisation rate
      p_S = p_M * V;                    % J/d, somatic maint rate during larval stage
      p_A = f * p_Amb * V;               % J/d, assimilation flux of larvae
      p_G = kap * p_C - p_S;            % J/d, growth rate (actual growth and overheads)
      p_D = p_J + p_S;                  % J/d, dissipation flux (somatic and maturity maint) not reprod overhead because it does not build eggs yet
      
     
      dE = p_A  - p_C;          % J/d, change in reserve
      dL = r * L/ 3;                    % cm/d, change in length
      dE_R = (1 - kap) * p_C - p_J;     % J/d, change in reprod buffer
      dER = dE_R/ V - r * E_R/ V;       % J/d.cm^3, change in [E_R]
      if E_R/V < E_Rs;
      dE_S  = 0; % %joules/d, silk flux
      else 
      dE_S  = kap_SA*p_A + kap_SG*p_G + kap_SD*p_D; % J/d, silk flux
      end
      daELRS = [1; dE; dL; dE_R; dE_S]/ dER;  
      
    end

    function daVELR0 = dget_aVELR0(E_H, aVELR0, k_J, k_E, v, E_m, L_m, kap, g, M_V, y_EV, mu_E, E_H0, kap_R)
      % pupal development
      % E_H: maturity of imago (J)
      %   y_EV: yield of reserve of imago on structure of larva
       
      % unpack variables
      %a  = aVELR0(1); % d, time since pupation
      V  = aVELR0(2); % cm^3, structural volume of larva
      E  = aVELR0(3); % J, reserve of pupa
      L  = aVELR0(4); % cm, structural length of imago
      E_R = aVELR0(5); % J,reproduction buffer, excluding eggs
      E_0 = aVELR0(6);% J, reserve in eggs
            
      dV = - V * k_E;                    % cm^3/d, change larval structure
      e = E/ L^3/ E_m;                   % -, scaled reserve density
      r = v * (e/ L - 1/ L_m)/ (e + g);  % 1/d, specific growth rate
      p_C = E * (v/ L - r);              % J/d, mobilisation rate
      dE = - dV * M_V * y_EV * mu_E - p_C; % J/d, change in reserve
      dL = r * L/ 3;                     % cm/d, change in length
      dE_H = (1 - kap) * p_C - k_J * E_H;    %J/d; change in maturity
     
      if E_H < E_H0;
      dE_R = 0; % J, cumulative investment in reprod buffer
      dE_0 = 0; % J, cumulative investment in eggs
      else        
        p_C2 = E_R * k_E;          % J/d, total drain from E_R to eggs
        dE_R = - p_C2;             % J, cumulative investment in reprod buffer
        dE_0 = kap_R * p_C2; % J, cumulative investment in eggs      
        
      end
      daVELR0 = [1; dV; dE; dL; dE_R; dE_0]/ dE_H; % pack derivatives
    end
                                      
 function dVHELR0 = dget_VHELR0(a, VHELR0, k_J, k_E, v, E_m, L_m, kap, g, M_V, y_EV, mu_E, E_H0, kap_R)
      % pupal development
      % a: d, time since pupation (before spinning)
      %   y_EV: yield of reserve of imago on structure of larva
      %   y_SV: yield of silk on structure of larva
       
      % unpack variables
      V  = VHELR0(1); % cm^3, structural volume of larva
      E_H  = VHELR0(2); % J, maturity 
      E  = VHELR0(3); % J, reserve of larva
      L  = VHELR0(4); % cm, structural length of imago
      E_R = VHELR0(5); % J,reprod buffer, excluding eggs
      E_0 = VHELR0(6);% J, reserve in eggs
            
      dV = - V * k_E;                     % cm^3/d, change larval structure
      e = E/ L^3/ E_m;                    % -, scaled reserve density
      r = v * (e/ L - 1/ L_m)/ (e + g);   % 1/d, specific growth rate
      p_C = E * (v/ L - r);               % J/d, mobilisation rate
      dE = - dV * M_V * y_EV * mu_E - p_C;  % J/d, change in reserve
      dL = r * L/ 3;                      % cm/d, change in length
      dE_H = (1 - kap) * p_C - k_J * E_H;     % J/d, change in maturity
      
       if E_H < E_H0;
      dE_R = 0; % J, cumulative investment in reprod buffer
      dE_0 = 0; % J, cumulative investment in eggs
      else        
        p_C2 = E_R * k_E;          % J/d, total drain from E_R to eggs
        dE_R = - p_C2;             % J, cumulative investment in reprod buffer
        dE_0 = kap_R * p_C2; % J, cumulative investment in eggs      
        
      end
      dVHELR0 = [dV; dE_H; dE; dL; dE_R; dE_0];
  end
  
  function dER0 = dget_ER0(a, ER0, p_A, p_J, p_M, k_E, L_e,E_m,kap_R)
        % change in reserves during the wasp life
        % a: d, time since wasp emergence
        
        %unpack variables
        E = ER0(1); % J, reserves of the wasp
        E_R = ER0(2);% J,reprod buffer, excluding eggs
        E_0 = ER0(3);% J, reserve in eggs
        
        V = L_e^3; 
        e = E/ V/ E_m; 
        p_C = E*k_E;
        p_S = p_M + p_J;
        p_R = p_C-p_S;
        p_C2 = E_R*k_E; %total drain from ER to eggs
        p_Ro = (1-kap_R)*p_C2; %eggs overhead
        p_R2 = p_C2-p_Ro;% flux of energy that goes to eggs
        
        dE = p_A - p_C;   % J/d, change in reserve
        dE_R = max(0,p_R)-p_C2; % J, cumulative investment in reprod buffer
        dE_0 = max (0,p_R2); % J, cumulative investment in eggs
        
        dER0 = [dE; dE_R; dE_0];
   end