analc_endsversion/AnalyticLC/AnalyticLC.m

function [LC,RV_o_r,Y_o_r,Z_o_r,O] = AnalyticLC(P,Tmid0,ror,aor,mu,ex0,ey0,I0,Omega0,t,u1,u2,varargin)
% AnalyticLC: generate an analytic light-curve, and optionally RV and astrometry data, from a set of initial (free) orbital
% elements.

% The full logic of AnalyticLC and the model derivation are described in the paper "An Accurate 3D
% Analytic Model for Exoplanetary Photometry, Radial Velocity and
% Astrometry" by Judkovsky, Ofir and Aharonson (2021), please cite this
% paper upon usage. Model and enjoy :)

%Inputs: 
% P - vector of orbital periods, 
% Tmid0 - vector of refernce times-of-mid-transit, 
% ror - vector of planets/star radii ratio, 
% aor - semi-major-axis in stellar radius of innermost planet, 
% mu - vector of planets/star mass ratio, 
% ex0 and ey0 - vectors of initial FREE eccentricities components of planets
% I0 = vector of initial FREE inclinations (radians), 
% Omega0 - vector of initial FREE longiudes of ascending nodes (radians), 
% t - time vector to calculate light curve model at, 
% u1,u2 - limb darkening coefficients. 
% Eccentricity components and inclination are defined with respect to the system plane, and NOT with respect to the sky plane. The system plane is
% defined such that the x axis points from the star to the observer, and y
% is on the sky plane such that the inclinations are measured with respect
% to the xy plane. Transits are expected to occur for small
% values of sin(I0)*sin(Omega), which is equal to the cosine of the
% sky-plane inclination.
% The time units are to the user's choice, and should be consistent in all
% parameters (t,P,Tmid etc.).

% Examples for usage:
% LC = AnalyticLC(P,Tmid0,ror,aor,mu,ex0,ey0,I0,Omega0,t,u1,u2,'OutputList','');
% LC = AnalyticLC(P,Tmid0,ror,aor,mu,ex0,ey0,I0,Omega0,t,u1,u2,'BinningTime',BinningTime,'OutputList','all');
% [LC,RV_o_r,Y_o_r,Z_o_r,O] = AnalyticLC(P,Tmid0,ror,aor,mu,ex0,ey0,I0,Omega0,t,u1,u2,'tRV',tRV);

    
BinningTime = range(t); %rough time steps for performing the binning to construct the light curve. A value per cadence type is allowed, or a single value for all cadence types. For cadences shorter the the binning time, no binning will be introcued.
CadenceType = ones(size(t)); %cadence type should label each time stamp. values must be consecutive natural integers: 1, 2,...
MaxTTV = inf(1,length(P)); %maximal TTV amplitude allowed for each planet. A value per planet is allowed, or a single value for all planets.
tRV = []; %times to evaluate radial velocity at
ExtraIters = 0; %extra iterations for the TTV calculations
OrdersAll = []; %orders for calculating the forced elements; see function ForcedElements1.
OutputList = 'all'; %a list of parameters to give as an additional output in the structure O. Upon "all", a default list is returned.
Calculate3PlanetsCrossTerms = 0; %Should the 3-planets cross terms be incorporated.
MaxPower = 4; %maximal (joint) power of eccentricities and inclinations used in the expansion
OrbitalElementsAsCellArrays = 0; %an option to return (within the structure O) cell arrays including the orbital elements variations. Each such cell array will contain one cell per planet.

%read optional inputs for the parameters above
for i=1:2:numel(varargin)
    switch lower(varargin{i})
        case 'binningtime'
            BinningTime=varargin{i+1};
        case 'cadencetype'
            CadenceType=varargin{i+1};
        case 'maxttv'
            MaxTTV=varargin{i+1};
        case 'trv'
            tRV=vertical(varargin{i+1});
        case 'extraiters'
            ExtraIters=(varargin{i+1});
        case 'ordersall'
            OrdersAll=(varargin{i+1});
        case 'maxpower'
            MaxPower=(varargin{i+1});
        case 'outputlist'
            OutputList=(varargin{i+1});
        case 'calculate3planetscrossterms'
            Calculate3PlanetsCrossTerms=(varargin{i+1});
        case 'orbitalelementsascellarrays'
            OrbitalElementsAsCellArrays=(varargin{i+1});
        otherwise
            fprintf('WARNING: Parameter %s is unknown and therefore ignored.\n',varargin{i})
    end
end

%number of planets and number of time stamps and number of cadence types
Npl = length(P);
Nt = length(t);
Ncad = max(CadenceType);
Nrv = length(tRV);

%prepare a time axis of one point per transit for the calculation of the
%orbital dynamics
[tT,BeginIdx,EndIdx,Ntr] = DynamicsTimeAxis(P,Tmid0,t);
NtT = length(tT);

%make sure that each cadence type is characterized by its own binning time
if length(BinningTime)==1, BinningTime = BinningTime*ones(1,Ncad); end

%make sure that each planet receives a maximal TTV value
if length(MaxTTV)==1, MaxTTV = MaxTTV*ones(1,Npl); end

%calculate mean motion
n = 2*pi./P;

%calculate time relative to first mid-transit
t0_T = vertical(tT)-Tmid0;
t0_RV = vertical(tRV)-Tmid0;

%calculate mean longitudes
Lambda_T = t0_T.*repmat(n,NtT,1);

%calculate period ratio
PeriodRatio = P(2:end)./P(1);

%calculate semi-major axis of outer planets, deduced from Period Ratio
if length(aor)==1
    aor(2:Npl) = aor(1)*PeriodRatio.^(2/3).*(((1+mu(2:end))/(1+mu(1))).^(1/3));
end

%calculate the secular motion of the eccentricities and inclinations
[A,B] = SecularInteractionsMatrix(n,mu);
free_e_0 = ex0+1i*ey0;
free_I_0 = I0.*exp(1i.*Omega0);
free_e_T = Vector1stOrderDE(A,free_e_0.',horizontal(tT)-t(1)); 
free_I_T = Vector1stOrderDE(B,free_I_0.',horizontal(tT)-t(1));

%calculate the near-resonant interactions and their resulting TTV's
[dz,dLambda,du,da_oa,dTheta,zT,uT,TTV] = ResonantInteractions2(P,mu,Lambda_T,free_e_T,free_I_T,tT,BeginIdx,EndIdx,ExtraIters,OrdersAll,Calculate3PlanetsCrossTerms,MaxPower);


%check if the solution is invalid. Reasons for invalidity:
%if the TTV's do not exceed the maximal value, and if it does - return nans
%innermost planet hits star: a(1-e)<Rs
%e>1 are invalid
if any( max(abs(TTV))>MaxTTV ) || any((1-abs(zT(:,1))).*aor(1)<1) || any(abs(zT(:))>1)
    LC = nan(size(t));
    O = [];
    RV_o_r = nan(size(tRV));
    Y_o_r = nan(size(tRV));
    Z_o_r = nan(size(tRV));
    return
end

%translate the orbital elements to transits properties
CalcDurations = nargout>=5;
[D,D1,Tau,Tau1,w,w1,d,d1,b,b1,AssociatedTmid,AssociatedInd,dbdt,b0] = TransitsProperties(Nt,Npl,BeginIdx,EndIdx,t,tT,TTV,P,n,da_oa,aor,zT,uT,ror,CalcDurations);

%construct the light-curve for each planet
AllLC = GenerateLightCurve(Nt,Npl,Ncad,CadenceType,t,BinningTime,AssociatedTmid,w,d,b,ror,u1,u2);

%sum all planets light-curves to get the total light curve. The possibility of mutual transits is neglected.
LC = sum(AllLC,2)-Npl+1;

%calculate the forced elements of the RV time stamps if the output requires doing so
if Nrv>0 
    
    Lambda_RV = t0_RV.*repmat(n,Nrv,1);
    
    %secular interactions at the RV time stamps
    free_e_RV = Vector1stOrderDE(A,free_e_0.',horizontal(tRV));
    free_I_RV = Vector1stOrderDE(B,free_I_0.',horizontal(tRV));
    
    %near-resonant interactions at the RV time stamps
    [forced_e_RV,dLambda_RV,forced_I_RV] = ForcedElements1(P,Lambda_RV,mu,free_e_RV,free_I_RV,OrdersAll,MaxPower);
    
    zRV = free_e_RV+forced_e_RV;
    uRV = free_I_RV+forced_I_RV;
    Lambda_RV = Lambda_RV+dLambda_RV;
    
    %calculate the radial velocity/astrometry model
    [RV_o_r,Y_o_r,Z_o_r] = CalcRV(Lambda_RV,zRV,uRV,mu,aor,n,nargout>2);
    
else
    
    forced_e_RV = zeros(Nrv,Npl);
    Lambda_RV = zeros(Nrv,Npl);
    free_e_RV = zeros(Nrv,Npl);
    free_I_RV = zeros(Nrv,Npl);
    forced_I_RV = zeros(Nrv,Npl); 
    zRV = free_e_RV+forced_e_RV;
    uRV = free_I_RV+forced_I_RV;
    RV_o_r = [];
    Y_o_r = [];
    Z_o_r  = [];
    
end

%if required, give more outputs then just the light-curve and the RV
%values. These outputs are specified in the optinal input variable
%OutputList.
if nargout>4
    
    %upon the value of "all", set OutputList to include a default list of
    %variables
    if strcmp(OutputList,'all')
        OutputList = {'tT','tRV','zT','zRV','uT','uRV','D','Tau','w','d','b','AssociatedTmid','AssociatedInd','AllLC','TmidActualCell','free_e_T','free_I_T','free_e_RV','free_I_RV','TTVCell','dzCell','duCell','dLambdaCell','zfreeCell','dbdt','b0','Lambda_T','P','BeginIdx','EndIdx','Ntr'}; %a list of parameters to give as an additional output, if required
        OrbitalElementsAsCellArrays = 1;
    end
    
    %generate cell arrays for the orbital elements
    if OrbitalElementsAsCellArrays 
        
        TmidActualCell = cell(1,Npl);
        TTVCell = cell(1,Npl);
        dzCell = cell(1,Npl);
        duCell = cell(1,Npl);
        dLambdaCell = cell(1,Npl);
        zfreeCell = cell(1,Npl);
        ufreeCell = cell(1,Npl);
        da_oaCell = cell(1,Npl);
        
        for j = 1:Npl
            TmidActualCell{j} = (tT(BeginIdx(j):EndIdx(j)))'+TTV(BeginIdx(j):EndIdx(j),j).';
            TTVCell{j} = TTV(BeginIdx(j):EndIdx(j),j).';
            dzCell{j} = dz(BeginIdx(j):EndIdx(j),j).';
            duCell{j} = du(BeginIdx(j):EndIdx(j),j).';
            dLambdaCell{j} = dLambda(BeginIdx(j):EndIdx(j),j).';
            zfreeCell{j} = free_e_T(BeginIdx(j):EndIdx(j),j).';
            ufreeCell{j} = free_I_T(BeginIdx(j):EndIdx(j),j).';
            da_oaCell{j} = da_oa(BeginIdx(j):EndIdx(j),j).';         
        end
        
    end
    
    %move all variables specified in OutputList to a single structure
    %variable O
    for j = 1:length(OutputList)
        try
            eval(['O.',OutputList{j},'=',OutputList{j},';']);
        catch
            fprintf('Requested variable %s does not exist \n',OutputList{j});
        end
    end
    
end


function [AA,BB] = SecularInteractionsMatrix(n,mu)
%calculate the interaction matrix that will describe the equation dz/dt=Az
%where z is a vector of the complex eccentricities and for the similar
%equation for I*exp(i*Omega) using the matrix BB

Npl = length(n);

AA = ones(Npl);
BB = ones(Npl);

for j1 = 1:Npl
    for j2 = j1:Npl
        
        if j1~=j2
            alph = (n(j2)/n(j1))^(2/3)*((1+mu(j1))/(1+mu(j2)))^(1/3);
            A1 = LaplaceCoeffs(alph,1/2,1);
            
            
            DA1 = LaplaceCoeffs(alph,1/2,1,1);
            D2A1 = LaplaceCoeffs(alph,1/2,1,2);
            
            f10 = 0.5*A1-0.5*alph*DA1-0.25*alph^2*D2A1; %calculate f2 and f10 for j=0
            AA(j1,j2) = n(j1)*mu(j2)/(1+mu(j1))*alph*f10;
            AA(j2,j1) = n(j2)*mu(j1)/(1+mu(j2))*f10;
            
            B1 = LaplaceCoeffs(alph,3/2,1);
            f14 = alph*B1;
            BB(j1,j2) = 0.25*n(j1)*mu(j2)/(1+mu(j1))*alph*f14;
            BB(j2,j1) = 0.25*n(j2)*mu(j1)/(1+mu(j2))*f14;
            
        else
            alphvec = zeros(1,Npl);
            alphvec(1:j1) = (n(1:j1)/n(j1)).^(-2/3).*((1+mu(1:j1))./(1+mu(j1))).^(1/3);
            alphvec(j1+1:end) = (n(j1+1:end)/n(j1)).^(2/3).*((1+mu(j1+1:end))./(1+mu(j1))).^(-1/3);
            DA0 = LaplaceCoeffs(alphvec,1/2,0,1);
            D2A0 = LaplaceCoeffs(alphvec,1/2,0,2);
            f2 = 0.25*alphvec.*DA0+1/8*alphvec.^2.*D2A0;
            InnerVec = mu(1:(j1-1))/(1+mu(j1)).*f2(1:(j1-1));
            OuterVec = alphvec((j1+1):end).*mu((j1+1):end)/(1+mu(j1)).*f2((j1+1):end);
            
            AA(j1,j2) = n(j1)*2*(sum(InnerVec)+sum(OuterVec));
            
            B1 = LaplaceCoeffs(alphvec,3/2,1);
            f3 = -0.5*alphvec.*B1;
            InnerVec = mu(1:(j1-1)).*f3(1:(j1-1));
            OuterVec = alphvec((j1+1):end).*mu((j1+1):end).*f3((j1+1):end);
            BB(j1,j2) = 0.25*n(j1)*2*(sum(InnerVec)+sum(OuterVec));
        end
    end
end

AA = 1i*AA;
BB = 1i*BB;


function [TmidActualCell,TTVCell,dzCell,duCell,da_oaCell,dLambdaCell,zfreeCell] = ResonantInteractions1(P,mu,Lambda,free_all,free_i,tT,BeginIdx,EndIdx,ExtraIters,OrdersAll,Calculate3PlanetsCrossTerms,MaxPower)
%calculate the near-resonant interactions arising from the appropriate
%disturbing function terms

%calculate mean motions
n = 2*pi./P;

% [dz,dLambda] = ForcedEccentricitiesAndMeanLongitudes1(P,Lambda,mu,free_all,OrdersAll);
% [dz,dLambda] = ForcedElements(P,Lambda,mu,free_all,free_i,OrdersAll);
[dz,dLambda,du,da_oa] = ForcedElements1(P,Lambda,mu,free_all,free_i,OrdersAll,MaxPower);

if Calculate3PlanetsCrossTerms
    dLambda = dLambda+ForcedMeanLongitudes3PlanetsCrossTerms(P,Lambda,mu,free_all);
end

for jj = 1:ExtraIters
    %     [dz,dLambda] = ForcedEccentricitiesAndMeanLongitudes1(P,Lambda+dLambda,mu,free_all+dz,OrdersAll);
    %     [dz,dLambda] = ForcedElements(P,Lambda+dLambda,mu,free_all+dz,free_i,OrdersAll);
    [dz,dLambda,du,da_oa] = ForcedElements1(P,Lambda+dLambda,mu,free_all+dz,free_i+du,OrdersAll,MaxPower);
   
    %add the 3-planets cross terms
    if Calculate3PlanetsCrossTerms
        dLambda = dLambda+ForcedMeanLongitudes3PlanetsCrossTerms(P,Lambda+dLambda,mu,free_all+dz);
    end
    
end

%calculate the total eccentricity for each nominal transit time
z = free_all+dz;

%calculate the true longitude variations
dTheta = dLambda+2*real((conj(z).*dLambda+conj(dz)/1i).*(1+5/4*conj(dz)));

%translate the variations in true longitude to variations in time via the
%planetary angular velocity
% TTV = -dTheta./repmat(n,size(z,1),1).*((1-abs(z).^2).^(3/2))./((1+real(z)).^2);
TTV = -dTheta./(repmat(n,size(z,1),1).*(1-1.5*da_oa)).*((1-abs(z).^2).^(3/2))./((1+real(z)).^2); %taking into account the variations in n - 8.2.2021

Npl = length(P);
%pre allocation
TmidActualCell = cell(1,Npl);
TTVCell = cell(1,Npl);
dzCell = cell(1,Npl);
dLambdaCell = cell(1,Npl);
zfreeCell = cell(1,Npl);
ufreeCell = cell(1,Npl);
da_oaCell = cell(1,Npl);


for j = 1:Npl
    %     TmidActualCell{j} = TmidActualCell{j}+TTV(BeginIdx(j):EndIdx(j),j).';
    TmidActualCell{j} = (tT(BeginIdx(j):EndIdx(j)))'+TTV(BeginIdx(j):EndIdx(j),j).';
    TTVCell{j} = TTV(BeginIdx(j):EndIdx(j),j).';
    dzCell{j} = dz(BeginIdx(j):EndIdx(j),j).';
    duCell{j} = du(BeginIdx(j):EndIdx(j),j).';
    dLambdaCell{j} = dLambda(BeginIdx(j):EndIdx(j),j).';
    zfreeCell{j} = free_all(BeginIdx(j):EndIdx(j),j).';
    ufreeCell{j} = free_i(BeginIdx(j):EndIdx(j),j).';
    da_oaCell{j} = da_oa(BeginIdx(j):EndIdx(j),j).';
end




function [dz,dLambda,du,da_oa,dTheta,z,u,TTV] = ResonantInteractions2(P,mu,Lambda,free_all,free_i,tT,BeginIdx,EndIdx,ExtraIters,OrdersAll,Calculate3PlanetsCrossTerms,MaxPower)
%calculate the forced orbital elements due to near-resonant interactions

n = 2*pi./P;

%pair-wise interactions
[dz,dLambda,du,da_oa] = ForcedElements1(P,Lambda,mu,free_all,free_i,OrdersAll,MaxPower);

%triplets interactions
if Calculate3PlanetsCrossTerms
    dLambda = dLambda+ForcedMeanLongitudes3PlanetsCrossTerms(P,Lambda,mu,free_all);
end

%iterative solution to the equations, if desired
for jj = 1:ExtraIters    
    [dz,dLambda,du,da_oa] = ForcedElements1(P,Lambda+dLambda,mu,free_all+dz,free_i+du,OrdersAll,MaxPower);
    if Calculate3PlanetsCrossTerms
        dLambda = dLambda+ForcedMeanLongitudes3PlanetsCrossTerms(P,Lambda+dLambda,mu,free_all+dz);
    end
end

%calculate the total eccentricity and inclination for each nominal transit time
z = free_all+dz;
u = free_i+du;

%translate to true longitude
dTheta = dLambda+...
2*real( conj(z).*dLambda+conj(dz)/1i )...
+5/2*real( conj(z).^2.*dLambda+conj(z).*conj(dz).*conj(z)/1i )...
+13/4*real( conj(z).^3.*dLambda+conj(z).^2.*conj(dz)/1i )...
-1/4*real( (conj(z).*dz+z.*conj(dz)).*conj(z)/1i + z.*conj(z).*(conj(z).*dLambda+conj(dz)/1i) )...
+103/24*real( conj(z).^4.*dLambda+conj(z).^3.*conj(dz)/1i )...
-11/24*real( (conj(z).*dz+z.*conj(dz)).*conj(z).^2/1i + 2*z.*conj(z).*(conj(z).^2.*dLambda+conj(z).*conj(dz)/1i) );

%translate the variations in true longitude to variations in time using the
%planetary angular velocity, taking into account the variations in n
TTV = -dTheta./(repmat(n,size(z,1),1).*(1-1.5*da_oa)).*((1-abs(z).^2).^(3/2))./((1+real(z)).^2); 



function [RV_o_r,Y_o_r,Z_o_r] = CalcRV(Lambda,z,u,mu,aor,n,CalcAstrometry)
%calculate the RV at the desired time stamps

Nt = size(z,1);
Pom = angle(z);
Om = angle(u);
e = abs(z);
I = abs(u);
om = Pom-Om;
M = Lambda-Pom;

%calculate sin(f) and cos(f) to 2nd order in e, Solar System Dynamics eq.
%2.84 and 2.85, page 40
sinf = sin(M)+e.*sin(2*M)+e.^2.*(9/8*sin(3*M)-7/8*sin(M));
cosf = cos(M)+e.*(cos(2*M)-1)+9/8*e.^2.*(cos(3*M)-cos(M));

%calculate the planets astrocentric velocities on the x axis
Xdot_o_r = -repmat(n,Nt,1).*repmat(aor,Nt,1)./sqrt(1-e.^2).*((cos(om).*cos(Om)-sin(om).*sin(Om).*cos(I)).*sinf+(sin(om).*cos(Om)+cos(om).*sin(Om).*cos(I)).*(e+cosf));

%sum the astrocentric velocities along x axis to get stellar
%barycentric velocity along x axis
Xsdot_o_r = -sum(repmat(mu,Nt,1).*Xdot_o_r,2)/(1+sum(mu));

%the radial velocity is defined positive when the star moves away from the
%observer, i.e. it is positive then Xsdot is negative and vice versa.
RV_o_r = -Xsdot_o_r;

if CalcAstrometry
    %calculate the planets astrocentric positions on the y and z axes
    y = repmat(aor,Nt,1).*(1-e.^2)./(1+e.*cosf).*(sin(Om).*(cos(om).*cosf-sin(om).*sinf)+cos(Om).*cos(I).*(sin(om).*cosf+cos(om).*sinf));
    z = repmat(aor,Nt,1).*(1-e.^2)./(1+e.*cosf).*sin(I).*(sin(om).*cosf+cos(om).*sinf);
    
    %sum the astrocentric positions along x axis to get the
    %stellar barycentric position along y and z axes
    Y_o_r = -sum(repmat(mu,Nt,1).*y,2)/(1+sum(mu));
    Z_o_r = -sum(repmat(mu,Nt,1).*z,2)/(1+sum(mu));
else
    Y_o_r = [];
    Z_o_r = [];
end


function [td,BeginIdx,EndIdx,Ntr] = DynamicsTimeAxis(P,Tmid0,t)
%generate a time axis that includes one point per transit event
td = [];
BeginIdx = zeros(1,length(P));
EndIdx = zeros(1,length(P));
Ntr = zeros(1,length(P));

for j = 1:length(P)

    td0 = ((min(t)+mod(Tmid0(j)-min(t),P(j))):P(j):max(t))'; %new version, taking into account the possibility that Tmid0 is not within the range of t
    td = [td;td0];
    EndIdx(j) = length(td);
    if j==1, BeginIdx(j) = 1; else, BeginIdx(j) = EndIdx(j-1)+1; end
    Ntr(j) = length(td0);
    
end


function [D,D1,Tau,Tau1,w,w1,d,d1,b,b1,AssociatedTmid,AssociatedInd,dbdt,b0] = TransitsProperties(Nt,Npl,BeginIdx,EndIdx,t,tT,TTV,P,n,da_oa,aor,zT,uT,ror,CalcDurations)
%translate orbital elements at transits to transits parameters

%pre-allocate memory for transit properties: duration, ingress/egress,
%angular velocity, planet-star separation, impact parameter.
D = zeros(Nt,Npl);
Tau = zeros(Nt,Npl);
w = zeros(Nt,Npl);
d = zeros(Nt,Npl);
b = zeros(Nt,Npl);

D1 = cell(1,Npl);
Tau1 = cell(1,Npl);
w1 = cell(1,Npl);
d1 = cell(1,Npl);
b1 = cell(1,Npl);

AssociatedTmid = zeros(Nt,Npl);
AssociatedInd = zeros(Nt,Npl);

dbdt = zeros(1,Npl);
b0 = zeros(1,Npl);


for j = 1:Npl
    
    %indices
    Idx = BeginIdx(j):EndIdx(j);    

    %Associate time-of-mid-transit to each time stamp
    [AssociatedTmid(:,j),AssociatedInd(:,j)] = AssociateTmid(t,tT(Idx)+TTV(Idx,j),P(j));
             
    %based on the presumed orbital elements as a function of time, calculate
    %the transit parameters - Duration, Ingress-Egress, angular velocity,
    %planet-star distance and impact parameter, and translate them to
    %light-curve
    if ~CalcDurations
        [w1{j},d1{j},b1{j}] = OrbitalElements2TransitParams(n(j)*(1-1.5*da_oa(Idx,j)),aor(j)*(1+da_oa(Idx,j)),real(zT(Idx,j)),imag(zT(Idx,j)),abs(uT(Idx,j)),angle(uT(Idx,j)),ror(j)); %treating I, Omega, aor, n as varying
        w(:,j) = w1{j}(AssociatedInd(:,j));
        d(:,j) = d1{j}(AssociatedInd(:,j));
        b(:,j) = b1{j}(AssociatedInd(:,j));
        
    else
        [w1{j},d1{j},b1{j},D1{j},Tau1{j}] = OrbitalElements2TransitParams(n(j)*(1-1.5*da_oa(Idx,j)),aor(j)*(1+da_oa(Idx,j)),real(zT(Idx,j)),imag(zT(Idx,j)),abs(uT(Idx,j)),angle(uT(Idx,j)),ror(j)); %treating I, Omega, aor, n as varying
        w(:,j) = w1{j}(AssociatedInd(:,j));
        d(:,j) = d1{j}(AssociatedInd(:,j));
        b(:,j) = b1{j}(AssociatedInd(:,j));
        D(:,j) = D1{j}(AssociatedInd(:,j));
        Tau(:,j) = Tau1{j}(AssociatedInd(:,j));
        
    end
    
    %calculate db_dt and b0
    [dbdt(j),b0(j)] = fit_linear(tT(Idx),b1{j});    
        
end


function AllLC = GenerateLightCurve(Nt,Npl,Ncad,CadenceType,t,BinningTime,AssociatedTmid,w,d,b,ror,u1,u2)
%construct the light-curve for each planet based on the transit properties
%of each individual event

%pre-allocate
AllLC = zeros(Nt,Npl);

%go over cadence types
for jc = 1:Ncad
    
    %index of data corresponding to the cadence type
    indc = CadenceType==jc;
    
    %construct the integration vector based on half integration time and number
    %of neighbors on each side, which is a function of the binning time
    dtHalfIntegration = median(diff(t(indc)))/2;
    nNeighbors = round(dtHalfIntegration/BinningTime(jc));
    IntegrationStep = 2*dtHalfIntegration/(2*nNeighbors+1);
    IntegrationVector = -dtHalfIntegration+IntegrationStep*(0.5:(2*nNeighbors+1));
    NI = length(IntegrationVector);
    
    %create a table of times around the original times
    t1 = vertical(t(indc))+IntegrationVector;
        
    %go over planets one by one and generate planetary light curve
    for j = 1:Npl
        
        %calculate the phase (for the specific cadence type jc)
        Phi = (t1-repmat(AssociatedTmid(indc,j),1,NI)).*repmat(w(indc,j),1,NI);
        
        %calculate position on the sky with respect to the stellar center in units
        %of stellar radii and using its own w, a, b (for the specific cadence type jc)
        x = repmat(d(indc,j),1,NI).*sin(Phi);
        y = repmat(b(indc,j),1,NI).*cos(Phi);
        
        %calculate the distance from stellar disk center. For phases larger than
        %0.25 null the Mandel-Agol model by using Rsky>1+Rp (for the specific cadence type jc)
        Rsky = sqrt(x.^2+y.^2);
        Rsky(abs(Phi)>0.25) = 5;
        
        %calculate the instanteneous Mandel-Agol model for each time stamp
        %in the binned time vector (for the specific cadence type jc)        
        %lc = occultquadVec(ror(j),Rsky(:),u1,u2);
        lc = MALookupTable(ror(j),Rsky(:),u1,u2);
                
        %bin the light-curve by averaging, and assign in the relevant
        %positions of AllLC
        AllLC(indc,j) = mean(reshape(lc,size(Rsky)),2);
        
    end
    
end