function PrepareMonteCarlo(NAME,N,RedNoise,WhiteNoise,UseJKTEBOP)
load([NAME '.mat'])
if UseJKTEBOP
    load([NAME '/Noise.mat'])
    E.Period= R(1).F_Period;
    E.MidTransit= R(1).F_MidT;
    E.Inclination= R(1).F_Inclination;
    E.LLD= R(1).F_LLD;
    E.QLD= R(1).F_QLD;
    E.ap= (R(1).F_RatioR+1)/(R(1).F_SumR);
    E.p= R(1).F_RatioR;
end
E.NAME= [NAME '_MC'];
E.N=N;
Y=E.YN(1,:);
[phi,F]=quadLimbDark(E.p,E.ap,E.Period,E.Inclination,E.LLD,E.QLD,E.t-E.MidTransit);
Noise=Y-F;
P=Noise*Noise';


M=length(F);
A=nan(M,N);

for i=1:N
    red=rednoise(M);
    white=wgn(1,M,1);
    red=red-mean(red);
    white=white-mean(white);
    red=red /sqrt(red*red');
    white=white/sqrt(white*white');
    total=sqrt(RedNoise)*red+sqrt(WhiteNoise)*white;
    total=total*sqrt(RedNoise+WhiteNoise)*sqrt(P)/sqrt(total*total');
    A(:,i)=F+total;
    
end

E.YN=A;

save([E.NAME '.mat'],'E')   
end

function [phi,F]=quadLimbDark(p,ap,P,i,gamma1,gamma2,t)

% Based on Mandel & Agol's Model for Quadratic Limb Darkening
% Mandel K. & Agol E. 2002, ApJ, 580, L171; please cite this
% paper if you make use of this in your research.  Also, a
% thanks to Gil Nachmani who wrote this routine would be appreciated.

% [phi,F] are the observed phase and relative flux
% p is the planet's radius in units of the star's radius (Rs)
% ap is the planet's orbital radius in units of Rs, assuming zero eccentricity
% P is the planet's period in days
% i is the inclination of the orbit in degrees
% gamma1, gamma2 are the Quadratic Limb Darkening coefficients:
% I(r) = 1-gamma1*(1-mu)-gamma2*(1-mu)^2, where: mu=cos(th)=sqrt(1-r^2);
% E.g. gamma1=0.296, gamma2=0.34 for HD209458
% n is the number of phase points in the resulting lightcurve
% percentOfOrbit is the percentage of orbital phase to be used for the flux
% calculations, e.g. for full orbit use 100(%).

% Gil Nachmani, April 2011

phi=t/P;
Z=ap*(sin(2*pi/P*t).^2+(cos(pi/180*i)*cos(2*pi/P*t)).^2).^(1/2);

c1=0; c2=gamma1+2*gamma2; c3=0; c4=-gamma2; % I(r) = 1-sum(cn*(1-mu^(n/2)))
c0=1-c1-c2-c3-c4; ohmega = c0/(0+4)+c1/(1+4)+c2/(2+4)+c3/(3+4)+c4/(4+4);

j=0; 
for z=Z; j=j+1; 
    a=(z-p)^2; b=(z+p)^2; k=sqrt((1-a)/4/z/p); q=p^2-z^2; 
    k1=acos((1-p^2+z^2)/2/z); k0=acos((p^2+z^2-1)/2/p/z);
       
    %Evaluating lam_e
    if 1+p<z || abs(phi(j))>(p+1)/ap/2/pi;
        lam_e = 0;
    elseif abs(1-p)<z && z<=1+p;
        lam_e = 1/pi*(p^2*k0+k1-1/2*sqrt(4*z^2-(1+z^2-p^2)^2));
    elseif z<=1-p && z>p-1;
        lam_e = p^2;
    elseif z<=p-1;
        lam_e = 1;
    end
    
    %Evaluating lam_d and eta_d
    if z>=1+p || p==0 || abs(phi(j))>(p+1)/ap/2/pi; %1
        lam_d = 0; eta_d=0; 
    elseif z>=1/2+abs(p-1/2) && z<1+p; %2
        lam_d = lam1(p,z,a,b,k,q); eta_d = eta1(p,z,a,b,k1,k0);
    elseif p<0.5 && z>p && z<1-p; %3
        lam_d = lam2(p,z,a,b,k,q); eta_d = eta2(p,z);
    elseif p<0.5 && z==1-p; %4
        lam_d = lam5(p); eta_d = eta2(p,z);
    elseif p<0.5 && z==p; %5
        lam_d = lam4(p); eta_d = eta2(p,z);
    elseif p==0.5 && z==0.5; %6
        lam_d = 1/3-4/pi/9; eta_d = 3/32; 
    elseif p>0.5 && z==p; %7
        lam_d = lam3(p); eta_d = eta1(p,z,a,b,k1,k0); 
    elseif p>0.5 && z>=abs(1-p) && z<p; %8
        lam_d = lam1(p,z,a,b,k,q); eta_d = eta1(p,z,a,b,k1,k0); 
    elseif p<1 && z>0 && z<=1/2-abs(p-1/2); %9
        lam_d = lam2(p,z,a,b,k,q); eta_d = eta2(p,z); 
    elseif p<1 && z==0; %10
        lam_d = lam6(p); eta_d=eta2(p,z); 
    elseif p>1 && z<=p-1; %11
        lam_d = 0; eta_d=1/2; 
    end
F(j) = 1 - 1/(4*ohmega)*( (1-c2)*lam_e + c2*(lam_d+2/3*heaviside(p-z)) - c4*eta_d );
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function lam=lam1(p,z,a,b,k,q)
lam = 1/9/pi/sqrt(p*z) * ( ((1-b)*(2*b+a-3)-3*q*(b-2))*K(k) + 4*p*z*(z^2+7*p^2-4)*E(k)-3*q/a*PI((a-1)/a,k) );
end
function lam=lam2(p,z,a,b,k,q)
lam = 2/9/pi/sqrt(1-a) * ( (1-5*z^2+p^2+q^2)*K(1/k) + (1-a)*(z^2+7*p^2-4)*E(1/k)-3*q/a*PI((a-b)/a,1/k) );
end
function lam=lam3(p)
lam = 1/3 + 16*p/9/pi*(2*p^2-1)*E(1/2/p) - (1-4*p^2)*(3-8*p^2)/9/pi/p*K(1/2/p);
end
function lam=lam4(p)
lam = 1/3 + 2/9/pi*(4*(2*p^2-1)*E(2*p)+(1-4*p^2)*K(2*p));
end
function lam=lam5(p)
lam = 2/3/pi*acos(1-2*p) - 4/9/pi*(3+2*p-8*p^2)*sqrt(p*(1-p))-2/3*heaviside(p-1/2);
end
function lam=lam6(p)
lam = -2/3*(1-p^2)^(3/2);
end
function eta=eta1(p,z,a,b,k1,k0)
eta = 1/2/pi*(k1+2*eta2(p,z)*k0-1/4*(1+5*p^2+z^2)*sqrt((1-a)*(b-1)));
end
function eta=eta2(p,z)
eta = p^2/2*(p^2+2*z^2);
end
function f=E(k)
f = mfun('EllipticE',k);
end
function f=PI(n,k) 
f = mfun('EllipticPi',n,k);
end
function f=K(k)
f = mfun('EllipticK',k);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Red Noise Generation with MATLAB Implementation    %
%                                                      %
% Author: M.Sc. Eng. Hristo Zhivomirov       07/31/13  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y = rednoise(N)

% function: y = rednoise(N) 
% N - number of samples to be returned in row vector
% y - row vector of red (brown) noise samples

% The function generates a sequence of red (brown) noise samples. 
% In terms of power at a constant bandwidth, red noise falls off at 6 dB per octave. 

% difine the length of the vector
% ensure that the M is even
if rem(N,2)
    M = N+1;
else
    M = N;
end

% generate white noise
x = randn(1, M);

% FFT
X = fft(x);

% prepare a vector for 1/(f^2) multiplication
NumUniquePts = M/2 + 1;
n = 1:NumUniquePts;

% multiplicate the left half of the spectrum so the power spectral density
% is proportional to the frequency by factor 1/(f^2), i.e. the
% amplitudes are proportional to 1/f
X(1:NumUniquePts) = X(1:NumUniquePts)./n;

% prepare a right half of the spectrum - a copy of the left one,
% except the DC component and Nyquist frequency - they are unique
X(NumUniquePts+1:M) = real(X(M/2:-1:2)) -1i*imag(X(M/2:-1:2));

% IFFT
y = ifft(X);

% prepare output vector y
y = real(y(1, 1:N));

% ensure unity standard deviation and zero mean value
y = y - mean(y);
yrms = sqrt(mean(y.^2));
y = y/yrms;

end


