function chi2 = chi2function(x, M_ni, E_51, kappa, kappa_pp, N_data, tvec_days, ydata_log, sigma_log, uselim, Nlimits, tvec_days_limits, ydata_log_limits, redshift, dotrap)

% COMPUTES THE CHI2 ERROR BETWEEN THE DATA (N_data, tvec_days, ydata_log, sigma_log), and
% a specified magnetar model (x, M_Ni, E_51, kappa)

M_ej = x(1);  % Ejecta mass (Msun)
B_14 = x(2);  % Magnetic field (10^14 G)
P_ms = x(3);  % Initial spin period (ms)
shift_days = 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

% 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 = 1500;  % Number of time steps to run integrals over

for i = 1:N_data         % TEMP Loop over time step i..for each time the solution is an integral over all previous time

    t = (tvec_days(i)/(1+redshift)+shift_days)*24*3600;          %
                                                                 % Get
                                                                 % time
                                                                 % since
                                                                 % explosion
                                                                 % in
                                                                 % seconds
                                                                 % (t)
                                                                 % from
                                                                 % i
                                                                 % in
                                                                 % restframe
                                                                 % CHECKED. Note
                                                                 % that
                                                                 % shift_days
                                                                 % is
                                                                 % in
                                                                 % rest frame
 
    % Calculate parameters for photospheric phase

    x = t/tau_m;     % Arnett 1982 Eq 32 CHECKED

    z = linspace(0.000001,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 = G.*(1 + 2*G.*(1-G).*(1-0.75*G));    % Arnett 1982 eq 50  CHECKED

    % 56Ni part: -------------------------------
  
    integrand1= exp(-x^2)*2*z.*exp((-2*z*y)+(z.^2)).*D_gamma;    % Define function A(z) Arnett 1982 Eq 31  CHECKED

    Lambda1 = sum(integrand1.*(x/Nintegrate));                 % Numerically integrate A(z)  Results checked against table 1 in Arnett 1982 CHECKED
   
    % 56Co part: --------------------------------
   
    integrand2 = 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);   % Define function B(z) (extension of Eq 31 in Arnett 1982)  CHECKED
 
    Lambda2 = sum(integrand2.*(x/Nintegrate));                 % Numerically integrate B(z)

    % Magnetar part: ------------------------------

    Ep = 1.97e52*P_ms^-2;  % Kasen Bildsten Eq 1  CHECKED

    tp = 0.013*365*24*3600*B_14^(-2)*P_ms^2;  % Kasen Bildsten Eq 2
                                              % (and converted to
                                              % seconds)  CHECKED
    if (dotrap == 1)

        % 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_pp = kappa_pp*rho.*R;  % CHECKED

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

        D_gamma = G.*(1 + 2*G.*(1-G).*(1-0.75*G));    % Arnett 1982
                                                      % eq 50
                                                      % CHECKED
    else

        D_gamma = 1.;   % No trapping

    end

    Cfunct = Ep/tp*(2-1)./(1 + z*tau_m/tp).^2.*D_gamma*exp(-x^2).*exp(z.^2).*2.*z;   % Modified Arnett 1982 Eq 31 CHECKED

    C = sum(Cfunct.*(x/Nintegrate));   % Integrate

    L_ph_1(i) = M_ni*m_sol*e_ni0*Lambda1 + M_ni*m_sol*e_co0*Lambda2 + C;   % Calculate luminosity 

    chi2 = chi2 + (log10(L_ph_1(i)) - ydata_log(i))^2/sigma_log(i)^2;   % Error in log  CHECKED
    %chi2 = chi2 + (L_ph_1(i) - 10^ydata_log(i))^2/(10^ydata_log(i)*(10^sigma_log(i) - 1))^2;  % Error in linear

end

if (uselim > 0)

    for i = 1:Nlimits

        t = (tvec_days_limits(i)/(1+redshift)+shift_days)*24*3600;          % Get time since explosion in seconds (t) from i in restframe  CHECKED
 
        % Calculate parameters for photospheric phase

        x = t/tau_m;     % Arnett 1982 Eq 32 CHECKED

        z = linspace(0.000001,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 = G.*(1 + 2*G.*(1-G).*(1-0.75*G));    % Arnett 1982 eq 50  CHECKED

        % 56Ni part: -------------------------------
  
        integrand1= exp(-x^2)*2*z.*exp((-2*z*y)+(z.^2)).*D_gamma;    % Define function A(z) Arnett 1982 Eq 31  CHECKED

        Lambda1 = sum(integrand1.*(x/Nintegrate));                 % Numerically integrate A(z)  Results checked against table 1 in Arnett 1982 CHECKED
        
        % 56Co part: --------------------------------
   
        integrand2 = 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);   % Define function B(z) (extension of Eq 31 in Arnett 1982)  CHECKED
        
        Lambda2 = sum(integrand2.*(x/Nintegrate));                 % Numerically integrate B(z)

        % Magnetar part: ------------------------------

        Ep = 1.97e52*P_ms^-2;  % Kasen Bildsten Eq 1  CHECKED

        tp = 0.013*365*24*3600*B_14^(-2)*P_ms^2;  % Kasen Bildsten Eq 2 (and converted to seconds)  CHECKED

        Cfunct = Ep/tp*(2-1)./(1 + z*tau_m/tp).^2*exp(-x^2).*exp(z.^2).*2.*z;   % Midified Arnett 1982 Eq 31 CHECKED

        C = sum(Cfunct.*(x/Nintegrate));   % Integrate

        L_ph_1(i) = M_ni*m_sol*e_ni0*Lambda1 + M_ni*m_sol*e_co0*Lambda2 + C;   % Calculate luminosity 
  
        if (log10(L_ph_1(i)) > ydata_log_limits(i))

            chi2 = chi2 + 1e7;%*abs(L_ph_1(i) - % TEMP
                              %ydata_log_limits(i));  % Upper bound
                              %violated
            
        end

    end

end

chi2 = chi2/(N_data-4.);  % CHECKED


