functions {

  real coupling(real e_inner, real e_outer, real A_inner, real A_outer) {
    
    real E;
    E = e_inner*e_outer / (e_outer + (A_inner/A_outer)*(e_inner - e_inner*e_outer));
    return E;
  }
  
  real Cp_Si(real T) {
    real cp = -82.94826 + 2.712235*T + 0.01404751*T^2 - 7.977691e-5*T^3 + 1.079905e-7*T^4;
    return cp;
  }

  real Cp_Al(real T) {
    real cp = 0.0000459067716*T^3 + -0.0385698040*T^2 + 11.5824787*T - 341.887005;
    return cp;
  }
  
  real k_Cu(real T) {
    num =  2.8075 - 1.2777*T^0.5 + 0.36444*T - 0.051727*T^1.5 + 0.0030964*T^2 
    den =  1 - 0.54074*T^0.5 + 0.15362*T - 0.02105*T^1.5 + 0.0012226*T^2
    return 10^(num/den)
  }

  real[] cooldown(real t,
                  real[] y,
                  real[] e_blackcoat,
                  real[] x_r,
                  int[] x_i) {

    real dydt[3];
    real P_is_tm;
    real P_os_is;
    real P_ch_is;
    real T_ch = 45.769;
    
    P_ch_is = k_Cu(T_ch)*0.000370967/0.38443215119969776 * (y[2] - T_ch);
    P_is_tm = coupling(e_blackcoat[1], e_blackcoat[1], 0.0486, 0.7509) * 5.670374419e-8 * (y[1]^4 - y[2]^4);

    dydt[1] = -1*P_is_tm/(Cp_Si(y[1])*1.92);
    dydt[2] = (-P_ch_is + P_is_tm)/(Cp_Al(y[2])*21.74);
    return dydt;
  }
}

data {
  int<lower=1> T;
  real y[T,4];
  real y0[2];
  real t0;
  real ts[T];
}

transformed data {
  real x_r[0];
  int x_i[0];
}

parameters {
  real<lower=0, upper=1> e_blackcoat[1];
  vector<lower=0>[4] sigma;
}


model {
  real y_hat[T,4];
  sigma ~ cauchy(0, 2.5);
  
  e_blackcoat ~ normal(0.7,0.05);
  
  y_hat = integrate_ode_rk45(cooldown, y0, t0, ts, e_blackcoat, x_r, x_i);
  for (t in 1:T) {
    y[t] ~ normal(y_hat[t], sigma);
  }
}