function [prdData, info] = predict_Bactrocera_oleae(par, data, auxData)
% file generated by prt_predict

% unpack par, data, auxData
cPar = parscomp_st(par); vars_pull(par);
v2struct(par); v2struct(cPar); v2struct(data); v2struct(auxData);

if T_H < T_ref || T_AH < 0 ||  T_AL < 0 || T_L > T_ref || f_3 >1 || f_4 >1
    prdData = []; info = 0; return
end

% compute temperature correction factors
pars_T = [T_A, T_L, T_H, T_AL, T_AH];

TC_ab = tempcorr(temp.ab, T_ref, pars_T);
TC_am = tempcorr(temp.am, T_ref, pars_T);
TC_tj = tempcorr(temp.tj, T_ref, pars_T);
TC_te = tempcorr(temp.te, T_ref, pars_T);
TC_Ri = tempcorr(temp.Ri, T_ref, pars_T);
TC_25 = tempcorr(C2K(25), T_ref, pars_T);
TC_26 = tempcorr(temp.tS_cold, T_ref, pars_T);
TC_14 = tempcorr(temp.tL_14, T_ref, pars_T);
TC_21 = tempcorr(temp.tL_21, T_ref, pars_T);

% life cycle
pars_tj = [g k 0  v_Hb v_Hp v_Hp + 1e-08 v_He kap kap_V];
[t_j, t_e, t_b, l_j, l_e, l_b, rho_j, u_Ee, info] = get_tj_dacos(pars_tj, f);
if info == 0
  prdData = []; return;
end

 % initial
  pars_UE0 = [V_Hb; g; k_J; k_M; v]; % compose parameter vector
  E_0 = p_Am * initial_scaled_reserve(f, pars_UE0); % d.cm^2, initial scaled reserve
  Ww_0 = 1e3 * w_E * E_0/ mu_E;     % mg, initial wet weight

% birth
  aT_b = t_b / k_M/ TC_ab;          % d, age at birth f and T
  L_b = L_m * l_b;                  % cm, structural length at birth at f
  Lw_b = L_b/del_M; 
  Ww_b = L_b^3 * (1 + f * w);       % g, wet weight at birth 
  
 % pupation 
  tT_j = t_j / k_M/ TC_tj;          % d, time since birth at puberty at f and T
  L_j = L_m * l_j; Lw_j = L_j/del_M;                 % cm, structural length at puberty at f
  Ww_j = L_j^3 * (1 + f * w);       % g, wet weight at puberty 

 % emergence 
  tT_e = (t_e-t_j) / k_M/ TC_te;          % d, time since birth at puberty at f and T
  L_e = L_m * l_e;  Lw_e = L_e/del_M;                % cm, structural length at emergence at f

 % body lenth at larval instars
  L_1 = L_b * s_1^0.5;  Lw_1 = L_1 / del_M; % cm, length at end of instar 1
  L_2 = L_1 * s_2^0.5; Lw_2 = L_2 / del_M; % cm,  length at end of instar 2
  L_3 = L_j; Lw_3 = L_3 / del_M; % cm,  length at end of instar 3

 % life span

  pars_tm = [g; k; v_Hb; v_Hp; 0; v_He; kap; kap_V; h_a/k_M^2; s_G]; % compose parameter vector at T_ref
  [t_m, a, b]= get_tm_mod_dacos('hax',pars_tm, f);      % -, scaled mean life span at T_ref
  aT_m = (t_m - t_j) / k_M/TC_am;         % d, mean life span as an imago at T

  
%  reproduction (kappa-rule)
 s_M = L_j/L_b;
 R_i = (kap_R * ((1-kap) * f * s_M * v * E_m * L_e^2 - k_J * E_Hp)/ E_0)* TC_Ri; % #/d, ultimate reproduction rate

  % pack to output
  prdData.ab = aT_b;
  prdData.tj = tT_j;
  prdData.te = tT_e;
  prdData.am = aT_m;
  prdData.Ww0 = Ww_0;
  prdData.Wwj = Ww_j;
  prdData.Ri = R_i;

  prdData.L1 = Lw_1;
  prdData.L2 = Lw_2;
  prdData.L3 = Lw_3;

  prdData.Lj = Lw_j;
  prdData.Le = Lw_e;

% data for different temperatures
  TC = tempcorr(C2K(T_ab(:,1)), T_ref, pars_T);
  ab_T = t_b / k_M./ TC;
  r_embryo = 1./ ab_T;
  tT = (t_j-t_b) / k_M./TC;
  r_larva = 1./tT;
  TC = tempcorr(C2K(T_te(:,1)), T_ref, pars_T);
  te_T = (t_e-t_j)/ k_M./TC;
  r_pupa = 1./ te_T;

  % reproduction
  pars_T = [p_M; k_J; E_G; E_Hp; v; p_Am; kap; h_a; s_G; E_Hb; 0; L_b; L_j; E_He; kap_V];
  ELHE_0 = [f * E_m * L_e^3; L_e; E_He; 0; 0; 0; 1]; % initial state vector
  [t ELHE] = ode45(@dget_ELHE_surv_imago, [t_e/k_M/TC_25; tR(:,1)+t_e/k_M/TC_25], ELHE_0, [], pars_T, f, TC_25);
  EN = kap_R * ELHE(:,4)/ E_0;     % #, cumulative reprod output
  EN(1) = []; t(1) = [];

  % time - survival
  h_T_egg = max(0,h_Tb * ( C2K(TS_b(:,1))/T_H - 1));
  h_T_pupa = max(0,h_Tp * (C2K(TS_p(:,1))/T_H - 1));
  h_T_imago = max(0,h_Ti * (C2K(TS_i(:,1))/T_H - 1));

  % 70, 875, 85 and 94 denote the survival at the lowest temp
  ES_egg = 70*exp (-h_T_egg  .* t_b);
  ES_pupa = 94*exp (-h_T_pupa .* (t_e-t_j));
  ES_imago = 100*exp (-h_T_imago .* (t_m-t_e));

  % adult survival 
  ELHE_0 = [f * E_m * L_b^3; L_b; E_Hb; 0; 0; 0; 0.89]; % initial state vector
  pars_T = [p_M; k_J; E_G; E_Hp; v; p_Am; kap; h_a; s_G; E_Hb; h_i; L_b; L_j; s_1; s_2];  
  [t ELHE] = ode45(@dget_ELHE_surv_imago, [t_e/k_M/TC_26; tS_cold(:,1)+t_e/k_M/TC_26], ELHE_0, [], pars_T, f, TC_26);
  ES_cold = ELHE(:,7) * 100; ES_cold(1)=[];    % #, survival 

  % length at various temperatures
  opts = odeset('RelTol', 1e-4, 'AbsTol', 1e-8);
  ELH_0 = [f * E_m * Lb_tL_T^3;Lb_tL_T; E_Hb]; % initial state vector

  [t ELH_14] = ode45(@dget_ELH,[0;tL_14(:,1)], ELH_0, opts, pars_T, f, TC_14);
  ELH_14(1,:) = []; t(1) = [];
  L_14 = 10 * ELH_14(:,2)/del_M; % mm, length

  [t ELH_21] = ode45(@dget_ELH,[0;tL_21(:,1)], ELH_0, opts, pars_T, f, TC_21);
  ELH_21(1,:) = []; t(1) = [];
  L_21 = 10 * ELH_21(:,2)/del_M; % mm, length

  [t ELH_26] = ode45(@dget_ELH,[0;tL_26(:,1)], ELH_0, opts, pars_T, f, TC_26);
  ELH_26(1,:) = []; t(1) = [];
  L_26 = 10 * ELH_26(:,2)/del_M; % mm, length

  % time - wet weight
  ELH_0 = [f * E_m * L_b^3; L_b; E_Hb]; % initial state vector
  opts = odeset('RelTol', 1e-4, 'AbsTol', 1e-8);
  [t1 ELH1] = ode45(@dget_ELH, [0;  tWw_f(:,1)], ELH_0, opts, pars_T, f, TC_25);
  ELH1(1,:) = []; t1(1) = [];
  EWw_1 = ELH1(:,2).^3 + ELH1(:,1)*(w_E/mu_E/d_E);  % g, wet weight 

  [t2 ELH2] = ode45(@dget_ELH, [0;  tWw_f(:,1)], ELH_0, opts, pars_T, f, TC_25);
  ELH2(1,:) = []; t2(1) = [];
  EWw_2 = ELH2(:,2).^3 + ELH2(:,1)*(w_E/mu_E/d_E);  % g, wet weight

  [t3 ELH3] = ode45(@dget_ELH, [0;  tWw_f(:,1)], ELH_0, opts, pars_T, f_3, TC_25);
  ELH3(1,:) = []; t3(1) = [];
  EWw_3 = ELH3(:,2).^3 + ELH3(:,1)*(w_E/mu_E/d_E);  % g, wet weight 

  [t4 ELH4] = ode45(@dget_ELH, [0;  tWw_f(:,1)], ELH_0, opts, pars_T, f_4, TC_25);
  ELH4(1,:) = []; t4(1) = [];
  EWw_4 = ELH4(:,2).^3 + ELH4(:,1)*(w_E/mu_E/d_E);  % g, wet weight


%  pack to output 
  prdData.T_ab = ab_T;
  prdData.T_tj = tT;
  prdData.T_te = te_T;

  prdData.T_re = r_embryo;
  prdData.T_rl = r_larva;
  prdData.T_rp = r_pupa;

  prdData.TS_b = ES_egg;
  prdData.TS_p = ES_pupa;
  prdData.TS_i = ES_imago;

  prdData.tR = EN;
  prdData.tS_cold = ES_cold;

  prdData.tL_14 = L_14;
  prdData.tL_21 = L_21;
  prdData.tL_26 = L_26;

  prdData.tWw_f = [EWw_1 EWw_2 EWw_3 EWw_4];

end

function dELH = dget_ELH(t, ELH, p, f, TC)
  % larval development
  % t: time since birth
      
  % unpack variables    
  E  = ELH(1); % J, reserve
  L  = ELH(2); % cm, structural length
  EH = ELH(3); % J, maturity

  % unpack parameters 
  p_M = p(1); k_J = p(2); E_G = p(3); v = p(5); 
  p_Am = p(6); kap = p(7); h_a = p(8); L_b = p(12); L_j = p(13);
  V = L^3;  % cm^3, structural volume

  pT_M = TC * p_M;  pT_Am = TC * p_Am;   kT_J = TC * k_J;   vT = TC * v; 

  % juvenile;  
  s_M = L / L_b;
      
  r = (E * vT * s_M/ L - pT_M*V/kap)/(E + E_G*V/kap);       % 1/d, specific growth rate 
  p_C = E * (E_G * vT* s_M / L  + pT_M)/ (kap*E/V + E_G);  % J/d, mobilisation rate
  % generate derivatives
  dE = f * pT_Am * s_M* L^2 - p_C;  % J/d, change in reserve
  dL = r * L/ 3 * (L<L_j);                             % cm/d, change in length
  dEH = (1 - kap) * p_C - kT_J * EH; % J/d, change in maturity
      
  dELH = [dE; dL; dEH];
end


function dELHE_imago = dget_ELHE_surv_imago(t, ELHE, p, f, TC)
  % t: time since birth
      
  % unpack variables    
  E  = ELHE(1); % J, reserve
  L  = ELHE(2); % cm, structural length
  EH = ELHE(3); % J, maturity
  ER = ELHE(4); % J, reproduction buffer
  q = ELHE(5);  % 1/d^2, acceleration
  h = ELHE(6);  % 1/d,  hazard rate
  S = ELHE(7);  % -, surviving fraction

  % unpack parameters 
  p_M = p(1); k_J = p(2); v = p(5); 
  p_Am = p(6); kap = p(7); h_a = p(8); s_G = p(9); h_t = p(11); L_b = p(12); L_p = p(13);
  L_m = kap*p_Am/p_M;
  V = L^3;  % cm^3, structural volume
  E_m = p_Am/v;
  e = E/V/E_m; %-. scaled reserve density

  pT_Am = TC * p_Am;   kT_J = TC * k_J;   vT = TC * v; hT = TC * h; hT_a = TC * TC * h_a;

  % adult
  s_M = L_p / L_b;
  pC  = E * vT * s_M/ L; % J/d, mobilized energy flux
  % generate derivatives
  dE    = f * pT_Am * L^2 * s_M  - pC;  % J/d, change in energy in reserve
  dL    = 0;  % cm/d, change in structural length
  dEH  = 0; % J/d, change in cumulated energy invested in maturation 
  dER  =  (1-kap)*pC - kT_J * EH; % J/d, change in reproduction buffer
  % aging
  dq = (q * L^3 * s_G/ L_m^3  + hT_a)* e * vT * s_M /L;
  dh = q;
  % survival
  dS = -(hT + h_t) * S;

  dELHE_imago = [dE; dL; dEH; dER; dq; dh; dS];
end
