function lcval = lc_kilonova(x, t0, kappa_pp, N_data, tvec_days, dotrap)

powermethod = 'metzger';
slopesteep = 1;
pscale = 1.2;

M_ej = x(1);  % Ejecta mass (Msun)
E50 = x(2);  % Magnetic field (10^14 G)
E_51 = E50*0.1
kappa = x(3);  % Initial spin period (ms)
slope = x(4);  % Shift in days of light curve (t_sinceexplosion = tvec_days + shift_days)

% Constants
c      = 2.998e10  ; % Speed of light (cm/s) CHECKED
e_ni0   = 3.90e10   ; % Energy release in Ni56 decay (erg/s/g)  CHECKED
e_co0   = 6.78e9    ; % Energy release in Co56 decay (erg/s/g)  CHECKED
m_sol  = 2e33   ; % Solar mass (g)  CHECKED
tau_ni = 6.077*24*60*60/log(2)   ; % Decay time (half life / ln2) for Ni56 (s) 6.077  CHECKED
tau_co = 77.27*24*60*60/log(2)   ; % Decay time (half life / ln2) for Co56 (s) 77.27d  CHECKED
kappa_gamma = 0.03;  % CHECKED

Ntimes = 100

% Calculate tau_m, y and s (for photospheric phase)

V_ej = sqrt(10./3*E_51*1e51/(M_ej*2e33)); % Ejecta velocity, constant density  % CHECKED
 
tau_m = 1.05*((kappa/(13.7*c))^0.5) * (((((M_ej*m_sol)^3))/(E_51*1e51))^0.25);    % Diffusion time (Arnett 1982) Eq 18, 19, 22, 23 CHECKED

y = tau_m/(2*tau_ni);   % Arnett 1982 Eq 33 CHECKED

yp = tau_m/(2*tau_co);    % Arnet 1982 Eq 33, modified to 56Co decay  CHECKED

chi2 = 0;  % Initialize chi2 sum

Nintegrate = 2000;  % Number of time steps to run integrals over

for i = 1:Ntimes

  t = tvec_days(i)*24*3600         % Time in seconds
                                                                
  x = t/tau_m;     % Arnett 1982 Eq 32 CHECKED

  z = linspace(0.0001,x,Nintegrate);        % Define limits of intergration for A(z)  CHECKED

  % Compute gamma ray deposition
  R = V_ej*z*tau_m;     % Radius vector...z*tau_m = time;  CHECKED

  rho = M_ej*2e33./(4*pi/3*R.^3); % CHECKED

  tau_56co_gamma = kappa_gamma*rho.*R;  % CHECKED

  G = tau_56co_gamma./(tau_56co_gamma + 1.6);    % Arnett 1982 Eq 51 CHECKED

  %D_gamma = 0.1 + 0.9*(0.25 + 0.75*G.*(1 + 2*G.*(1-G).*(1-0.75*G)));    % Arnett 1982 eq 50  CHECKED
  %D_gamma=ones(Nintegrate);

  % Kilonova part

  power = zeros(Nintegrate,1);

  slopeuse = 0.;

  if (powermethod == 'anders1')

      ind = find((z*tau_m > 0.001*24*3600) & (z*tau_m < 10*24*3600)); % One interval is 0.0001-10d

      slopeuse = slope;

      power(ind) = pscale*1.6e10*(M_ej*m_sol)*(z(ind)*tau_m/(t0*24*3600)).^(slopeuse);%.*D_gamma(ind);

      ind = find(z*tau_m >= 10*24*3600); % Late time break

      slopeuse = slope*slopesteep;
      %slopeuse = slope;

      power(ind) = 10^(slope-slopeuse)*pscale*1.6e10*(M_ej*m_sol)*(z(ind)*tau_m/(t0*24*3600)).^(slopeuse);%.*D_gamma(ind);

elseif (powermethod == 'metzger')
 
      slopeuse = slope;

      %eth = ones(1,Nintegrate);
      %ts = 1.3; % 1.3 % s
      %sigma = 0.11; % s

      eth = 0.36*(exp(-0.56*z*tau_m/(24*3600)) + (log(1 + 2*0.17*(z*tau_m/(24*3600)).^0.74))/(2*0.17*(z*tau_m/(24*3600)).^0.74));

      %power = 4e18*eth.*(0.5 - 1/3.14*atan((z*tau_m-ts)/sigma)).^1.3
	
      ind = find(z*tau_m/(24*3600) > 0.0001 & (z*tau_m/(24*3600) < 10));
      power(ind) = eth(ind).*pscale*1.6e10.*(M_ej*m_sol).*(z(ind)*tau_m/(t0*24*3600)).^(slopeuse);
      %power(ind) =  4e18*eth(ind).*(0.5 - 1./3.1415*atan((z(ind)*tau_m-ts)/sigma)).^1.3;

      ind = find(z*tau_m/(24*3600) >= 10);

      slopeuse = slopesteep*slope

      power(ind) = 10^(slope-slopeuse)*eth(ind).*pscale*1.6e10*(M_ej*m_sol).*(z(ind)*tau_m/(t0*24*3600)).^(slopeuse);
      %power(ind) =  4e18*eth(ind).*(0.5 - 1./3.1415*atan((z(ind)*tau_m-ts)/sigma)).^1.3;

  end

  if (tvec_days(i) <= 2.5*tau_m/(24*3600)) % TEMP

      integrand_rprocess = power'.*exp(z.^2-x^2).*2.*z;

      %integrand_56ni= M_ej*m_sol*e_ni0*exp(-x^2)*2*z.*exp((-2*z*y)+(z.^2)).*D_gamma +  M_ej*m_sol*e_co0*exp(-x^2)*2*z.*(exp(-2*z*yp) - exp(-2*z*y))/(1 - tau_ni/tau_co).*exp(z.^2).*(0.966*D_gamma + 0.034*D_gamma);

      Lambda_kilonova = sum(integrand_rprocess.*(x/Nintegrate));


   else

     Lambda_kilonova = power(Nintegrate);

     %Lambda_kilonova = M_ej*m_sol*e_ni0*exp(-tvec_days(i)/8.8)*D_gamma(Nintegrate) +  M_ej*m_sol*e_co0*(exp(-tvec_days(i)*24*3600/tau_co) - exp(-tvec_days(i)*24*3600/tau_ni))/(1 - tau_ni/tau_co)*(0.966*D_gamma(Nintegrate) + 0.034*D_gamma(Nintegrate));
    
   end


    L_ph_1(i) = Lambda_kilonova;   % Calculate luminosity
     
end

lcval = L_ph_1
