% fct calc_epsilon. calculates strain with model 'WibCreep V5.1'    
function solu_calc_epsilon = calc_epsilon( time, sigma, epsilon0, pE0, pE1, ...
    peta1, peta2, peta3, peta4, pE01d, peta11d, peta21d, peta31d, ...
    peta41d, kkorr)

  % estimate material parameters
  E0 = e_aging(pE0, pE01d, time(1));
  E2 = pE1;
  E3 = pE1;
  E4 = pE1;
  eta1 = eta_aging(peta1, peta11d, time(1));
  eta2 = eta_aging(peta2, peta21d, time(1));
  eta3 = eta_aging(peta3, peta31d, time(1));
  eta4 = eta_aging(peta4, peta41d, time(1));

  % estimate coefficients
  a0 = 1/(E0*eta1*eta2*eta3*eta4);
  a1 = 1/(E0*pE1*eta2*eta3*eta4) + 1/(E0*E2*eta1*eta3*eta4) + ...
      1/(E0*E3*eta1*eta2*eta4) + 1/(E0*E4*eta1*eta2*eta3);
  a2 = 1/(E0*pE1*E2*eta3*eta4) + 1/(E0*pE1*E3*eta2*eta4) + ...
      1/(E0*pE1*E4*eta2*eta3) + 1/(E0*E2*E3*eta1*eta4) + ...
      1/(E0*E2*E4*eta1*eta3) + 1/(E0*E3*E4*eta1*eta2);
  a3 = 1/(E0*pE1*E2*E3*eta4) + 1/(E0*pE1*E2*E4*eta3) + ...
      1/(E0*pE1*E3*E4*eta2) + 1/(E0*E2*E3*E4*eta1);
  a4 = 1/(E0*pE1*E2*E3*E4);
  b0 = 1/(eta1*eta2*eta3*eta4);
  b1 = 1/(E0*eta1*eta2*eta3) + 1/(E0*eta1*eta2*eta4) + ...
      1/(E0*eta1*eta3*eta4) + 1/(E0*eta2*eta3*eta4) + ...
      1/(pE1*eta2*eta3*eta4) + 1/(E2*eta1*eta3*eta4) + ...
      1/(E3*eta1*eta2*eta4) + 1/(E4*eta1*eta2*eta3);
  b2 = 1/(E0*pE1*eta2*eta3) + 1/(E0*E2*eta1*eta3) + ...
      1/(E0*E3*eta1*eta2) + 1/(E0*pE1*eta2*eta4) + ...
      1/(E0*E2*eta1*eta4) + 1/(E0*E4*eta1*eta2) + ...
      1/(E0*pE1*eta3*eta4) + 1/(E0*E3*eta1*eta4) + 1/(E0*E4*eta1*eta3) + ...
      1/(E0*E2*eta3*eta4) + 1/(E0*E3*eta2*eta4) + 1/(E0*E4*eta2*eta3) + ...
      1/(pE1*E2*eta3*eta4) + 1/(pE1*E3*eta2*eta4) + ...
      1/(pE1*E4*eta2*eta3) + 1/(E2*E3*eta1*eta4) + 1/(E2*E4*eta1*eta3) + ...
      1/(E3*E4*eta1*eta2);
  b3 = 1/(E0*pE1*E2*eta3) + 1/(E0*pE1*E3*eta2) + 1/(E0*E2*E3*eta1) + ...
      1/(E0*pE1*E2*eta4) + 1/(E0*pE1*E4*eta2) + 1/(E0*E2*E4*eta1) + ...
      1/(E0*pE1*E3*eta4) + 1/(E0*pE1*E4*eta3) + 1/(E0*E3*E4*eta1) + ...
      1/(E0*E2*E3*eta4) + 1/(E0*E2*E4*eta3) + 1/(E0*E3*E4*eta2) + ...
      1/(pE1*E2*E3*eta4) + 1/(pE1*E2*E4*eta3) + 1/(pE1*E3*E4*eta2) + ...
      1/(E2*E3*E4*eta1);
  b4 = 1/(E0*pE1*E2*E3) + 1/(E0*pE1*E2*E4) + 1/(E0*pE1*E3*E4) + ...
      1/(E0*E2*E3*E4) + 1/(pE1*E2*E3*E4);

  % estimate initial conditions
  eps0    = a4/b4*sigma;
  deps0   = (a3*sigma - b3*eps0)/b4;
  ddeps0  = (a2*sigma - b2*eps0 - b3*deps0)/b4;
  dddeps0 = (a1*sigma - b1*eps0 - b2*deps0 - b3*ddeps0)/b4;

  eps0    = kkorr*eps0 + epsilon0(1);
  deps0   = kkorr*deps0 + epsilon0(2);
  ddeps0  = kkorr*ddeps0 + epsilon0(3);
  dddeps0 = kkorr*dddeps0 + epsilon0(4);

  % solve ode
  [t_calc, eps_calc] = ode23s(@(t,eps) odefun_materialmodel(t, eps, pE0, ...
    pE1, peta1, peta2, peta3, peta4, pE01d, peta11d, peta21d, peta31d, ...
    peta41d, sigma), time, [dddeps0 ddeps0 deps0 eps0]);

  % return solution
  solu_calc_epsilon(:,1) = eps_calc(:,4);
end
