commit 8e602ca6bf4e1557a06179c88eeb8b28e60a8f80 Author: Tim Blume Date: Mon Apr 11 15:28:22 2022 +0200 initial version diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5ceb386 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +venv diff --git a/AnalyticLC/AnalyticLC.m b/AnalyticLC/AnalyticLC.m new file mode 100644 index 0000000..c5846c7 --- /dev/null +++ b/AnalyticLC/AnalyticLC.m @@ -0,0 +1,556 @@ +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)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 \ No newline at end of file diff --git a/AnalyticLC/AssociateTmid.m b/AnalyticLC/AssociateTmid.m new file mode 100644 index 0000000..289b3ac --- /dev/null +++ b/AnalyticLC/AssociateTmid.m @@ -0,0 +1,44 @@ +function [AssociatedTmid,AssociatedInd] = AssociateTmid(t,Tmid,P) +% AssociateTmid: Associate a time-of-mid-transit to each time stamp t. + +%Inputs: t - list of times +% Tmid - list of times of mid transit +% P - orbital period (optional input) + +% Yair Judkovsky, 2.5.2020 + +%estimate the period if it is not given +if ~exist('P','var') +P = median(diff(Tmid)); +end + +%transit indices - the first transit in the data is one, and going up +try + Tr_Ind = 1+round((Tmid-Tmid(1))/P); +catch + AssociatedTmid = nan(size(t)); + AssociatedInd = nan(size(t)); + return +end + +%associate an index of transit per time stamp. Any event before the first +%Tmid will be associated with the first Tmid. Any event after the last Tmid +%will be associated with the last Tmid. +AssociatedInd = 1+round((t-Tmid(1))/P); +AssociatedInd(AssociatedInd<1) = 1; +AssociatedInd(AssociatedInd>Tr_Ind(end)) = Tr_Ind(end); + +%Associated transit index +try +Tmid_Of_Tr_Ind = inf(1,max(Tr_Ind)); +Tmid_Of_Tr_Ind(Tr_Ind) = Tmid; +%next few lines for debugging purposes +catch + save('debug_tmid'); + AssociatedTmid = nan(size(t)); + AssociatedInd = nan(size(t)); + return +end + +%associate a value of Tmid for each time stamp +AssociatedTmid = Tmid_Of_Tr_Ind(AssociatedInd); \ No newline at end of file diff --git a/AnalyticLC/DisturbingFunctionMultiTermVariations.m b/AnalyticLC/DisturbingFunctionMultiTermVariations.m new file mode 100644 index 0000000..4ffa007 --- /dev/null +++ b/AnalyticLC/DisturbingFunctionMultiTermVariations.m @@ -0,0 +1,154 @@ +function [dLambda1,dLambda2,dz1,dz2,du1,du2,da_oa1,da_oa2] = DisturbingFunctionMultiTermVariations(n1,n2,alph,Lambda1,Lambda2,mu1,mu2,j,k,e1,e2,Pom1,Pom2,I1,I2,s1,s2,Om1,Om2,A,Ap,B,Bp,C,Cp,D,Dp,f,fE,fI,Df,DfE,DfI) + +%get the number of time stamps and number of terms, in order to replicate +%the vectors and perform the calculation once for all terms +% Ntime = length(Lambda1); +% Nterm = length(k); +% +% +% Lambda1 = repmat(Lambda1,1,Nterm); +% Lambda2 = repmat(Lambda2,1,Nterm); +% e1 = repmat(e1,1,Nterm); +% e2 = repmat(e2,1,Nterm); +% Pom1 = repmat(Pom1,1,Nterm); +% Pom2 = repmat(Pom2,1,Nterm); +% I1 = repmat(I1,1,Nterm); +% I2 = repmat(I2,1,Nterm); +% s1 = repmat(s1,1,Nterm); +% s2 = repmat(s2,1,Nterm); +% Om1 = repmat(Om1,1,Nterm); +% Om2 = repmat(Om2,1,Nterm); +% +% if length(j)==length(k) +% j = repmat(j,Ntime,1); +% end +% k = repmat(k,Ntime,1); +% A = repmat(A,Ntime,1); +% Ap = repmat(Ap,Ntime,1); +% B = repmat(B,Ntime,1); +% Bp = repmat(Bp,Ntime,1); +% C = repmat(C,Ntime,1); +% Cp = repmat(Cp,Ntime,1); +% D = repmat(D,Ntime,1); +% Dp = repmat(Dp,Ntime,1); + +f_inner = f+fE; +f_outer = f+fI; +Df_inner = Df+DfE; +Df_outer = Df+DfI; + +% f = repmat(f,Ntime,1); +% fE = repmat(fE,Ntime,1); +% fI = repmat(fI,Ntime,1); +% Df = repmat(Df,Ntime,1); +% DfE = repmat(DfE,Ntime,1); +% DfI = repmat(DfI,Ntime,1); + +% f_inner = repmat(f_inner,Ntime,1); +% f_outer = repmat(f_outer,Ntime,1); +% Df_inner = repmat(Df_inner,Ntime,1); +% Df_outer = repmat(Df_outer,Ntime,1); + +SQ1 = sqrt(1-e1.^2); +SQ2 = sqrt(1-e2.^2); + +%calculate the cosine argument +Phi = j.*Lambda2+(k-j).*Lambda1-C.*Pom1-Cp.*Pom2-D.*Om1-Dp.*Om2; + +%calculate Cjk and Sjk +PowersFactor = e1.^A.*e2.^Ap.*s1.^B.*s2.^Bp; +Cjk = PowersFactor.*cos(Phi); +Sjk = PowersFactor.*sin(Phi); + +%calculate the mean-motion of the longitude of conjunction +njk = j*n2+(k-j)*n1; + +%orbital elements variations - inner planet +da_oa1 = 2*mu2/(1+mu1)*alph*(f_inner).*(k-j)*n1./njk.*Cjk; + +dLambda1 = mu2/(1+mu1)*alph*n1./njk.*Sjk... + .*(... + 3.*(f_inner).*(j-k)*n1./njk... + -2*alph*(Df_inner)... + +(f_inner).*A.*(SQ1.*(1-SQ1))./(e1.^2)... + +(f_inner).*B/2./SQ1... +); + +de1 = mu2/(1+mu1)*alph*(f_inner).*SQ1./e1*n1./njk.*Cjk... + .*((j-k).*(1-SQ1)+C); + +dPom1 = mu2/(1+mu1)*alph*(f_inner)*n1./njk.*Sjk... + .*(A.*SQ1./(e1.^2)+B/2./SQ1); + +dI1 = mu2/(1+mu1)*alph*(f_inner)./SQ1*n1./njk.*Cjk... + .*((j-k+C).*tan(I1/2)+D./sin(I1)); + +dOm1 = mu2/(1+mu1)*alph*(f_inner).*B/2.*cot(I1/2)./(SQ1.*sin(I1))*n1./njk.*Sjk; + +% dI1 = 0; dOm1 = 0; + + + +%orbital elements variations - outer planet +da_oa2 = 2*mu1/(1+mu2)*(f_outer).*j*n2./njk.*Cjk; + +dLambda2 = mu1/(1+mu2)*n2./njk.*Sjk.*... + (... + -3*(f_outer).*j*n2./njk... + +2*(f_outer+alph*Df_outer)... + +(f_outer).*Ap.*SQ2.*(1-SQ2)./(e2.^2)... + +(f_outer).*Bp/2./SQ2... +); + +de2 = mu1/(1+mu2)*(f_outer).*SQ2./e2*n2./njk.*Cjk.*... + (... + -j.*(1-SQ2)+Cp... +); + +dPom2 = mu1/(1+mu2)*(f_outer)*n2./njk.*Sjk.*... + (... + Ap.*SQ2./(e2.^2)+Bp/2./SQ2... +); + +dI2 = mu1/(1+mu2)*(f_outer)./SQ2*n2./njk.*Cjk.*... + (... + (Cp-j).*tan(I2/2)+Dp./sin(I2)... +); + +dOm2 = mu1/(1+mu2)*(f_outer).*Bp/2.*cot(I2/2)./(SQ2.*sin(I2))*n2./njk.*Sjk; + +% dI2 = 0; dOm2 = 0; + +NewVersion = 1; + + +if ~NewVersion +%translate to the complex form +dz1 = (de1+1i*dPom1.*e1).*exp(1i*Pom1); +dz2 = (de2+1i*dPom2.*e2).*exp(1i*Pom2); + +du1 = (dI1+1i*dOm1.*I1).*exp(1i*Om1); +du2 = (dI2+1i*dOm2.*I2).*exp(1i*Om2); + + +%sum the distribution of all terms to get the total effect +dz1 = sum(dz1,2); +dz2 = sum(dz2,2); +du1 = sum(du1,2); +du2 = sum(du2,2); +da_oa1 = sum(da_oa1,2); +da_oa2 = sum(da_oa2,2); +dLambda1 = sum(dLambda1,2); +dLambda2 = sum(dLambda2,2); + +else + %translate to the complex form +dz1 = (sum(de1,2)+1i*sum(dPom1,2).*e1).*exp(1i*Pom1); +dz2 = (sum(de2,2)+1i*sum(dPom2,2).*e2).*exp(1i*Pom2); +du1 = (sum(dI1,2)+1i*sum(dOm1,2).*I1).*exp(1i*Om1); +du2 = (sum(dI2,2)+1i*sum(dOm2,2).*I2).*exp(1i*Om2); +da_oa1 = sum(da_oa1,2); +da_oa2 = sum(da_oa2,2); +dLambda1 = sum(dLambda1,2); +dLambda2 = sum(dLambda2,2); +end \ No newline at end of file diff --git a/AnalyticLC/ForcedElements1.m b/AnalyticLC/ForcedElements1.m new file mode 100644 index 0000000..e24ed7e --- /dev/null +++ b/AnalyticLC/ForcedElements1.m @@ -0,0 +1,1230 @@ +function [dz,dLambda,du,da_oa] = ForcedElements1(P,Lambda,mu,z,u,OrdersAll,MaxPower) +% ForcedElements1: calculate the forced eccentricities and mean longitudes +% by a disturbing function expansion. + +%Difference from ForcedElements: +% rather than writing down a series of terms, use a general expression for +% the disturbing function terms contributions and sum them up. + +%Inputs: P - orbital periods, Lambda - mean longitudes, mu - planets/star +%mass ratio, z - free eccentricities (if they are not given, zero values +%are assumed), u - free inclinations (=I*exp(i*Omega) (if they are not given, zero values +%are assumed). Labmda, z and u can be given as columns of multiple values. +%OrdersAll: the number of j values taken in the expansion. If given an +%empty value [], the code will decide by itself what is the required +%number. MaxPower: the maximal (joint) power of eccentricities and +%inclinations taken in the expansion. Default is 4. + +% Yair Judkovsky, 6.12.2020. Coded with lots of blood, swear and ink :) + +if ~exist('z','var') + z = zeros(size(P)); +end + +if ~exist('u','var') + u = zeros(size(P)); +end + +if ~exist('OrdersAll','var') + OrdersAll = []; +end + +if ~exist('MaxPower','var') + MaxPower = 4; +end + +%mean motions +n = 2*pi./P; + +%absolute values and angles of eccentricities +e = abs(z); +pom = angle(z); + +%absolute values and angles of inclinations +I = abs(u); +Om = angle(u); + +%allocate space for 3d matrices that will hold the values of dLambda and dz +%according to these dimensios: orders, data points, number of planets +dLambdaMatrix = zeros(size(z,1),length(P)); +dzMatrix = zeros(size(z,1),length(P)); +duMatrix = zeros(size(z,1),length(P)); +da_oaMatrix = zeros(size(z,1),length(P)); + +%for each pair of planets, calculate the 2nd order forced +%eccentricity. secular terms are ignored. +for j1 = 1:length(P) + for j2 = (j1+1):length(P) + + %ratio of semi-major axes + alph = (n(j2)/n(j1))^(2/3)*((1+mu(j1))/(1+mu(j2)))^(1/3); + + %write the orbital elements for this specific pair + Lambda1 = Lambda(:,j1); + Lambda2 = Lambda(:,j2); + + %write the complex eccentricities, magnitude and angle of the + %currently analyzed pair + e1 = e(:,j1); + e2 = e(:,j2); + pom1 = pom(:,j1); + pom2 = pom(:,j2); + + %write the complex inclinations, magnitude and angle of the + %currently analyzed pair + I1 = I(:,j1); + I2 = I(:,j2); + Om1 = Om(:,j1); + Om2 = Om(:,j2); + + %translate from I to s - as defined in the disturbing function + %expansion (Solar System Dynamics) + s1 = sin(I1/2); + s2 = sin(I2/2); + + %write the relative masses and the mean motions of the currently analyzed pair + mu1 = mu(j1); + mu2 = mu(j2); + n1 = n(j1); + n2 = n(j2); + + %set number of orders, if not given, by the proximity to closest + %MMRs + if isempty(OrdersAll) + J = round((1:MaxPower)*P(j2)/(P(j2)-P(j1))); + Orders = max([J 6])+1; + else + Orders = OrdersAll; + end + + %prepare the structure Terms, which will include the table elements + Terms = []; + clear AllTerms + + %here, instead of the loop on j, make a loop on the different + %terms, to be continued + for j = -Orders:Orders + + + %% CalculateSMAFunctions; + + %Calculate Laplace coefficients and their semi-major axis functions - Solar + %System Dynamics (Murray & Dermott 1999), Appendix B + + for DummyIdx = 1:1 + + + +% Aj = LaplaceCoeffs(alph,0.5,j,0); +% DAj = LaplaceCoeffs(alph,0.5,j,1); +% D2Aj = LaplaceCoeffs(alph,0.5,j,2); +% D3Aj = LaplaceCoeffs(alph,0.5,j,3); +% D4Aj = LaplaceCoeffs(alph,0.5,j,4); +% D5Aj = LaplaceCoeffs(alph,0.5,j,5); +% +% Ajp1 = LaplaceCoeffs(alph,0.5,j+1,0); +% DAjp1 = LaplaceCoeffs(alph,0.5,j+1,1); +% D2Ajp1 = LaplaceCoeffs(alph,0.5,j+1,2); +% D3Ajp1 = LaplaceCoeffs(alph,0.5,j+1,3); +% D4Ajp1 = LaplaceCoeffs(alph,0.5,j+1,4); +% D5Ajp1 = LaplaceCoeffs(alph,0.5,j+1,5); +% +% Ajp2 = LaplaceCoeffs(alph,0.5,j+2,0); +% DAjp2 = LaplaceCoeffs(alph,0.5,j+2,1); +% D2Ajp2 = LaplaceCoeffs(alph,0.5,j+2,2); +% D3Ajp2 = LaplaceCoeffs(alph,0.5,j+2,3); +% D4Ajp2 = LaplaceCoeffs(alph,0.5,j+2,4); +% D5Ajp2 = LaplaceCoeffs(alph,0.5,j+2,5); +% +% Ajm1 = LaplaceCoeffs(alph,0.5,j-1,0); +% DAjm1 = LaplaceCoeffs(alph,0.5,j-1,1); +% D2Ajm1 = LaplaceCoeffs(alph,0.5,j-1,2); +% D3Ajm1 = LaplaceCoeffs(alph,0.5,j-1,3); +% D4Ajm1 = LaplaceCoeffs(alph,0.5,j-1,4); +% D5Ajm1 = LaplaceCoeffs(alph,0.5,j-1,5); +% +% Ajm2 = LaplaceCoeffs(alph,0.5,j-2,0); +% DAjm2 = LaplaceCoeffs(alph,0.5,j-2,1); +% D2Ajm2 = LaplaceCoeffs(alph,0.5,j-2,2); +% D3Ajm2 = LaplaceCoeffs(alph,0.5,j-2,3); +% D4Ajm2 = LaplaceCoeffs(alph,0.5,j-2,4); +% D5Ajm2 = LaplaceCoeffs(alph,0.5,j-2,5); +% +% Ajm3 = LaplaceCoeffs(alph,0.5,j-3,0); +% DAjm3 = LaplaceCoeffs(alph,0.5,j-3,1); +% D2Ajm3 = LaplaceCoeffs(alph,0.5,j-3,2); +% D3Ajm3 = LaplaceCoeffs(alph,0.5,j-3,3); +% D4Ajm3 = LaplaceCoeffs(alph,0.5,j-3,4); +% D5Ajm3 = LaplaceCoeffs(alph,0.5,j-3,5); +% +% Ajm4 = LaplaceCoeffs(alph,0.5,j-4,0); +% DAjm4 = LaplaceCoeffs(alph,0.5,j-4,1); +% D2Ajm4 = LaplaceCoeffs(alph,0.5,j-4,2); +% D3Ajm4 = LaplaceCoeffs(alph,0.5,j-4,3); +% D4Ajm4 = LaplaceCoeffs(alph,0.5,j-4,4); +% D5Ajm4 = LaplaceCoeffs(alph,0.5,j-4,5); +% +% +% Bj = LaplaceCoeffs(alph,1.5,j,0); +% DBj = LaplaceCoeffs(alph,1.5,j,1); +% D2Bj = LaplaceCoeffs(alph,1.5,j,2); +% D3Bj = LaplaceCoeffs(alph,1.5,j,3); +% +% Bjm1 = LaplaceCoeffs(alph,1.5,j-1,0); +% DBjm1 = LaplaceCoeffs(alph,1.5,j-1,1); +% D2Bjm1 = LaplaceCoeffs(alph,1.5,j-1,2); +% D3Bjm1 = LaplaceCoeffs(alph,1.5,j-1,3); +% +% Bjm2 = LaplaceCoeffs(alph,1.5,j-2,0); +% DBjm2 = LaplaceCoeffs(alph,1.5,j-2,1); +% D2Bjm2 = LaplaceCoeffs(alph,1.5,j-2,2); +% D3Bjm2 = LaplaceCoeffs(alph,1.5,j-2,3); +% +% Bjm3 = LaplaceCoeffs(alph,1.5,j-3,0); +% DBjm3 = LaplaceCoeffs(alph,1.5,j-3,1); +% D2Bjm3 = LaplaceCoeffs(alph,1.5,j-3,2); +% D3Bjm3 = LaplaceCoeffs(alph,1.5,j-3,3); +% +% +% +% Bjp1 = LaplaceCoeffs(alph,1.5,j+1,0); +% DBjp1 = LaplaceCoeffs(alph,1.5,j+1,1); +% D2Bjp1 = LaplaceCoeffs(alph,1.5,j+1,2); +% D3Bjp1 = LaplaceCoeffs(alph,1.5,j+1,3); +% +% Bjp2 = LaplaceCoeffs(alph,1.5,j+2,0); +% DBjp2 = LaplaceCoeffs(alph,1.5,j+2,1); +% D2Bjp2 = LaplaceCoeffs(alph,1.5,j+2,2); +% D3Bjp2 = LaplaceCoeffs(alph,1.5,j+2,3); +% +% Cj = LaplaceCoeffs(alph,2.5,j,0); +% DCj = LaplaceCoeffs(alph,2.5,j,1); +% +% Cjm2 = LaplaceCoeffs(alph,2.5,j-2,0); +% DCjm2 = LaplaceCoeffs(alph,2.5,j-2,1); +% +% Cjp2 = LaplaceCoeffs(alph,2.5,j+2,0); +% DCjp2 = LaplaceCoeffs(alph,2.5,j+2,1); + + + + sVec = [0.5*ones(1,42) 1.5*ones(1,24) 2.5*ones(1,6)]; + jVec = [j*ones(1,6) (j+1)*ones(1,6) (j+2)*ones(1,6) (j-1)*ones(1,6) (j-2)*ones(1,6) (j-3)*ones(1,6) (j-4)*ones(1,6) (j)*ones(1,4) (j-1)*ones(1,4) (j-2)*ones(1,4) (j-3)*ones(1,4) (j+1)*ones(1,4) (j+2)*ones(1,4) j j j-2 j-2 j+2 j+2]; + DVec = [0:5 0:5 0:5 0:5 0:5 0:5 0:5 0:3 0:3 0:3 0:3 0:3 0:3 0:1 0:1 0:1]; + CoeffsVec = LaplaceCoeffs(alph,sVec,jVec,DVec); + Aj = CoeffsVec(1); + DAj = CoeffsVec(2); + D2Aj = CoeffsVec(3); + D3Aj = CoeffsVec(4); + D4Aj = CoeffsVec(5); + D5Aj = CoeffsVec(6); + + Ajp1 = CoeffsVec(7); + DAjp1 = CoeffsVec(8); + D2Ajp1 = CoeffsVec(9); + D3Ajp1 = CoeffsVec(10); + D4Ajp1 = CoeffsVec(11); + D5Ajp1 = CoeffsVec(12); + + Ajp2 = CoeffsVec(13); + DAjp2 = CoeffsVec(14); + D2Ajp2 = CoeffsVec(15); + D3Ajp2 = CoeffsVec(16); + D4Ajp2 = CoeffsVec(17); + D5Ajp2 = CoeffsVec(18); + + Ajm1 = CoeffsVec(19); + DAjm1 = CoeffsVec(20); + D2Ajm1 = CoeffsVec(21); + D3Ajm1 = CoeffsVec(22); + D4Ajm1 = CoeffsVec(23); + D5Ajm1 = CoeffsVec(24); + + Ajm2 = CoeffsVec(25); + DAjm2 = CoeffsVec(26); + D2Ajm2 = CoeffsVec(27); + D3Ajm2 = CoeffsVec(28); + D4Ajm2 = CoeffsVec(29); + D5Ajm2 = CoeffsVec(30); + + Ajm3 = CoeffsVec(31); + DAjm3 = CoeffsVec(32); + D2Ajm3 = CoeffsVec(33); + D3Ajm3 = CoeffsVec(34); + D4Ajm3 = CoeffsVec(35); + D5Ajm3 = CoeffsVec(36); + + Ajm4 = CoeffsVec(37); + DAjm4 = CoeffsVec(38); + D2Ajm4 = CoeffsVec(39); + D3Ajm4 = CoeffsVec(40); + D4Ajm4 = CoeffsVec(41); + D5Ajm4 = CoeffsVec(42); + + + Bj = CoeffsVec(43); + DBj = CoeffsVec(44); + D2Bj = CoeffsVec(45); + D3Bj = CoeffsVec(46); + + Bjm1 = CoeffsVec(47); + DBjm1 = CoeffsVec(48); + D2Bjm1 = CoeffsVec(49); + D3Bjm1 = CoeffsVec(50); + + Bjm2 = CoeffsVec(51); + DBjm2 = CoeffsVec(52); + D2Bjm2 = CoeffsVec(53); + D3Bjm2 = CoeffsVec(54); + + Bjm3 = CoeffsVec(55); + DBjm3 = CoeffsVec(56); + D2Bjm3 = CoeffsVec(57); + D3Bjm3 = CoeffsVec(58); + + + + Bjp1 = CoeffsVec(59); + DBjp1 = CoeffsVec(60); + D2Bjp1 = CoeffsVec(61); + D3Bjp1 = CoeffsVec(62); + + Bjp2 = CoeffsVec(63); + DBjp2 = CoeffsVec(64); + D2Bjp2 = CoeffsVec(65); + D3Bjp2 = CoeffsVec(66); + + Cj = CoeffsVec(67); + DCj = CoeffsVec(68); + + Cjm2 = CoeffsVec(69); + DCjm2 = CoeffsVec(70); + + Cjp2 = CoeffsVec(71); + DCjp2 = CoeffsVec(72); + + + + + + %0th order semi-major axis functions + f1 = 0.5*Aj; + Df1 = 0.5*DAj; + + f2 = 1/8*(-4*j^2*Aj +2*alph*DAj +alph^2*D2Aj); + Df2 = 1/8*(-4*j^2*DAj+2*DAj+2*alph*D2Aj+2*alph*D2Aj+alph^2*D3Aj); + + f3 = -0.25*alph*Bjm1 -0.25*alph*Bjp1; + Df3 = -0.25*Bjm1-0.25*alph*DBjm1-0.25*Bjp1-0.25*alph*DBjp1; + + f4 = 1/128*((-9*j^2+16*j^4)*Aj +(-8*j^2)*alph*DAj +(-8*j^2)*alph^2*D2Aj +4*alph^3*D3Aj +alph^4*D4Aj); + Df4 = 1/128*((-9*j^2+16*j^4)*DAj+(-8*j^2)*(DAj+alph*D2Aj)+(-8*j^2)*(2*alph*D2Aj+alph^2*D3Aj)+4*(3*alph^2*D3Aj+alph^3*D4Aj)+(4*alph^3*D4Aj+alph^4+D5Aj)); + + f5 = 1/32*((16*j^4)*Aj +(4-16*j^2)*alph*DAj +(14-8*j^2)*alph^2*D2Aj +(8)*alph^3*D3Aj +(1)*alph^4*D4Aj); + Df5 = 1/32*((16*j^4)*DAj +(4-16*j^2)*(DAj+alph*D2Aj) +(14-8*j^2)*(2*alph*D2Aj+alph^2*D3Aj) +(8)*(3*alph^2*D3Aj+alph^3*D4Aj) +(1)*(4*alph^3*D4Aj+alph^4*D5Aj)); + + f6 = 1/128*((-17*j^2+16*j^4)*Aj +(24-24*j^2)*alph*DAj +(36-8*j^2)*alph^2*D2Aj +(12)*alph^3*D3Aj +(1)*alph^4*D4Aj); + Df6 = 1/128*((-17*j^2+16*j^4)*DAj+(24-24*j^2)*(DAj+alph*D2Aj)+(36-8*j^2)*(2*alph*D2Aj+alph^2*D3Aj)+(12)*(3*alph^2*D3Aj+alph^3*D4Aj)+(1)*(4*alph^3*D4Aj+alph^4*D5Aj)); + + f7 = 1/16*((-2+4*j^2)*alph*Bjm1 +(-4)*alph^2*DBjm1 +(-1)*alph^3*D2Bjm1)... + +1/16*((-2+4*j^2)*alph*Bjp1 +(-4)*alph^2*DBjp1 +(-1)*alph^3*D2Bjp1); + Df7 = 1/16*((-2+4*j^2)*(Bjm1+alph*DBjm1)+(-4)*(2*alph*DBjm1+alph^2*D2Bjm1)+(-1)*(3*alph^2*D2Bjm1+alph^3*D3Bjm1))... + +1/16*((-2+4*j^2)*(Bjp1+alph*DBjp1) +(-4)*(2*alph*DBjp1+alph^2*D2Bjp1) +(-1)*(3*alph^2*D2Bjp1+alph^3*D3Bjp1)); + + f8 = alph^2*(3/16*Cjm2+0.75*Cj+3/16*Cjp2); + Df8 = 2*alph*(3/16*Cjm2+0.75*Cj+3/16*Cjp2)+alph^2*(3/16*DCjm2+0.75*DCj+3/16*DCjp2); + + f9 = 0.25*alph*(Bjm1+Bjp1) +3/8*alph^2*Cjm2 +15/4*alph^2*Cj +3/8*alph^2*Cjp2; + Df9 = 0.25*((Bjm1+Bjp1)+alph*(DBjm1+DBjp1))+3/8*(2*alph*Cjm2+alph^2*DCjm2)+15/4*(2*alph*Cj+alph^2*DCj)+3/8*(2*alph*Cjp2+alph^2*DCjp2); + + f10 = 0.25*((2+6*j+4*j^2)*Ajp1 -2*alph*DAjp1 -alph^2*D2Ajp1); + Df10 = 0.25*((2+6*j+4*j^2)*DAjp1-2*(DAjp1+alph*D2Ajp1)-(2*alph*D2Ajp1+alph^2*D3Ajp1)); + + f11 = 1/32*((-6*j-26*j^2-36*j^3-16*j^4)*Ajp1 +(6*j+12*j^2)*alph*DAjp1 +(-4+7*j)*alph^2*D2Ajp1 -6*alph^3*D3Ajp1 -alph^4*D4Ajp1); + Df11 = 1/32*((-6*j-26*j^2-36*j^3-16*j^4)*DAjp1+(6*j+12*j^2)*(DAjp1+alph*D2Ajp1)+(-4+7*j)*(2*alph*D2Ajp1+alph^2*D3Ajp1)-6*(3*alph^2*D3Ajp1+alph^3*D4Ajp1)-(4*alph^3*D4Ajp1+alph^4*D5Ajp1)); + + f12 = 1/32*((4+2*j-22*j^2-36*j^3-16*j^4)*Ajp1 +(-4+22*j+20*j^2)*alph*DAjp1 +(-22+7*j+8*j^2)*alph^2*D2Ajp1 -10*alph^3*D3Ajp1 -alph^4*D4Ajp1); + Df12 = 1/32*((4+2*j-22*j^2-36*j^3-16*j^4)*DAjp1+(-4+22*j+20*j^2)*(DAjp1+alph*D2Ajp1)+(-22+7*j+8*j^2)*(2*alph*D2Ajp1+alph^2*D3Ajp1)-10*(3*alph^2*D3Ajp1+alph^3*D4Ajp1)-(4*alph^3*D4Ajp1+alph^4*D5Ajp1)); + + f13 = 1/8*((-6*j-4*j^2)*alph*(Bj+Bjp2) +4*alph^2*(DBj+DBjp2) +alph^3*(D2Bj+D2Bjp2)); + Df13 = 1/8*((-6*j-4*j^2)*((Bj+Bjp2)+alph*(DBj+DBjp2))+4*(2*alph*(DBj+DBjp2)+alph^2*(D2Bj+D2Bjp2))+(3*alph^2*(D2Bj+D2Bjp2)+alph^3*(D3Bj+D3Bjp2))); + + f14 = alph*Bjp1; + Df14 = Bjp1+alph*DBjp1; + + f15 = 0.25*((2-4*j^2)*alph*Bjp1 +4*alph^2*DBjp1 +alph^3*D2Bjp1); + Df15 = 0.25*((2-4*j^2)*(Bjp1+alph*DBjp1)+4*(2*alph*DBjp1+alph^2*D2Bjp1)+(3*alph^2*D2Bjp1+alph^3*D3Bjp1)); + + f16 = 0.5*(-alph)*Bjp1 +3*(-alph^2)*Cj +1.5*(-alph^2)*Cjp2; + Df16 = 0.5*((-1)*Bjp1+(-alph)*DBjp1)+3*((-2*alph)*Cj+(-alph^2)*DCj)+1.5*((-2*alph)*Cjp2+(-alph^2)*DCjp2); + + f17 = 1/64*((12+64*j+109*j^2+72*j^3+16*j^4)*Ajp2 +(-12-28*j-16*j^2)*alph*DAjp2 +(6-14*j-8*j^2)*alph^2*D2Ajp2 +8*alph^3*D3Ajp2 +alph^4*D4Ajp2); + Df17 = 1/64*((12+64*j+109*j^2+72*j^3+16*j^4)*DAjp2+(-12-28*j-16*j^2)*(DAjp2+alph*D2Ajp2)+(6-14*j-8*j^2)*(2*alph*D2Ajp2+alph^2*D3Ajp2)+8*(3*alph^2*D3Ajp2+alph^3*D4Ajp2)+(4*alph^3*D4Ajp2+alph^4*D5Ajp2)); + + f18 = 1/16*((12-15*j+4*j^2)*alph*Bjm1 +(8-4*j)*alph^2*DBjm1 +alph^3*D2Bjm1); + Df18 = 1/16*((12-15*j+4*j^2)*(Bjm1+alph*DBjm1)+(8-4*j)*(2*alph*DBjm1+alph^2*D2Bjm1)+(3*alph^2*D2Bjm1+alph^3*D3Bjm1)); + + f19 = 1/8*((6*j-4*j^2)*alph*Bj +(-4+4*j)*alph^2*DBj -alph^3*D2Bj); + Df19 = 1/8*((6*j-4*j^2)*(Bj+alph*DBj)+(-4+4*j)*(2*alph*DBj+alph^2*D2Bj)-(3*alph^2*D2Bj+alph^3*D3Bj)); + + f20 = 1/16*((3*j+4*j^2)*alph*Bjp1 -4*j*alph^2*DBjp1 +alph^3*D2Bjp1); + Df20 = 1/16*((3*j+4*j^2)*(Bjp1+alph*DBjp1)-4*j*(2*alph*DBjp1+alph^2*D2Bjp1)+(3*alph^2*D2Bjp1+alph^3*D3Bjp1)); + + f21 = 1/8*((-12+15*j-4*j^2)*alph*Bjm1 +(-8+4*j)*alph^2*DBjm1 -alph^3*D2Bjm1); + Df21 = 1/8*((-12+15*j-4*j^2)*(Bjm1+alph*DBjm1)+(-8+4*j)*(2*alph*DBjm1+alph^2*D2Bjm1)-(3*alph^2*D2Bjm1+alph^3*D3Bjm1)); + + f22 = 0.25*((6*j+4*j^2)*alph*Bj -4*alph^2*DBj -alph^3*D2Bj); + Df22 = 0.25*((6*j+4*j^2)*(Bj+alph*DBj)-4*(2*alph*DBj+alph^2*D2Bj)-(3*alph^2*D2Bj+alph^3*D3Bj)); + + f23 = 0.25*((6*j+4*j^2)*alph*Bjp2 -4*alph^2*DBjp2 -alph^3*D2Bjp2); + Df23 = 0.25*((6*j+4*j^2)*(Bjp2+alph*DBjp2)-4*(2*alph*DBjp2+alph^2*D2Bjp2)-(3*alph^2*D2Bjp2+alph^3*D3Bjp2)); + + f24 = 0.25*((-6*j+4*j^2)*alph*Bj +(4-4*j)*alph^2*DBj +alph^3*D2Bj); + Df24 = 0.25*((-6*j+4*j^2)*(Bj+alph*DBj)+(4-4*j)*(2*alph*DBj+alph^2*D2Bj)+(3*alph^2*D2Bj+alph^3*D3Bj)); + + f25 = 1/8*((-3*j-4*j^2)*alph*Bjp1 +4*j*alph^2*DBjp1 -alph^3*D2Bjp1); + Df25 = 1/8*((-3*j-4*j^2)*(Bjp1+alph*DBjp1)+4*j*(2*alph*DBjp1+alph^2*D2Bjp1)-(3*alph^2*D2Bjp1+alph^3*D3Bjp1)); + + f26 = 0.5*alph*Bjp1 +0.75*alph^2*Cj +1.5*alph^2*Cjp2; + Df26 = 0.5*(Bjp1+alph*DBjp1)+0.75*(2*alph*Cj+alph^2*DCj)+1.5*(2*alph*Cjp2+alph^2*DCjp2); + + + %1st order semi-major axis functions + f27 = 0.5*(-2*j*Aj-alph*DAj); + Df27 = 0.5*(-2*j*DAj-DAj-alph*D2Aj); + + f28 = 1/16*((2*j-10*j^2+8*j^3)*Aj +(3-7*j+4*j^2)*alph*DAj +(-2-2*j)*alph^2*D2Aj -alph^3*D3Aj); + Df28 = 1/16*((2*j-10*j^2+8*j^3)*DAj+(3-7*j+4*j^2)*(DAj+alph*D2Aj)+(-2-2*j)*(2*alph*D2Aj+alph^2*D3Aj)-(3*alph^2*D3Aj+alph^3*D4Aj)); + + f29 = 1/8*((8*j^3)*Aj +(-2-4*j+4*j^2)*alph*DAj +(-4-2*j)*alph^2*D2Aj -alph^3*D3Aj); + Df29 = 1/8*((8*j^3)*DAj+(-2-4*j+4*j^2)*(DAj+alph*D2Aj)+(-4-2*j)*(2*alph*D2Aj+alph^2*D3Aj)-(3*alph^2*D3Aj+alph^3*D4Aj)); + + f30 = 0.25*((1+2*j)*alph*(Bjm1+Bjp1) +alph^2*(DBjm1+DBjp1)); + Df30 = 0.25*((1+2*j)*((Bjm1+Bjp1)+alph*(DBjm1+DBjp1))+(2*alph*(DBjm1+DBjp1)+alph^2*(D2Bjm1+D2Bjp1))); + + f31 = 0.5*((-1+2*j)*Ajm1+alph*DAjm1); + Df31 = 0.5*((-1+2*j)*DAjm1+DAjm1+alph*D2Ajm1); + + f32 = 1/8*((4-16*j+20*j^2-8*j^3)*Ajm1 +(-4+12*j-4*j^2)*alph*DAjm1 +(3+2*j)*alph^2*D2Ajm1 +alph^3*D3Ajm1); + Df32 = 1/8*((4-16*j+20*j^2-8*j^3)*DAjm1+(-4+12*j-4*j^2)*(DAjm1+alph*D2Ajm1)+(3+2*j)*(2*alph*D2Ajm1+alph^2*D3Ajm1)+(3*alph^2*D3Ajm1+alph^3*D4Ajm1)); + + f33 = 1/16*((-2-j+10*j^2-8*j^3)*Ajm1 +(2+9*j-4*j^2)*alph*DAjm1 +(5+2*j)*alph^2*D2Ajm1 +alph^3*D3Ajm1); + Df33 = 1/16*((-2-j+10*j^2-8*j^3)*DAjm1+(2+9*j-4*j^2)*(DAjm1+alph*D2Ajm1)+(5+2*j)*(2*alph*D2Ajm1+alph^2*D3Ajm1)+(3*alph^2*D3Ajm1+alph^3*D4Ajm1)); + + f34 = 0.25*(-2*j*alph*(Bjm2+Bj) -alph^2*(DBjm2+DBj)); + Df34 = 0.25*(-2*j*((Bjm2+Bj)+alph*(DBjm2+DBj))-(2*alph*(DBjm2+DBj)+alph^2*(D2Bjm2+D2Bj))); + + f35 = 1/16*((1-j-10*j^2-8*j^3)*Ajp1 +(-1-j-4*j^2)*alph*DAjp1 +(3+2*j)*alph^2*D2Ajp1 +alph^3*D3Ajp1); + Df35 = 1/16*((1-j-10*j^2-8*j^3)*DAjp1+(-1-j-4*j^2)*(DAjp1+alph*D2Ajp1)+(3+2*j)*(2*alph*D2Ajp1+alph^2*D3Ajp1)+(3*alph^2*D3Ajp1+alph^3*D4Ajp1)); + + f36 = 1/16*((-8+32*j-30*j^2+8*j^3)*Ajm2 +(8-17*j+4*j^2)*alph*DAjm2 +(-4-2*j)*alph^2*D2Ajm2 -alph^3*D3Ajm2); + Df36 = 1/16*((-8+32*j-30*j^2+8*j^3)*DAjm2+(8-17*j+4*j^2)*(DAjm2+alph*D2Ajm2)+(-4-2*j)*(2*alph*D2Ajm2+alph^2*D3Ajm2)-(3*alph^2*D3Ajm2+alph^3*D4Ajm2)); + + f37 = 0.25*((-5+2*j)*alph*Bjm1 -alph^2*DBjm1); + Df37 = 0.25*((-5+2*j)*(Bjm1+alph*DBjm1)-(2*alph*DBjm1+alph^2*D2Bjm1)); + + f38 = 0.25*(-2*j*alph*Bj +alph^2*DBj); + Df38 = 0.25*(-2*j*(Bj+alph*DBj)+(2*alph*DBj+alph^2*D2Bj)); + + f39 = 0.5*((-1-2*j)*alph*Bjm1 -alph^2*DBjm1); + Df39 = 0.5*((-1-2*j)*(Bjm1+alph*DBjm1)-(2*alph*DBjm1+alph^2*D2Bjm1)); + + f40 = 0.5*((-1-2*j)*alph*Bjp1 -alph^2*DBjp1); + Df40 = 0.5*((-1-2*j)*(Bjp1+alph*DBjp1)-(2*alph*DBjp1+alph^2*D2Bjp1)); + + f41 = 0.5*((5-2*j)*alph*Bjm1 +alph^2*DBjm1); + Df41 = 0.5*((5-2*j)*(Bjm1+alph*DBjm1)+(2*alph*DBjm1+alph^2*D2Bjm1)); + + f42 = 0.5*(2*j*alph*Bjm2 +alph^2*DBjm2); + Df42 = 0.5*(2*j*(Bjm2+alph*DBjm2)+(2*alph*DBjm2+alph^2*D2Bjm2)); + + f43 = 0.5*(2*j*alph*Bj +alph^2*DBj); + Df43 = 0.5*(2*j*(Bj+alph*DBj)+(2*alph*DBj+alph^2*D2Bj)); + + f44 = 0.5*(2*j*alph*Bj -alph^2*DBj); + Df44 = 0.5*(2*j*(Bj+alph*DBj)-(2*alph*DBj+alph^2*D2Bj)); + + %2nd order semi-major axis functions + f45 = 1/8*((-5*j+4*j^2)*Aj +(-2+4*j)*alph*DAj +alph^2*D2Aj); + Df45 = 1/8*((-5*j+4*j^2)*DAj+(-2+4*j)*(DAj+alph*D2Aj)+(2*alph*D2Aj+alph^2*D3Aj)); + + f46 = 1/96*((22*j-64*j^2+60*j^3-16*j^4)*Aj +(16-46*j+48*j^2-16*j^3)*alph*DAj +(-12+9*j)*alph^2*D2Aj +4*j*alph^3*D3Aj +alph^4*D4Aj); + Df46 = 1/96*((22*j-64*j^2+60*j^3-16*j^4)*DAj+(16-46*j+48*j^2-16*j^3)*(DAj+alph*D2Aj)+(-12+9*j)*(2*alph*D2Aj+alph^2*D3Aj)+4*j*(3*alph^2*D3Aj+alph^3*D4Aj)+(4*alph^3*D4Aj+alph^4*D5Aj)); + + f47 = 1/32*((20*j^3-16*j^4)*Aj +(-4-2*j+16*j^2-16*j^3)*alph*DAj +(-2+11*j)*alph^2*D2Aj +(4+4*j)*alph^3*D3Aj +alph^4*D4Aj); + Df47 = 1/32*((20*j^3-16*j^4)*DAj+(-4-2*j+16*j^2-16*j^3)*(DAj+alph*D2Aj)+(-2+11*j)*(2*alph*D2Aj+alph^2*D3Aj)+(4+4*j)*(3*alph^2*D3Aj+alph^3*D4Aj)+(4*alph^3*D4Aj+alph^4*D5Aj)); + + f48 = 1/16*((2+j-4*j^2)*alph*(Bjm1+Bjp1) -4*j*alph^2*(DBjm1+DBjp1) -alph^3*(D2Bjm1+D2Bjp1)); + Df48 = 1/16*((2+j-4*j^2)*((Bjm1+Bjp1)+alph*(DBjm1+DBjp1))-4*j*(2*alph*(DBjm1+DBjp1)+alph^2*(D2Bjm1+D2Bjp1))-(3*alph^2*(D2Bjm1+D2Bjp1)+alph^3*(D3Bjm1+D3Bjp1))); + + f49 = 0.25*((-2+6*j-4*j^2)*Ajm1 +(2-4*j)*alph*DAjm1 -alph^2*D2Ajm1); + Df49 = 0.25*((-2+6*j-4*j^2)*DAjm1+(2-4*j)*(DAjm1+alph*D2Ajm1)-(2*alph*D2Ajm1+alph^2*D3Ajm1)); + + f50 = 1/32*((20-86*j+126*j^2-76*j^3+16*j^4)*Ajm1 +(-20+74*j-64*j^2+16*j^3)*alph*DAjm1 +(14-17*j)*alph^2*D2Ajm1 +(-2-4*j)*alph^3*D3Ajm1 -alph^4*D4Ajm1); %in this term there appears Aj- instead of Aj-1; this is explained in the erratum http://ssdbook.maths.qmul.ac.uk/errors.pdf + Df50 = 1/32*((20-86*j+126*j^2-76*j^3+16*j^4)*DAjm1+(-20+74*j-64*j^2+16*j^3)*(DAjm1+alph*D2Ajm1)+(14-17*j)*(2*alph*D2Ajm1+alph^2*D3Ajm1)+(-2-4*j)*(3*alph^2*D3Ajm1+alph^3*D4Ajm1)-(4*alph^3*D4Ajm1+alph^4*D5Ajm1)); + + f51 = 1/32*((-4+2*j+22*j^2-36*j^3+16*j^4)*Ajm1 +(4+6*j-32*j^2+16*j^3)*alph*DAjm1 +(-2-19*j)*alph^2*D2Ajm1 +(-6-4*j)*alph^3*D3Ajm1 -alph^4*D4Ajm1); + Df51 = 1/32*((-4+2*j+22*j^2-36*j^3+16*j^4)*DAjm1+(4+6*j-32*j^2+16*j^3)*(DAjm1+alph*D2Ajm1)+(-2-19*j)*(2*alph*D2Ajm1+alph^2*D3Ajm1)+(-6-4*j)*(3*alph^2*D3Ajm1+alph^3*D4Ajm1)-(4*alph^3*D4Ajm1+alph^4*D5Ajm1)); + + f52 = 1/8*((-2*j+4*j^2)*alph*(Bjm2+Bj) +4*j*alph^2*(DBjm2+DBj) +alph^3*(D2Bjm2+D2Bj)); + Df52 = 1/8*((-2*j+4*j^2)*((Bjm2+Bj)+alph*(DBjm2+DBj))+4*j*(2*alph*(DBjm2+DBj)+alph^2*(D2Bjm2+D2Bj))+(3*alph^2*(D2Bjm2+D2Bj)+alph^3*(D3Bjm2+D3Bj))); + + f53 = 1/8*((2-7*j+4*j^2)*Ajm2 +(-2+4*j)*alph*DAjm2 +alph^2*D2Ajm2); + Df53 = 1/8*((2-7*j+4*j^2)*DAjm2+(-2+4*j)*(DAjm2+alph*D2Ajm2)+(2*alph*D2Ajm2+alph^2*D3Ajm2)); + + f54 = 1/32*((-32+144*j-184*j^2+92*j^3-16*j^4)*Ajm2 +(32-102*j+80*j^2-16*j^3)*alph*DAjm2 +(-16+25*j)*alph^2*D2Ajm2 +(4+4*j)*alph^3*D3Ajm2 +alph^4*D4Ajm2); + Df54 = 1/32*((-32+144*j-184*j^2+92*j^3-16*j^4)*DAjm2+(32-102*j+80*j^2-16*j^3)*(DAjm2+alph*D2Ajm2)+(-16+25*j)*(2*alph*D2Ajm2+alph^2*D3Ajm2)+(4+4*j)*(3*alph^2*D3Ajm2+alph^3*D4Ajm2)+(4*alph^3*D4Ajm2+alph^4*D5Ajm2)); + + f55 = 1/96*((12-14*j-40*j^2+52*j^3-16*j^4)*Ajm2 +(-12-10*j+48*j^2-16*j^3)*alph*DAjm2 +(6+27*j)*alph^2*D2Ajm2 +(8+4*j)*alph^3*D3Ajm2 +alph^4*D4Ajm2); + Df55 = 1/96*((12-14*j-40*j^2+52*j^3-16*j^4)*DAjm2+(-12-10*j+48*j^2-16*j^3)*(DAjm2+alph*D2Ajm2)+(6+27*j)*(2*alph*D2Ajm2+alph^2*D3Ajm2)+(8+4*j)*(3*alph^2*D3Ajm2+alph^3*D4Ajm2)+(4*alph^3*D4Ajm2+alph^4*D5Ajm2)); + + f56 = 1/16*((3*j-4*j^2)*alph*(Bjm3+Bjm1) -4*j*alph^2*(DBjm3+DBjm1) -alph^3*(D2Bjm3+D2Bjm1)); + Df56 = 1/16*((3*j-4*j^2)*((Bjm3+Bjm1)+alph*(DBjm3+DBjm1))-4*j*(2*alph*(DBjm3+DBjm1)+alph^2*(D2Bjm3+D2Bjm1))-(3*alph^2*(D2Bjm3+D2Bjm1)+alph^3*(D3Bjm3+D3Bjm1))); + + f57 = 0.5*alph*Bjm1; + Df57 = 0.5*(Bjm1+alph*DBjm1); + + f58 = 1/8*((-14+16*j-4*j^2)*alph*Bjm1 +4*alph^2*DBjm1 +alph^3*D2Bjm1); + Df58 = 1/8*((-14+16*j-4*j^2)*(Bjm1+alph*DBjm1)+4*(2*alph*DBjm1+alph^2*D2Bjm1)+(3*alph^2*D2Bjm1+alph^3*D3Bjm1)); + + f59 = 1/8*((2-4*j^2)*alph*Bjm1 +4*alph^2*DBjm1 +alph^3*D2Bjm1); + Df59 = 1/8*((2-4*j^2)*(Bjm1+alph*DBjm1)+4*(2*alph*DBjm1+alph^2*D2Bjm1)+(3*alph^2*D2Bjm1+alph^3*D3Bjm1)); + + f60 = 0.75*(-alph^2)*(Cjm2+Cj); + Df60 = 0.75*((-2*alph)*(Cjm2+Cj)+(-alph^2)*(DCjm2+DCj)); + + f61 = 0.5*(-alph)*Bjm1 +0.75*(-alph^2)*Cjm2 +15/4*(-alph^2)*Cj; + Df61 = 0.5*((-1)*Bjm1+(-alph)*DBjm1)+0.75*((-2*alph)*Cjm2+(-alph^2)*DCjm2)+15/4*((-2*alph)*Cj+(-alph^2)*DCj); + + f62 = -alph*Bjm1; + Df62 = -Bjm1-alph*DBjm1; + + f63 = 0.25*((14-16*j+4*j^2)*alph*Bjm1 -4*alph^2*DBjm1 -alph^3*D2Bjm1); + Df63 = 0.25*((14-16*j+4*j^2)*(Bjm1+alph*DBjm1)-4*(2*alph*DBjm1+alph^2*D2Bjm1)-(3*alph^2*D2Bjm1+alph^3*D3Bjm1)); + + f64 = 0.25*((-2+4*j^2)*alph*Bjm1 -4*alph^2*DBjm1 -alph^3*D2Bjm1); + Df64 = 0.25*((-2+4*j^2)*(Bjm1+alph*DBjm1)-4*(2*alph*DBjm1+alph^2*D2Bjm1)-(3*alph^2*D2Bjm1+alph^3*D3Bjm1)); + + f65 = 0.5*alph*Bjm1 +3*alph^2*Cjm2 +1.5*alph^2*Cj; + Df65 = 0.5*(Bjm1+alph*DBjm1)+3*(2*alph*Cjm2+alph^2*DCjm2)+1.5*(2*alph*Cj+alph^2*DCj); + + f66 = 0.5*alph*Bjm1 +1.5*alph^2*Cjm2 +3*alph^2*Cj; + Df66 = 0.5*(Bjm1+alph*DBjm1)+1.5*(2*alph*Cjm2+alph^2*DCjm2)+3*(2*alph*Cj+alph^2*DCj); + + f67 = 0.5*(-alph)*Bjm1 +15/4*(-alph^2)*Cjm2 +0.75*(-alph^2)*Cj; + Df67 = 0.5*((-1)*Bjm1+(-alph)*DBjm1)+15/4*((-2*alph)*Cjm2+(-alph^2)*DCjm2)+0.75*((-2*alph)*Cj+(-alph^2)*DCj); + + f68 = 1/96*((4-2*j-26*j^2-4*j^3+16*j^4)*Ajp1 +(-4-2*j+16*j^3)*alph*DAjp1 +(6-3*j)*alph^2*D2Ajp1 +(-2-4*j)*alph^3*D3Ajp1 -alph^4*D4Ajp1); + Df68 = 1/96*((4-2*j-26*j^2-4*j^3+16*j^4)*DAjp1+(-4-2*j+16*j^3)*(DAjp1+alph*D2Ajp1)+(6-3*j)*(2*alph*D2Ajp1+alph^2*D3Ajp1)+(-2-4*j)*(3*alph^2*D3Ajp1+alph^3*D4Ajp1)-(4*alph^3*D4Ajp1+alph^4*D5Ajp1)); + + f69 = 1/96*((36-186*j+238*j^2-108*j^3+16*j^4)*Ajm3 +(-36+130*j-96*j^2+16*j^3)*alph*DAjm3 +(18-33*j)*alph^2*D2Ajm3 +(-6-4*j)*alph^3*D3Ajm3 -alph^4*D4Ajm3); + Df69 = 1/96*((36-186*j+238*j^2-108*j^3+16*j^4)*DAjm3+(-36+130*j-96*j^2+16*j^3)*(DAjm3+alph*D2Ajm3)+(18-33*j)*(2*alph*D2Ajm3+alph^2*D3Ajm3)+(-6-4*j)*(3*alph^2*D3Ajm3+alph^3*D4Ajm3)-(4*alph^3*D4Ajm3+alph^4*D5Ajm3)); + + f70 = 1/8*((-14*j+4*j^2)*alph*Bjm2 -8*alph^2*DBjm2 -alph^3*D2Bjm2); + Df70 = 1/8*((-14*j+4*j^2)*(Bjm2+alph*DBjm2)-8*(2*alph*DBjm2+alph^2*D2Bjm2)-(3*alph^2*D2Bjm2+alph^3*D3Bjm2)); + + f71 = 1/8*((-2*j+4*j^2)*alph*Bj -alph^3*D2Bj); + Df71 = 1/8*((-2*j+4*j^2)*(Bj+alph*DBj)-(3*alph^2*D2Bj+alph^3*D3Bj)); + + f72 = 1/8*((-2-j+4*j^2)*alph*Bjm1 +4*j*alph^2*DBjm1 +alph^3*D2Bjm1); + Df72 = 1/8*((-2-j+4*j^2)*(Bjm1+alph*DBjm1)+4*j*(2*alph*DBjm1+alph^2*D2Bjm1)+(3*alph^2*D2Bjm1+alph^3*D3Bjm1)); + + f73 = 1/8*((-2-j+4*j^2)*alph*Bjp1 +4*j*alph^2*DBjp1 +alph^3*D2Bjp1); + Df73 = 1/8*((-2-j+4*j^2)*(Bjp1+alph*DBjp1)+4*j*(2*alph*DBjp1+alph^2*D2Bjp1)+(3*alph^2*D2Bjp1+alph^3*D3Bjp1)); + + f74 = 0.25*((2*j-4*j^2)*alph*Bjm2 -4*j*alph^2*DBjm2 -alph^3*D2Bjm2); + Df74 = 0.25*((2*j-4*j^2)*(Bjm2+alph*DBjm2)-4*j*(2*alph*DBjm2+alph^2*D2Bjm2)-(3*alph^2*D2Bjm2+alph^3*D3Bjm2)); + + f75 = 0.25*((2*j-4*j^2)*alph*Bj -4*j*alph^2*DBj -alph^3*D2Bj); + Df75 = 0.25*((2*j-4*j^2)*(Bj+alph*DBj)-4*j*(2*alph*DBj+alph^2*D2Bj)-(3*alph^2*D2Bj+alph^3*D3Bj)); + + f76 = 0.25*((14*j-4*j^2)*alph*Bjm2 +8*alph^2*DBjm2 +alph^3*D2Bjm2); + Df76 = 0.25*((14*j-4*j^2)*(Bjm2+alph*DBjm2)+8*(2*alph*DBjm2+alph^2*D2Bjm2)+(3*alph^2*D2Bjm2+alph^3*D3Bjm2)); + + f77 = 0.25*((2*j-4*j^2)*alph*Bj +alph^3*D2Bj); + Df77 = 0.25*((2*j-4*j^2)*(Bj+alph*DBj)+(3*alph^2*D2Bj+alph^3*D3Bj)); + + f78 = 1/8*((-3*j+4*j^2)*alph*Bjm3 +4*j*alph^2*DBjm3 +alph^3*D2Bjm3); + Df78 = 1/8*((-3*j+4*j^2)*(Bjm3+alph*DBjm3)+4*j*(2*alph*DBjm3+alph^2*D2Bjm3)+(3*alph^2*D2Bjm3+alph^3*D3Bjm3)); + + f79 = 1/8*((-3*j+4*j^2)*alph*Bjm1 +4*j*alph^2*DBjm1 +alph^3*D2Bjm1); + Df79 = 1/8*((-3*j+4*j^2)*(Bjm1+alph*DBjm1)+4*j*(2*alph*DBjm1+alph^2*D2Bjm1)+(3*alph^2*D2Bjm1+alph^3*D3Bjm1)); + + f80 = 1.5*alph^2*Cj; + Df80 = 1.5*(2*alph*Cj+alph^2*DCj); + + f81 = 1.5*alph^2*Cjm2; + Df81 = 1.5*(2*alph*Cjm2+alph^2*DCjm2); + + %3rd order semi-major axis functions + f82 = 1/48*((-26*j+30*j^2-8*j^3)*Aj +(-9+27*j-12*j^2)*alph*DAj +(6-6*j)*alph^2*D2Aj -alph^3*D3Aj); + Df82 = 1/48*((-26*j+30*j^2-8*j^3)*DAj+(-9+27*j-12*j^2)*(DAj+alph*D2Aj)+(6-6*j)*(2*alph*D2Aj+alph^2*D3Aj)-(3*alph^2*D3Aj+alph^3*D4Aj)); + + f83 = 1/16*((-9+31*j-30*j^2+8*j^3)*Ajm1 +(9-25*j+12*j^2)*alph*DAjm1 +(-5+6*j)*alph^2*D2Ajm1 +alph^3*D3Ajm1); + Df83 = 1/16*((-9+31*j-30*j^2+8*j^3)*DAjm1+(9-25*j+12*j^2)*(DAjm1+alph*D2Ajm1)+(-5+6*j)*(2*alph*D2Ajm1+alph^2*D3Ajm1)+(3*alph^2*D3Ajm1+alph^3*D4Ajm1)); + + f84 = 1/16*((8-32*j+30*j^2-8*j^3)*Ajm2 +(-8+23*j-12*j^2)*alph*DAjm2 +(4-6*j)*alph^2*D2Ajm2 -alph^3*D3Ajm2); + Df84 = 1/16*((8-32*j+30*j^2-8*j^3)*DAjm2+(-8+23*j-12*j^2)*(DAjm2+alph*D2Ajm2)+(4-6*j)*(2*alph*D2Ajm2+alph^2*D3Ajm2)-(3*alph^2*D3Ajm2+alph^3*D4Ajm2)); + + f85 = 1/48*((-6+29*j-30*j^2+8*j^3)*Ajm3 +(6-21*j+12*j^2)*alph*DAjm3 +(-3+6*j)*alph^2*D2Ajm3 +alph^3*D3Ajm3); + Df85 = 1/48*((-6+29*j-30*j^2+8*j^3)*DAjm3+(6-21*j+12*j^2)*(DAjm3+alph*D2Ajm3)+(-3+6*j)*(2*alph*D2Ajm3+alph^2*D3Ajm3)+(3*alph^2*D3Ajm3+alph^3*D4Ajm3)); + + f86 = 0.25*((3-2*j)*alph*Bjm1 -alph^2*DBjm1); + Df86 = 0.25*((3-2*j)*(Bjm1+alph*DBjm1)-(2*alph*DBjm1+alph^2*D2Bjm1)); + + f87 = 0.25*(2*j*alph*Bjm2 +alph^2*DBjm2); + Df87 = 0.25*(2*j*(Bjm2+alph*DBjm2)+(2*alph*DBjm2+alph^2*D2Bjm2)); + + f88 = 0.5*((-3+2*j)*alph*Bjm1 +alph^2*DBjm1); + Df88 = 0.5*((-3+2*j)*(Bjm1+alph*DBjm1)+(2*alph*DBjm1+alph^2*D2Bjm1)); + + f89 = 0.5*(-2*j*alph*Bjm2 -alph^2*DBjm2); + Df89 = 0.5*(-2*j*(Bjm2+alph*DBjm2)-(2*alph*DBjm2+alph^2*D2Bjm2)); + + %4th order semi-major axis functions + f90 = 1/384*((-206*j+283*j^2-120*j^3+16*j^4)*Aj +(-64+236*j-168*j^2+32*j^3)*alph*DAj +(48-78*j+24*j^2)*alph^2*D2Aj +(-12+8*j)*alph^3*D3Aj +alph^4*D4Aj); + Df90 = 1/384*((-206*j+283*j^2-120*j^3+16*j^4)*DAj+(-64+236*j-168*j^2+32*j^3)*(DAj+alph*D2Aj)+(48-78*j+24*j^2)*(2*alph*D2Aj+alph^2*D3Aj)+(-12+8*j)*(3*alph^2*D3Aj+alph^3*D4Aj)+(4*alph^3*D4Aj+alph^4*D5Aj)); + + f91 = 1/96*((-64+238*j-274*j^2+116*j^3-16*j^4)*Ajm1 +(64-206*j+156*j^2-32*j^3)*alph*DAjm1 +(-36+69*j-24*j^2)*alph^2*D2Ajm1 +(10-8*j)*alph^3*D3Ajm1 -alph^4*D4Ajm1); + Df91 = 1/96*((-64+238*j-274*j^2+116*j^3-16*j^4)*DAjm1+(64-206*j+156*j^2-32*j^3)*(DAjm1+alph*D2Ajm1)+(-36+69*j-24*j^2)*(2*alph*D2Ajm1+alph^2*D3Ajm1)+(10-8*j)*(3*alph^2*D3Ajm1+alph^3*D4Ajm1)-(4*alph^3*D4Ajm1+alph^4*D5Ajm1)); + + f92 = 1/64*((52-224*j+259*j^2-112*j^3+16*j^4)*Ajm2 +(-52+176*j-144*j^2+32*j^3)*alph*DAjm2 +(26-60*j+24*j^2)*alph^2*D2Ajm2 +(-8+8*j)*alph^3*D3Ajm2 +alph^4*D4Ajm2); + Df92 = 1/64*((52-224*j+259*j^2-112*j^3+16*j^4)*DAjm2+(-52+176*j-144*j^2+32*j^3)*(DAjm2+alph*D2Ajm2)+(26-60*j+24*j^2)*(2*alph*D2Ajm2+alph^2*D3Ajm2)+(-8+8*j)*(3*alph^2*D3Ajm2+alph^3*D4Ajm2)+(4*alph^3*D4Ajm2+alph^4*D5Ajm2)); + + f93 = 1/96*((-36+186*j-238*j^2+108*j^3-16*j^4)*Ajm3 +(36-146*j+132*j^2-32*j^3)*alph*DAjm3 +(-18+51*j-24*j^2)*alph^2*D2Ajm3 +(6-8*j)*alph^3*D3Ajm3 -alph^4*D4Ajm3); + Df93 = 1/96*((-36+186*j-238*j^2+108*j^3-16*j^4)*DAjm3+(36-146*j+132*j^2-32*j^3)*(DAjm3+alph*D2Ajm3)+(-18+51*j-24*j^2)*(2*alph*D2Ajm3+alph^2*D3Ajm3)+(6-8*j)*(3*alph^2*D3Ajm3+alph^3*D4Ajm3)-(4*alph^3*D4Ajm3+alph^4*D5Ajm3)); + + f94 = 1/384*((24-146*j+211*j^2-104*j^3+16*j^4)*Ajm4 +(-24+116*j-120*j^2+32*j^3)*alph*DAjm4 +(12-42*j+24*j^2)*alph^2*D2Ajm4 +(-4+8*j)*alph^3*D3Ajm4 +alph^4*D4Ajm4); + Df94 = 1/384*((24-146*j+211*j^2-104*j^3+16*j^4)*DAjm4+(-24+116*j-120*j^2+32*j^3)*(DAjm4+alph*D2Ajm4)+(12-42*j+24*j^2)*(2*alph*D2Ajm4+alph^2*D3Ajm4)+(-4+8*j)*(3*alph^2*D3Ajm4+alph^3*D4Ajm4)+(4*alph^3*D4Ajm4+alph^4*D5Ajm4)); + + f95 = 1/16*((16-17*j+4*j^2)*alph*Bjm1 +(-8+4*j)*alph^2*DBjm1 +alph^3*D2Bjm1); + Df95 = 1/16*((16-17*j+4*j^2)*(Bjm1+alph*DBjm1)+(-8+4*j)*(2*alph*DBjm1+alph^2*D2Bjm1)+(3*alph^2*D2Bjm1+alph^3*D3Bjm1)); + + f96 = 1/8*((10*j-4*j^2)*alph*Bjm2 +(4-4*j)*alph^2*DBjm2 -alph^3*D2Bjm2); + Df96 = 1/8*((10*j-4*j^2)*(Bjm2+alph*DBjm2)+(4-4*j)*(2*alph*DBjm2+alph^2*D2Bjm2)-(3*alph^2*D2Bjm2+alph^3*D3Bjm2)); + + f97 = 1/16*((-3*j+4*j^2)*alph*Bjm3 +4*j*alph^2*DBjm3 +alph^3*D2Bjm3); + Df97 = 1/16*((-3*j+4*j^2)*(Bjm3+alph*DBjm3)+4*j*(2*alph*DBjm3+alph^2*D2Bjm3)+(3*alph^2*D2Bjm3+alph^3*D3Bjm3)); + + f98 = 3/8*alph^2*Cjm2; + Df98 = 3/8*(2*alph*Cjm2+alph^2*DCjm2); + + f99 = 1/8*((-16+17*j-4*j^2)*alph*Bjm1 +(8-4*j)*alph^2*DBjm1 -alph^3*D2Bjm1); + Df99 = 1/8*((-16+17*j-4*j^2)*(Bjm1+alph*DBjm1)+(8-4*j)*(2*alph*DBjm1+alph^2*D2Bjm1)-(3*alph^2*D2Bjm1+alph^3*D3Bjm1)); + + f100 = 0.25*((-10*j+4*j^2)*alph*Bjm2 +(-4+4*j)*alph^2*DBjm2 +alph^3*D2Bjm2); + Df100 = 0.25*((-10*j+4*j^2)*(Bjm2+alph*DBjm2)+(-4+4*j)*(2*alph*DBjm2+alph^2*D2Bjm2)+(3*alph^2*D2Bjm2+alph^3*D3Bjm2)); + + f101 = 1/8*((3*j-4*j^2)*alph*Bjm3 -4*j*alph^2*DBjm3 -alph^3*D2Bjm3); + Df101 = 1/8*((3*j-4*j^2)*(Bjm3+alph*DBjm3)-4*j*(2*alph*DBjm3+alph^2*D2Bjm3)-(3*alph^2*D2Bjm3+alph^3*D3Bjm3)); + + f102 = 1.5*(-alph^2)*Cjm2; + Df102 = 1.5*((-2*alph)*Cjm2+(-alph^2)*DCjm2); + + f103 = 9/4*(alph^2*Cjm2); + Df103 = 9/4*(2*alph*Cjm2+alph^2*DCjm2); + + + end + + + %% PrepareExpansionTerms; + + %Write down the disturbing function terms powers and resonant orders - Solar + %System Dynamics (Murray & Dermott 1999), Appendix B + + for DummyIdx = 1:1 + + if ~isfield(Terms,'k') %constant fields, not depending on the value of alpha or j, should be written once and for all + % if isempty(Terms) %constant fields, not depending on the value of alpha or j, should be written once and for all + + Terms.k = [zeros(1,38) ... %0th order + ones(1,22) ... %1st order + 2*ones(1,46) ... %2nd order + 3*ones(1,10) ... %3rd order + 4*ones(1,19)]; %4th order + + Terms.A = [0 2 0 0 0 4 2 0 2 0 2 0 0 0 0 ... + 1 3 1 1 1 0 2 0 0 0 2 2 1 0 2 1 1 1 0 2 1 0 0 ... %0th order + 1 3 1 1 1 0 2 0 0 0 2 1 1 0 1 1 1 0 0 0 1 0 ... %1st order + 2 4 2 2 2 1 3 1 1 1 0 2 0 0 0 0 2 0 0 0 ... + 0 2 0 0 0 0 2 0 0 0 3 1 1 1 2 2 1 1 1 1 0 0 0 1 1 0 ... %2nd order + 3 2 1 0 1 0 1 0 1 0 ... %3rd order + 4 3 2 1 0 2 1 0 0 2 1 0 0 2 1 0 0 0 0 ... %4th order + ]; + + Terms.Ap = [0 0 2 0 0 0 2 4 0 2 0 2 0 0 0 ... + 1 1 3 1 1 0 0 2 0 0 2 0 1 2 0 1 1 1 2 0 1 2 0 ... %0th order + 0 0 2 0 0 1 1 3 1 1 1 2 0 1 0 0 0 1 1 1 0 1 ... %1st order + 0 0 2 0 0 1 1 3 1 1 2 2 4 2 2 0 0 2 0 0 ... + 0 0 2 0 0 0 0 2 0 0 1 3 1 1 0 0 1 1 1 1 2 2 0 1 1 0 ... %2nd order + 0 1 2 3 0 1 0 1 0 1 ... %3rd order + 0 1 2 3 4 0 1 2 0 0 1 2 0 0 1 2 0 0 0 ... %4th order + ]; + + Terms.B = [0 0 0 2 0 0 0 0 2 2 0 0 4 0 2 ... + 0 0 0 2 0 1 1 1 3 1 0 2 2 2 1 1 1 1 1 0 0 0 2 ... %0th order + 0 0 0 2 0 0 0 0 2 0 0 0 2 2 1 1 1 1 1 1 0 0 ... %1st order + 0 0 0 2 0 0 0 0 2 0 0 0 0 2 0 2 2 2 4 2 ... + 1 1 1 3 1 0 0 0 2 0 0 0 2 2 1 1 1 1 1 1 1 1 3 0 0 1 ... %2nd order + 0 0 0 0 2 2 1 1 0 0 ... %3rd order + 0 0 0 0 0 2 2 2 4 1 1 1 3 0 0 0 2 1 0 ... %4th order + ]; + + Terms.Bp = [0 0 0 0 2 0 0 0 0 0 2 2 0 4 2 ... + 0 0 0 0 2 1 1 1 1 3 0 0 0 0 1 1 1 1 1 2 2 2 2 ... %0th order + 0 0 0 0 2 0 0 0 0 2 0 0 0 0 1 1 1 1 1 1 2 2 ... %1st order + 0 0 0 0 2 0 0 0 0 2 0 0 0 0 2 0 0 0 0 2 ... + 1 1 1 1 3 2 2 2 2 4 0 0 0 0 1 1 1 1 1 1 1 1 1 2 2 3 ... %2nd order + 0 0 0 0 0 0 1 1 2 2 ... %3rd order + 0 0 0 0 0 0 0 0 0 1 1 1 1 2 2 2 2 3 4 ... %4th order + ]; + + Terms.C = [zeros(1,15) ... + ones(1,5) zeros(1,5) 2 -2 -1 0 -2 1 1 -1 0 -2 -1 0 0 ... %0th order + 1 1 1 1 1 0 0 0 0 0 2 -1 -1 0 1 1 -1 0 0 0 -1 0 ... %1st order + 2 2 2 2 2 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 ... + 0 0 0 0 0 0 0 0 0 0 3 -1 -1 1 2 2 1 1 -1 1 0 0 0 -1 1 0 ... %2nd order + 3 2 1 0 1 0 1 0 1 0 ... %3rd order + 4 3 2 1 0 2 1 0 0 2 1 0 0 2 1 0 0 0 0 ... %4th order + ]; + + Terms.Cp = [zeros(1,15) ... + -ones(1,5) zeros(1,5) -2 0 -1 -2 0 -1 -1 -1 -2 0 -1 -2 0 ... %0th order + 0 0 0 0 0 1 1 1 1 1 -1 2 0 -1 0 0 0 1 1 -1 0 -1 ... %1st order + 0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 0 0 0 0 0 ... + 0 0 0 0 0 0 0 0 0 0 -1 3 1 -1 0 0 1 1 1 -1 2 2 0 1 -1 0 ... %2nd order + 0 1 2 3 0 1 0 1 0 1 ... %3rd order + 0 1 2 3 4 0 1 2 0 0 1 2 0 0 1 2 0 0 0 .... %4th order + ]; + + Terms.D = [zeros(1,15) ... + zeros(1,5) ones(1,5) 0 2 2 2 1 -1 1 1 1 0 0 0 2 ... %0th order + zeros(1,10) 0 0 2 2 -1 1 1 -1 1 1 0 0 ... %1st order + zeros(1,15) 2 2 2 2 2 1 1 1 1 1 0 0 0 0 0 ... + 0 0 2 2 -1 1 -1 1 1 1 -1 1 3 0 0 -1 ... %2nd order + 0 0 0 0 2 2 1 1 0 0 ... %3rd order + 0 0 0 0 0 2 2 2 4 1 1 1 3 0 0 0 2 1 0 ... %4th order + ]; + + Terms.Dp = [zeros(1,15) ... + zeros(1,5) -ones(1,5) 0 0 0 0 1 1 -1 1 1 2 2 2 -2 ... %0th order + zeros(1,10) 0 0 0 0 1 -1 1 1 -1 1 2 2 ... %1st order + zeros(1,20) 1 1 1 1 1 2 2 2 2 2 ... + 0 0 0 0 1 -1 1 -1 1 1 1 -1 -1 2 2 3 ... %2nd order + 0 0 0 0 0 0 1 1 2 2 ... %3rd order + 0 0 0 0 0 0 0 0 0 1 1 1 1 2 2 2 2 3 4 ... %4th order + ]; + + Terms.N = length(Terms.A); + end + + Terms.f = [f1 f2 f2 f3 f3 f4 f5 f6 f7 f7 f7 f7 f8 f8 f9 ... + f10 f11 f12 f13 f13 f14 f15 f15 f16 f16 f17 f18 f19 f20 f21 f22 f23 f24 f25 f18 f19 f20 f26 ... %0th order + f27 f28 f29 f30 f30 f31 f32 f33 f34 f34 f35 f36 f37 f38 f39 f40 f41 f42 f43 f44 f37 f38 ... %1st order + f45 f46 f47 f48 f48 f49 f50 f51 f52 f52 f53 f54 f55 f56 f56 f57 f58 f59 f60 f61 f62 f63 f64 f65 f66 f57 f58 f59 f67 f60 ... + f68 f69 f70 f71 f72 f73 f74 f75 f76 f77 f78 f79 f80 f70 f71 f81 ... %2nd order + f82 f83 f84 f85 f86 f87 f88 f89 f86 f87 ... %3rd order + f90 f91 f92 f93 f94 f95 f96 f97 f98 f99 f100 f101 f102 f95 f96 f97 f103 f102 f98 ... %4th order + + ]; + + %indirect terms due to an external perturber. The prefactor of alph + %comes from the definition of the disturbing function (Solar System + %Dynamics equation 6.44). + + Terms.fE = alph*[(j==1)*[-1 0.5 0.5 1 1 1/64 -0.25 1/64 -0.5 -0.5 -0.5 -0.5 0 0 -1] ... + (j==-2)*[-1 0.75 0.75 1 1] (j==-1)*[-2 1 1 1 1] (j==-1)*(-1/64)+(j==-3)*(-81/64) ... + (j==1)*(-1/8) 0 (j==-1)*(-1/8) (j==1)*0.25 0 (j==-2)*(-2) 0 (j==-1)*0.25 ... + (j==1)*(-1/8) 0 (j==-1)*(-1/8) (j==-1)*(-1) ... %up to here - 0th order indirect + (j==-1)*[-0.5 3/8 0.25 0.5 0.5]+(j==1)*[1.5 0 -0.75 -1.5 -1.5] ... + (j==2)*[-2 1 1.5 2 2] (j==-2)*(-0.75) (j==1)*3/16+(j==3)*(-27/16)... + (j==1)*1.5 0 (j==1)*3 (j==-1)*(-1) (j==1)*(-3) (j==2)*(-4) 0 0 (j==1)*1.5 0 ... %up to here - 1st order indirect + (j==-1)*[-3/8 3/8 3/16 3/8 3/8]+(j==1)*[-1/8 -1/24 1/16 1/8 1/8] (j==2)*[3 0 -9/4 -3 -3] ... + (j==1)*[-1/8 1/16 -1/24 1/8 1/8]+(j==3)*[-27/8 27/16 27/8 27/8 27/8]... + (j==1)*[-1 0.5 0.5 0 1] (j==1)*[2 -1 -1 -1 -1] (j==1)*[-1 0.5 0.5 1 0]... + (j==-2)*(-2/3) (j==2)*(0.25)+(j==4)*(-8/3) (j==2)*3 0 ... + (j==1)*(-0.25) (j==-1)*(-0.75) (j==2)*6 0 (j==2)*(-6) 0 ... + (j==3)*(-27/4) (j==1)*(-0.25) 0 (j==2)*3 0 0 ... %up to here - 2nd order indirect + (j==-1)*(-1/3)+(j==1)*(-1/24) (j==2)*(-0.25) (j==1)*(-1/16)+(j==3)*(81/16) ... + (j==2)*(-1/6)+(j==4)*(-16/3) (j==1)*(-0.5) (j==2)*(-2) (j==1)*1 (j==2)*4 (j==1)*(-0.5) (j==2)*(-2) ... %up to here - 3rd order indirect + (j==-1)*(-125/384)+(j==1)*(-3/128) (j==2)*(-1/12) (j==1)*(-3/64)+(j==3)*(-27/64) ... + (j==2)*(-1/12)+(j==4)*(8) (j==3)*(-27/128)+(j==5)*(-3125/384) (j==1)*(-3/8) (j==2)*(-1) (j==3)*(-27/8) 0 ... + (j==1)*(0.75) (j==2)*(2) (j==3)*(27/4) 0 (j==1)*(-3/8) (j==2)*(-1) (j==3)*(-27/8) 0 0 0 ... %up to here - 4th order indirect + ]; + + %Indirect terms due to an internal perturber. The prefactor of 1/alph^2 + %comes from the definition of the disturbing function (Solar System + %Dynamics equation 6.45). + + + Terms.fI = 1/alph^2*[(j==1)*[-1 0.5 0.5 1 1 1/64 -0.25 1/64 -0.5 -0.5 -0.5 -0.5 0 0 -1] ... + +(j==-2)*[-1 0.75 0.75 1 1] (j==-1)*[-2 1 1 1 1] (j==-1)*(-1/64)+(j==-3)*(-81/64) ... + (j==1)*(-1/8) 0 (j==-1)*(-1/8) (j==1)*0.25 0 (j==-2)*(-2) 0 (j==-1)*0.25 ... + (j==1)*(-1/8) 0 (j==-1)*(-1/8) (j==-1)*(-1) ... %up to here - 0th order indirect + (j==-1)*[-2 1.5 1 2 2] ... + (j==0)*[1.5 -0.75 0 -1.5 -1.5]+(j==2)*[-0.5 0.25 3/8 0.5 0.5] ... + (j==0)*3/16+(j==-2)*(-27/16)... + (j==3)*(-0.75) 0 (j==0)*1.5 0 (j==-1)*(-4) 0 (j==2)*(-1) (j==0)*3 (j==0)*(-3) 0 (j==0)*1.5 ... %up to here - 1st order indirect + (j==-1)*[-27/8 27/8 27/16 27/8 27/8]+(j==1)*[-1/8 -1/24 1/16 1/8 1/8]... + (j==2)*[3 -9/4 0 -3 -3] (j==1)*[-1/8 1/16 -1/24 1/8 1/8]+(j==3)*[-3/8 3/16 3/8 3/8 3/8] ... + (j==1)*[-1 0.5 0.5 0 1] (j==1)*[2 -1 -1 -1 -1] (j==1)*[-1 0.5 0.5 1 0] ... + (j==0)*(0.25)+(j==-2)*(-8/3) (j==4)*(-2/3) 0 (j==0)*3 (j==1)*(-0.25) (j==-1)*(-27/4)... + 0 (j==0)*6 0 (j==0)*(-6) (j==3)*(-0.75) (j==1)*(-0.25) 0 0 (j==0)*3 0 ... %up to here - 2nd order indirect + (j==-1)*(-16/3)+(j==1)*(-1/6) (j==0)*(81/16)+(j==2)*(-1/16) (j==1)*(-0.25) ... + (j==2)*(-1/24)+(j==4)*(-1/3) (j==1)*(-2) (j==2)*(-0.5) (j==1)*4 (j==2)*1 (j==1)*(-2) (j==2)*(-0.5) ... %up to here - 3rd order indirect + (j==-1)*(-3125/384)+(j==1)*(-27/128) (j==0)*(8)+(j==2)*(-1/12) (j==1)*(-27/64)+(j==3)*(-3/64) ... + (j==2)*(-1/12) (j==3)*(-3/128)+(j==5)*(-125/384) (j==1)*(-27/8) (j==2)*(-1) (j==3)*(-3/8) 0 ... + (j==1)*(27/4) (j==2)*(2) (j==3)*(0.75) 0 (j==1)*(-27/8) (j==2)*(-1) (j==3)*(-3/8) 0 0 0 ... %up to here - 4th order indirect + ]; + + + Terms.Df = [Df1 Df2 Df2 Df3 Df3 Df4 Df5 Df6 Df7 Df7 Df7 Df7 Df8 Df8 Df9 ... + Df10 Df11 Df12 Df13 Df13 Df14 Df15 Df15 Df16 Df16 Df17 Df18 Df19 Df20 Df21 Df22 Df23 Df24 Df25 Df18 Df19 Df20 Df26... %0th order + Df27 Df28 Df29 Df30 Df30 Df31 Df32 Df33 Df34 Df34 Df35 Df36 Df37 Df38 Df39 Df40 Df41 Df42 Df43 Df44 Df37 Df38 ... %1st order + Df45 Df46 Df47 Df48 Df48 Df49 Df50 Df51 Df52 Df52 Df53 Df54 Df55 Df56 Df56 Df57 Df58 Df59 Df60 Df61 Df62 Df63 Df64 Df65 Df66 Df57 Df58 Df59 Df67 Df60 ... + Df68 Df69 Df70 Df71 Df72 Df73 Df74 Df75 Df76 Df77 Df78 Df79 Df80 Df70 Df71 Df81 ... %2nd order + Df82 Df83 Df84 Df85 Df86 Df87 Df88 Df89 Df86 Df87 ... %3rd order + Df90 Df91 Df92 Df93 Df94 Df95 Df96 Df97 Df98 Df99 Df100 Df101 Df102 Df95 Df96 Df97 Df103 Df102 Df98 ... %4th order + ]; + + Terms.DfE = Terms.fE/alph; %deriving alph*constant is the same as dividing by alph + Terms.DfI = Terms.fI*(-2/alph); %deriving alph^(-2)*constant is the same as multiplying by (-2/alph) + + end + + + + + + %% + Terms.j = ones(1,Terms.N)*j; + + if exist('AllTerms','var') + AllTerms = CatStructFields(AllTerms,Terms,2); + else + AllTerms = Terms; + end + + end + + %% one vectorized operation for all terms, including vectorizing on j + TermsInd = find((AllTerms.j~=0 | AllTerms.k~=0) & (AllTerms.A+AllTerms.Ap+AllTerms.B+AllTerms.Bp)<=MaxPower); + [dLambda1_pair,dLambda2_pair,dz1_pair,dz2_pair,du1_pair,du2_pair,da_oa1_pair,da_oa2_pair] = DisturbingFunctionMultiTermVariations(n1,n2,alph,Lambda1,Lambda2,mu1,mu2,AllTerms.j(TermsInd),AllTerms.k(TermsInd),e1,e2,pom1,pom2,I1,I2,s1,s2,Om1,Om2,AllTerms.A(TermsInd),AllTerms.Ap(TermsInd),AllTerms.B(TermsInd),AllTerms.Bp(TermsInd),AllTerms.C(TermsInd),AllTerms.Cp(TermsInd),AllTerms.D(TermsInd),AllTerms.Dp(TermsInd),AllTerms.f(TermsInd),AllTerms.fE(TermsInd),AllTerms.fI(TermsInd),AllTerms.Df(TermsInd),AllTerms.DfE(TermsInd),AllTerms.DfI(TermsInd)); + + + %% add values of dz and dLambda of the currently analyzed pair + dLambdaMatrix(:,j1) = dLambdaMatrix(:,j1)+dLambda1_pair; + dzMatrix(:,j1) = dzMatrix(:,j1)+dz1_pair; + duMatrix(:,j1) = duMatrix(:,j1)+du1_pair; + da_oaMatrix(:,j1) = da_oaMatrix(:,j1)+da_oa1_pair; + + dLambdaMatrix(:,j2) = dLambdaMatrix(:,j2)+dLambda2_pair; + dzMatrix(:,j2) = dzMatrix(:,j2)+dz2_pair; + duMatrix(:,j2) = duMatrix(:,j2)+du2_pair; + da_oaMatrix(:,j2) = da_oaMatrix(:,j2)+da_oa2_pair; + end +end + +%% +dLambda = dLambdaMatrix; +dz = dzMatrix; +du = duMatrix; +da_oa = da_oaMatrix; + + + + + + + +function [dLambda1_j,dz1_j,dLambda2_j,dz2_j] = VariationsPair(n1,n2,alph,Lambda1,Lambda2,mu1,mu2,j,z1,z2,e1,e2,pom1,pom2,u1,u2,s1,s2,Om1,Om2) +%This function is the algebraic "heart" of the code, as it calculates the +%integrals of the rhs of Lagrange's planetary equations to yield variations +%in the eccentricities and mean longitudes. + +%the following two flags, for adding additional terms, remain here for +%debugging and backward compatibility purposes, though for operational usage they will remain hard-coded as 1. +AddSecondOrderIncTerms = 1; +AddThirdOrder = 1; + + +%calculate laplace coefficients +Aj = LaplaceCoeffs(alph,0.5,j,0); +DAj = LaplaceCoeffs(alph,0.5,j,1); +D2Aj = LaplaceCoeffs(alph,0.5,j,2); +D3Aj = LaplaceCoeffs(alph,0.5,j,3); + +Ajp1 = LaplaceCoeffs(alph,0.5,j+1,0); +DAjp1 = LaplaceCoeffs(alph,0.5,j+1,1); +D2Ajp1 = LaplaceCoeffs(alph,0.5,j+1,2); +D3Ajp1 = LaplaceCoeffs(alph,0.5,j+1,3); + +Ajm1 = LaplaceCoeffs(alph,0.5,j-1,0); +DAjm1 = LaplaceCoeffs(alph,0.5,j-1,1); +D2Ajm1 = LaplaceCoeffs(alph,0.5,j-1,2); +D3Ajm1 = LaplaceCoeffs(alph,0.5,j-1,3); + +Ajm2 = LaplaceCoeffs(alph,0.5,j-2,0); +DAjm2 = LaplaceCoeffs(alph,0.5,j-2,1); +D2Ajm2 = LaplaceCoeffs(alph,0.5,j-2,2); +D3Ajm2 = LaplaceCoeffs(alph,0.5,j-2,3); + +f1 = 0.5*Aj; +f2 = -0.5*j^2*Aj+0.25*alph*DAj+1/8*alph^2*D2Aj; +f10 = 0.25*(2+6*j+4*j^2)*Ajp1-0.5*alph*DAjp1-0.25*alph^2*D2Ajp1; + +f27 = -j*Aj-0.5*alph*DAj; +f31 = (j-0.5)*Ajm1+0.5*alph*DAjm1; + +f45 = 1/8*(-5*j+4*j^2)*Aj+0.5*j*alph*DAj+1/8*alph^2*D2Aj; +f49 = 0.25*(-2+6*j-4*j^2)*Ajm1+0.25*alph*(2-4*j)*DAjm1-0.25*alph^2*D2Ajm1; +f53 = 1/8*(2-7*j+4*j^2)*Ajm2+1/8*(-2*alph+4*j*alph)*DAjm2+1/8*alph^2*D2Ajm2; + +Df1 = 0.5*DAj; +Df2 = -0.5*j^2*DAj+0.25*DAj+0.25*alph*D2Aj+1/4*alph*D2Aj+1/8*alph^2*D3Aj; +Df10 = 0.25*(2+6*j+4*j^2)*DAjp1-0.5*DAjp1-0.5*alph*D2Ajp1-0.5*alph*D2Ajp1-0.25*alph^2*D3Ajp1; + +Df27 = -j*DAj-0.5*DAj-0.5*alph*D2Aj; +Df31 = (j-0.5)*DAjm1+0.5*DAjm1+0.5*alph*D2Ajm1; + +Df45 = 1/8*(-5*j+4*j^2)*DAj+0.5*j*DAj+0.5*j*alph*D2Aj+1/4*alph*D2Aj+1/8*alph^2*D3Aj; +Df49 = 0.25*(-2+6*j-4*j^2)*DAjm1+0.25*(2-4*j)*DAjm1+0.25*alph*(2-4*j)*D2Ajm1-0.5*alph*D2Ajm1-0.25*alph^2*D3Ajm1; +Df53 = 1/8*(2-7*j+4*j^2)*DAjm2+1/8*(-2+4*j)*DAjm2+1/8*(-2*alph+4*j*alph)*D2Ajm2+1/4*alph*D2Ajm2+1/8*alph^2*D3Ajm2; + +%second order inclination-related and third-order inclination related terms +Bjm1 = LaplaceCoeffs(alph,1.5,j-1,0); +DBjm1 = LaplaceCoeffs(alph,1.5,j-1,1); +f57 = 0.5*alph*Bjm1; +Df57 = 0.5*Bjm1+0.5*alph*DBjm1; +f62 = -2*f57; +Df62 = -2*Df57; + +%define relative angle (as used at Hadden et al. 2019) +Psi = j*Lambda2-j*Lambda1; + +%define resonant angles +Lambda_j1 = j*Lambda2+(1-j)*Lambda1; +Lambda_j2 = j*Lambda2+(2-j)*Lambda1; +Lambda_j3 = j*Lambda2+(3-j)*Lambda1; + +%define "super mean motions" +nj1 = j*n2+(1-j)*n1; +nj2 = j*n2+(2-j)*n1; +nj3 = j*n2+(3-j)*n1; + +%da/a term - inner planet +if j~=0 + Order0Term = ((f1-alph*(j==1)+(f2+0.5*alph*(j==1))*(e1.^2+e2.^2)).*sin(Psi)... + +(f10-alph*(j==-2))*e1.*e2.*sin(Psi+pom2-pom1))/(-j)/((n2-n1).^2); +else + Order0Term = 0; +end + +Order1Term = (... + (f27+3/2*alph*(j==1)-0.5*alph*(j==-1)).*e1.*sin(Lambda_j1-pom1)+... + (f31-2*alph*(j==2)).*e2.*sin(Lambda_j1-pom2)... + )... + *(1-j)/((nj1).^2); + +Order2Term = (... + (f45-3/8*alph*(j==-1)-1/8*alph*(j==1))*e1.^2.*sin(Lambda_j2-2*pom1)+... + (f49+3*alph*(j==2)).*e1.*e2.*sin(Lambda_j2-pom1-pom2)... + +(f53-1/8*alph*(j==1)-27/8*alph*(j==3))*e2.^2.*sin(Lambda_j2-2*pom2)... + )... + .*(2-j)/((nj2).^2); + +Int_da_over_a_term = -3*alph*n1^2*(mu2/(1+mu1))*(Order0Term+Order1Term+Order2Term); + +%dR/dalpha terms - inner planet +if j~=0 + Order0Term = (... + (Df1-(j==1)+(Df2+0.5*(j==1))*(e1.^2+e2.^2)).*sin(Psi)... + +(Df10-(j==-2))*e1.*e2.*sin(Psi+pom2-pom1))/j./(n2-n1); +else + Order0Term = 0; +end + +Order1Term = (... + (Df27+3/2*(j==1)-0.5*(j==-1)).*e1.*sin(Lambda_j1-pom1)... + +(Df31-2*(j==2)).*e2.*sin(Lambda_j1-pom2)... + )... + ./(nj1); + +Order2Term = (... + (Df45-3/8*(j==-1)-1/8*(j==1)).*e1.^2.*sin(Lambda_j2-2*pom1)... + +(Df49+3*(j==2)).*e1.*e2.*sin(Lambda_j2-pom1-pom2)... + +(Df53-1/8*(j==1)-27/8*(j==3)).*e2.^2.*sin(Lambda_j2-2*pom2)... + )... + ./(nj2); + +Int_dR_dalph_term = -2*alph^2*n1*(mu2/(1+mu1))*(Order0Term+Order1Term+Order2Term); + +%e-derivative term - inner planet +if j~=0 + Order0Term = (... + (2*f2+alph*(j==1))*e1.*sin(Psi)+... + (f10-alph*(j==-2))*e2.*sin(Psi-pom1+pom2)... + )... + ./(j*(n2-n1)); +else + Order0Term = 0; +end + +Order1Term = (f27+3/2*alph*(j==1)-0.5*alph*(j==-1))*sin(Lambda_j1-pom1)./(nj1); + +Order2Term = (... + 2*(f45-3/8*alph*(j==-1)-1/8*alph*(j==1))*e1.*sin(Lambda_j2-2*pom1)+... + (f49+3*alph*(j==2))*e2.*sin(Lambda_j2-pom1-pom2)... + )... + ./(nj2); + +e_term = alph/2.*e1.*n1*(mu2/(1+mu1)).*(Order0Term+Order1Term+Order2Term); + + +dLambda1_j = Int_da_over_a_term+Int_dR_dalph_term+e_term; + +if AddSecondOrderIncTerms + %inclination terms of the inner planet + dl1Order2IncTerm1 = 3*mu2/(1+mu1)*alph*n1^2/(nj2)^2*(j-2)*(f57*s1.^2.*sin(Lambda_j2-2*Om1)+f62*s1.*s2.*sin(Lambda_j2-Om1-Om2)+f57*s2.^2.*sin(Lambda_j2-2*Om2)); + dl1Order2IncTerm2 = -2*mu2/(1+mu1)*alph^2*n1/(nj2)*(Df57*s1.^2.*sin(Lambda_j2-2*Om1)+Df62*s1.*s2.*sin(Lambda_j2-Om1-Om2)+Df57*s2.^2.*sin(Lambda_j2-2*Om2)); + dl1Order2IncTerm3 = mu2/(1+mu1)*alph*n1/nj2*(2*f57*s1.^2.*sin(Lambda_j2-2*Om1)+f62*s1.*s2.*sin(Lambda_j2-Om1-Om2)); + dLambda1_j = dLambda1_j+dl1Order2IncTerm1+dl1Order2IncTerm2+dl1Order2IncTerm3; +end + +%z variations - inner planet +if j~=0 + Order0Term = (... + ((2*f2+alph*(j==1)).*sin(Psi)+1i*j*(0.5*f1-0.5*alph*(j==1))*(-1)*cos(Psi)).*z1... + +(f10-alph*(j==-2))*exp(1i*(Psi)).*z2/1i... + )... + ./(j*(n2-n1)); +else + Order0Term = 0; +end + +Order1Term = (... + (f27+3/2*alph*(j==1)-0.5*alph*(j==-1))... + *(... + (1-0.5*e1.^2).*exp(1i*(Lambda_j1))./(1i)... + +0.5i*(j-1).*e1.^2.*exp(1i*pom1).*(-1).*cos(Lambda_j1-pom1)... + )... + +... + 0.5i*(f31-2*alph*(j==2))*(j-1)*e1.*e2.*exp(1i*pom1).*(-1).*cos(Lambda_j1-pom2)... + )... + ./(nj1); + +Order2Term = (... + 2*(f45-3/8*alph*(j==-1)-1/8*alph*(j==1)).*exp(1i*(Lambda_j2))/1i.*conj(z1)+... + (f49+3*alph*(j==2)).*exp(1i*(Lambda_j2))/1i.*conj(z2)... + )... + ./(nj2); + + +dz1_j = 1i*n1*alph*mu2/(1+mu1).*(Order0Term+Order1Term+Order2Term); + + +%da'/a' term - outer planet +if j~=0 + Order0Term = (... + -(f1-1/alph^2*(j==1)+(f2+0.5/alph^2*(j==1))*(e1.^2+e2.^2)).*sin(Psi)... + -(f10-1/alph^2*(j==-2))*e1.*e2.*sin(Psi+pom2-pom1)... + )... + /(-j)/((n2-n1).^2); +else + Order0Term = 0; +end + +Order1Term = (... + -(f27-2/alph^2*(j==-1))*e1.*sin(Lambda_j1-pom1)... + -(f31+3/2/alph^2*(j==0)-1/2/alph^2*(j==2)).*e2.*sin(Lambda_j1-pom2)... + )... + *(-j)/((nj1).^2); + +Order2Term = (... + -(f45-27/8/alph^2*(j==-1)-1/8/alph^2*(j==1))*e1.^2.*sin(Lambda_j2-2*pom1)... + -(f49+3/alph^2*(j==2)).*e1.*e2.*sin(Lambda_j2-pom1-pom2)... + -(f53-1/8/alph^2*(j==1)-3/8/alph^2*(j==3))*e2.^2.*sin(Lambda_j2-2*pom2)... + )... + .*(-j)./((nj2).^2); + +Int_da_over_a_term = -3*n2^2*(mu1/(1+mu2)).*(Order0Term+Order1Term+Order2Term); + +%dR'/dalpha terms - outer planet - corrected +if j~=0 + Order0Term = (... + (f1+alph*Df1+(j==1)/alph^2+(f2+alph*Df2-(j==1)/(2*alph^2)).*(e1.^2+e2.^2)).*sin(Psi)... + +(f10+alph*Df10+(j==-2)/alph^2).*e1.*e2.*sin(Psi+pom2-pom1)... + )... + /j./(n2-n1); +else + Order0Term = 0; +end + +Order1Term = (... + (f27+alph*Df27+2/(alph^2)*(j==-1)).*e1.*sin(Lambda_j1-pom1)... + +(f31+alph*Df31-3/alph^2*(j==0)+1/(2*alph^2)*(j==2)).*e2.*sin(Lambda_j1-pom2)... + )... + ./(nj1); + +Order2Term = (... + (f45+alph*Df45+27/8/alph^2*(j==-1)+1/8/alph^2*(j==1)).*e1.^2.*sin(Lambda_j2-2*pom1)... + +(f49+alph*Df49-3/alph^2*(j==2)).*e1.*e2.*sin(Lambda_j2-pom1-pom2)... + +(f53+alph*Df53+1/8/alph^2*(j==1)+3/8/alph^2*(j==3)).*e2.^2.*sin(Lambda_j2-2*pom2)... + )... + ./(nj2); + + +Int_dR_dalph_term = +2*n2*(mu1/(1+mu2))*(Order0Term+Order1Term+Order2Term); + +%e'-derivative term - outer planet +if j~=0 + Order0Term = (... + (2*f2+1/alph^2*(j==1))*e2.*sin(Psi)+... + (f10-1/alph^2*(j==-2))*e1.*sin(Psi+pom2-pom1)... + )... + ./(j*(n2-n1)); +else + Order0Term = 0; +end + +Order1Term = (f31+3/2/alph^2*(j==0)-1/2/alph^2*(j==2))*sin(Lambda_j1-pom2)... + ./(nj1); + +Order2Term = (... + (f49+3/alph^2*(j==2))*e1.*sin(Lambda_j2-pom1-pom2)+... + 2*(f53-1/8/alph^2*(j==1)-3/8/alph^2*(j==3))*e2.*sin(Lambda_j2-2*pom2)... + )... + ./(nj2); + +e_term = e2/2*n2*(mu1/(1+mu2)).*(Order0Term+Order1Term+Order2Term); + +dLambda2_j = (Int_da_over_a_term+Int_dR_dalph_term+e_term); + + +if AddSecondOrderIncTerms + %inclination terms of the outer planet + dl2Order2IncTerm1 = -3*mu1/(1+mu2)*n2^2/(nj2)^2*(j)*(f57*s1.^2.*sin(Lambda_j2-2*Om1)+f62*s1.*s2.*sin(Lambda_j2-Om1-Om2)+f57*s2.^2.*sin(Lambda_j2-2*Om2)); + dl2Order2IncTerm2 = +2*mu1/(1+mu2)*n2/(nj2)*((f57+alph*Df57)*s1.^2.*sin(Lambda_j2-2*Om1)+(f62+alph*Df62)*s1.*s2.*sin(Lambda_j2-Om1-Om2)+(f57+alph*Df57)*s2.^2.*sin(Lambda_j2-2*Om2)); + dl2Order2IncTerm3 = mu1/(1+mu2)*n2/nj2*(2*f57*s2.^2.*sin(Lambda_j2-2*Om2)+f62*s1.*s2.*sin(Lambda_j2-Om1-Om2)); + dLambda2_j = dLambda2_j+dl2Order2IncTerm1+dl2Order2IncTerm2+dl2Order2IncTerm3; +end + +%z' variations - outer planet +if j~=0 + Order0Term = (... + ((2*f2+1/alph^2*(j==1)).*sin(Psi)+1i*j*(-0.5*f1+0.5/alph^2*(j==1))*(-1)*cos(Psi)).*z2... + +(f10-1/alph^2*(j==-2))*exp(-1i*j*(Lambda2-Lambda1))/(-1i).*z1... + )... + ./(j*(n2-n1)); +else + Order0Term = 0; +end + +Order1Term = (... + (f31+3/2/alph^2*(j==0)-0.5/alph^2*(j==2))*(1-0.5*e2.^2).*exp(1i*(Lambda_j1))./(1i)... + -0.5i*j*(f27-1/alph^2*(j==-1)).*e1.*e2.*exp(1i*pom2).*(-1).*cos(Lambda_j1-pom1)... + -0.5i*j*(f31+3/2/alph^2*(j==0)-1/2/alph^2*(j==2))*e2.^2.*exp(1i*pom2).*(-1).*cos(Lambda_j1-pom2)... + )... + ./(nj1); + +Order2Term = (... + (f49+3/alph^2*(j==2))*exp(1i*(Lambda_j2))/1i.*conj(z1)+... + 2*(f53-1/8/alph^2*(j==1)-3/8/alph^2*(j==3)).*exp(1i*(Lambda_j2))/1i.*conj(z2)... + )... + ./(nj2); + + +dz2_j = 1i*n2*mu1/(1+mu2)*(Order0Term+Order1Term+Order2Term); + + + + +if AddThirdOrder + + %calculate Laplace Coefficients required for 3rd order resonant terms + Ajm3 = LaplaceCoeffs(alph,0.5,j-3,0); + DAjm3 = LaplaceCoeffs(alph,0.5,j-3,1); + D2Ajm3 = LaplaceCoeffs(alph,0.5,j-3,2); + D3Ajm3 = LaplaceCoeffs(alph,0.5,j-3,3); + + D4Aj = LaplaceCoeffs(alph,0.5,j,4); + D4Ajm1 = LaplaceCoeffs(alph,0.5,j-1,4); + D4Ajm2 = LaplaceCoeffs(alph,0.5,j-2,4); + D4Ajm3 = LaplaceCoeffs(alph,0.5,j-3,4); + + %calculate the coefficients from Murray&Dermott appendix B related to + %3rd order resonant terms + f82 = 1/48*(-26*j*Aj+30*j^2*Aj-8*j^3*Aj-9*alph*DAj+27*j*alph*DAj-12*j^2*alph*DAj+6*alph^2*D2Aj-6*j*alph^2*D2Aj-alph^3*D3Aj); + f83 = 1/16*(-9*Ajm1+31*j*Ajm1-30*j^2*Ajm1+8*j^3*Ajm1+9*alph*DAjm1-25*j*alph*DAjm1+12*j^2*alph*DAjm1-5*alph^2*D2Ajm1+6*j*alph^2*D2Ajm1+alph^3*D3Ajm1); + f84 = 1/16*(8*Ajm2-32*j*Ajm2+30*j^2*Ajm2-8*j^3*Ajm2-8*alph*DAjm2+23*j*alph*DAjm2-12*j^2*alph*DAjm2+4*alph^2*D2Ajm2-6*j*alph^2*D2Ajm2-alph^3*D3Ajm2); + f85 = 1/48*(-6*Ajm3+29*j*Ajm3-30*j^2*Ajm3+8*j^3*Ajm3+6*alph*DAjm3-21*j*alph*DAjm3+12*j^2*alph*DAjm3-3*alph^2*D2Ajm3+6*j*alph^2*D2Ajm3+alph^3*D3Ajm3); + + Df82 = 1/48*(-26*j*DAj+30*j^2*DAj-8*j^3*DAj-9*DAj-9*alph*D2Aj+27*j*DAj+27*j*alph*D2Aj-12*j^2*DAj-12*j^2*alph*D2Aj+12*alph*D2Aj+6*alph^2*D3Aj-12*j*alph*D2Aj-6*j*alph^2*D3Aj-3*alph^2*D3Aj-alph^3*D4Aj); + Df83 = 1/16*(-9*DAjm1+31*j*DAjm1-30*j^2*DAjm1+8*j^3*DAjm1+9*DAjm1+9*alph*D2Ajm1-25*j*DAjm1-25*j*alph*D2Ajm1+12*j^2*DAjm1+12*j^2*alph*D2Ajm1-10*alph*D2Ajm1-5*alph^2*D3Ajm1+12*j*alph*D2Ajm1+6*j*alph^2*D3Ajm1+3*alph^2*D3Ajm1+alph^3*D4Ajm1); + Df84 = 1/16*(8*DAjm2-32*j*DAjm2+30*j^2*DAjm2-8*j^3*DAjm2-8*DAjm2-8*alph*D2Ajm2+23*j*DAjm2+23*j*alph*D2Ajm2-12*j^2*DAjm2-12*j^2*alph*D2Ajm2+8*alph*D2Ajm2+4*alph^2*D3Ajm2-12*j*alph*D2Ajm2-6*j*alph^2*D3Ajm2-3*alph^2*D3Ajm2-alph^3*D4Ajm2); + Df85 = 1/48*(-6*DAjm3+29*j*DAjm3-30*j^2*DAjm3+8*j^3*DAjm3+6*DAjm3+6*alph*D2Ajm3-21*j*DAjm3-21*j*alph*D2Ajm3+12*j^2*DAjm3+12*j^2*alph*D2Ajm3-6*alph*D2Ajm3-3*alph^2*D3Ajm3+12*j*alph*D2Ajm3+6*j*alph^2*D3Ajm3+3*alph^2*D3Ajm3+alph^3*D4Ajm3); + + + D2Bjm1 = LaplaceCoeffs(alph,1.5,j-1,2); + Bjm2 = LaplaceCoeffs(alph,1.5,j-2,0); + DBjm2 = LaplaceCoeffs(alph,1.5,j-2,1); + D2Bjm2 = LaplaceCoeffs(alph,1.5,j-2,2); + + + f86 = 0.75*alph*Bjm1-0.5*j*alph*Bjm1-0.25*alph^2*DBjm1; + f87 = 0.5*j*alph*Bjm2+0.25*alph^2*DBjm2; + f88 = -2*f86; + f89 = -2*f87; + Df86 = 0.75*Bjm1+0.75*alph*DBjm1 -0.5*j*Bjm1-0.5*j*alph*DBjm1 -0.5*alph*DBjm1-0.25*alph^2*D2Bjm1; + Df87 = 0.5*j*Bjm2+0.5*j*alph*DBjm2+ 0.5*alph*DBjm2+0.25*alph^2*D2Bjm2; + Df88 = -2*Df86; + Df89 = -2*Df87; + + %calculate the 3rd order resonant terms contribution to the forced eccentricity + dz1Order3Term1 = mu2/(1+mu1)*alph*n1/(nj3)*exp(1i*pom1).*(3*f82*e1.^2.*exp(1i*(Lambda_j3-3*pom1)) +2*f83.*e1.*e2.*exp(1i*(Lambda_j3-2*pom1-pom2))+1*f84*e2.^2.*exp(1i*(Lambda_j3-pom1-2*pom2))); + dz2Order3Term1 = mu1/(1+mu2) *n2/(nj3)*exp(1i*pom2).*(1*f83*e1.^2.*exp(1i*(Lambda_j3-2*pom1-pom2))+2*f84*e1.*e2.*exp(1i*(Lambda_j3-pom1-2*pom2))+3*(f85-(j==4)/(3*alph^2))*e2.^2.*exp(1i*(Lambda_j3-3*pom2))); + + %add the contribution of the terms involving both e and I (or s in + %Murray&Dermott formulae) - these with the coefficients f86-f89 + dz1Order3Term2 = mu2/(1+mu1)*alph*n1/(nj3).*(f86*s1.^2.*exp(1i*(Lambda_j3-2*Om1))+f88*s1.*s2.*exp(1i*(Lambda_j3-Om1-Om2))+f86*s2.^2.*exp(1i*(Lambda_j3-2*Om2))); + dz2Order3Term2 = mu1/(1+mu2) *n2/(nj3).*(f87*s1.^2.*exp(1i*(Lambda_j3-2*Om1))+f89*s1.*s2.*exp(1i*(Lambda_j3-Om1-Om2))+f87*s2.^2.*exp(1i*(Lambda_j3-2*Om2))); + + %calculate the 3rd order resonant terms contribution to the forced + %mean longitude + + %dn/n contribution + dl1Order3Term1 = 3*(j-3)*mu2/(1+mu1)*alph*n1^2/(nj3)^2*(f82*e1.^3.*sin(Lambda_j3-3*pom1)+f83*e1.^2.*e2.*sin(Lambda_j3-2*pom1-pom2)+f84*e1.*e2.^2.*sin(Lambda_j3-pom1-2*pom2)+(f85-(j==4)*16/3*alph)*e2.^3.*sin(Lambda_j3-3*pom2)); + dl2Order3Term1 = 3*(-j)*mu1/(1+mu2)*n2^2/(nj3)^2*(f82*e1.^3.*sin(Lambda_j3-3*pom1)+f83*e1.^2.*e2.*sin(Lambda_j3-2*pom1-pom2)+f84*e1.*e2.^2.*sin(Lambda_j3-pom1-2*pom2)+(f85-(j==4)/(3*alph^2))*e2.^3.*sin(Lambda_j3-3*pom2)); + + %dR/da and dR/da' contribution + dl1Order3Term2 = -2*mu2/(1+mu1)*alph^2*n1/(nj3)*(Df82*e1.^3.*sin(Lambda_j3-3*pom1)+Df83*e1.^2.*e2.*sin(Lambda_j3-2*pom1-pom2)+Df84*e1.*e2.^2.*sin(Lambda_j3-pom1-2*pom2)+(Df85-16/3*(j==4))*e2.^3.*sin(Lambda_j3-3*pom2)); + dl2Order3Term2 = 2*mu1/(1+mu2)*n2/(nj3)*((f82+alph*Df82)*e1.^3.*sin(Lambda_j3-3*pom1)+(f83+alph*Df83)*e1.^2.*e2.*sin(Lambda_j3-2*pom1-pom2)+(f84+alph*Df84)*e1.*e2.^2.*sin(Lambda_j3-pom1-2*pom2)+(f85+alph*Df85+2/(3*alph^3)*(j==4))*e2.^3.*sin(Lambda_j3-3*pom2)); + + %add the contribution of the terms involving both e and I (or s in + %Murray&Dermott formulae) - these with the coefficients f86-f89 + + %dn/n contribution + dl1Order3Term3 = 3*(j-3)*mu2/(1+mu1)*alph*n1^2/(nj3)^2*(f86*e1.*s1.^2.*sin(Lambda_j3-pom1-2*Om1) +f87*e2.*s1.^2.*sin(Lambda_j3-pom2-2*Om1) +f88*e1.*s1.*s2.*sin(Lambda_j3-pom1-Om1-Om2) +f89*e2.*s1.*s2.*sin(Lambda_j3-pom2-Om1-Om2) + f86*e1.*s2.^2.*sin(Lambda_j3-pom1-2*Om2) +f87*e2.*s2.^2.*sin(Lambda_j3-pom2-2*Om2)); + dl2Order3Term3 = 3*(-j) *mu1/(1+mu2)*n2^2 ./(nj3)^2*(f86*e1.*s1.^2.*sin(Lambda_j3-pom1-2*Om1) +f87*e2.*s1.^2.*sin(Lambda_j3-pom2-2*Om1) +f88*e1.*s1.*s2.*sin(Lambda_j3-pom1-Om1-Om2) +f89*e2.*s1.*s2.*sin(Lambda_j3-pom2-Om1-Om2) + f86*e1.*s2.^2.*sin(Lambda_j3-pom1-2*Om2) +f87*e2.*s2.^2.*sin(Lambda_j3-pom2-2*Om2)); + + %dR/da and dR/da' contribution + dl1Order3Term4 = -2*mu2/(1+mu1)*alph^2*n1/(nj3)* (Df86*e1.*s1.^2.*sin(Lambda_j3-pom1-2*Om1) + Df87*e2.*s1.^2.*sin(Lambda_j3-pom2-2*Om1) +Df88*e1.*s1.*s2.*sin(Lambda_j3-pom1-Om1-Om2) +Df89*e2.*s1.*s2.*sin(Lambda_j3-pom2-Om1-Om2) + Df86*e1.*s2.^2.*sin(Lambda_j3-pom1-2*Om2) + Df87*e2.*s2.^2.*sin(Lambda_j3-pom2-2*Om2)); + dl2Order3Term4 = 2*mu1/(1+mu2)*n2/(nj3)* ((f86+alph*Df86)*e1.*s1.^2.*sin(Lambda_j3-pom1-2*Om1) + (f87+alph*Df87)*e2.*s1.^2.*sin(Lambda_j3-pom2-2*Om1) +(f88+alph*Df88)*e1.*s1.*s2.*sin(Lambda_j3-pom1-Om1-Om2) +(f89+alph*Df89)*e2.*s1.*s2.*sin(Lambda_j3-pom2-Om1-Om2) + (f86+alph*Df86)*e1.*s2.^2.*sin(Lambda_j3-pom1-2*Om2) + (f87+alph*Df87)*e2.*s2.^2.*sin(Lambda_j3-pom2-2*Om2)); + + %sum all the contributions to delta-lambda + dl1Order3 = dl1Order3Term1+dl1Order3Term2+dl1Order3Term3+dl1Order3Term4; + dl2Order3 = dl2Order3Term1+dl2Order3Term2+dl2Order3Term3+dl2Order3Term4; + + %sum all the contributions to delta-z + dz1Order3 = dz1Order3Term1+dz1Order3Term2; + dz2Order3 = dz2Order3Term1+dz2Order3Term2; + + %add the 3rd order contribution to the total forced eccentricities and + %mean longitudes + dLambda1_j = dLambda1_j+dl1Order3; + dz1_j = dz1_j+dz1Order3; + dLambda2_j = dLambda2_j+dl2Order3; + dz2_j = dz2_j+dz2Order3; + +end diff --git a/AnalyticLC/ForcedMeanLongitudes3PlanetsCrossTerms.m b/AnalyticLC/ForcedMeanLongitudes3PlanetsCrossTerms.m new file mode 100644 index 0000000..d5eacb6 --- /dev/null +++ b/AnalyticLC/ForcedMeanLongitudes3PlanetsCrossTerms.m @@ -0,0 +1,82 @@ +function dLambda = ForcedMeanLongitudes3PlanetsCrossTerms(P,Lambda,mu,z) +% ForcedMeanLongitudes3Planets: calculate the forced mean longitudes of a set of +% planets due to the cross-terms among different resonances. + +%Inputs: P - orbital periods, Lambda - mean longitudes, mu - planets/star +%mass ratio, z - free eccentricities (if they are not given, zero values +%are assumed). Labmda and z can be given as columns of multiple values. + +% Yair Judkovsky, 2.11.2020 + +if ~exist('z','var') + z = zeros(size(P)); +end + +%mean motions +n = 2*pi./P; + +%absolute values and angles of eccentricities +pom = angle(z); +e = abs(z); +zx = real(z); +zy = imag(z); + +%allocate space for 3d matrices that will hold the values of dLambda +%according to these dimensios: orders, data points, number of planets +dLambda = zeros(size(z,1),length(P)); + +%go over all possible planets triplets +for j1 = 1:length(P) + for j2 = (j1+1):length(P) + for j3 = (j2+1):length(P) + + %calculate the nearest 1st order MMR for each pair in the + %triplet + j = round(P(j2)/(P(j2)-P(j1))); + k = round(P(j3)/(P(j3)-P(j2))); + nj = j*n(j2)+(1-j)*n(j1); + nk = k*n(j3)+(1-k)*n(j2); + alph12 = (n(j2)/n(j1))^(2/3)*((1+mu(j1))/(1+mu(j2)))^(1/3); + alph23 = (n(j3)/n(j2))^(2/3)*((1+mu(j2))/(1+mu(j3)))^(1/3); + + %longitudes of conjunction + Lambda_j = j*Lambda(:,j2)+(1-j)*Lambda(:,j1); + Lambda_k = k*Lambda(:,j3)+(1-k)*Lambda(:,j2); + + %laplace coefficients + Ak = LaplaceCoeffs(alph23,0.5,k,0); + DAk = LaplaceCoeffs(alph23,0.5,k,1); + Ajm1 = LaplaceCoeffs(alph12,0.5,j-1,0); + DAjm1 = LaplaceCoeffs(alph12,0.5,j-1,1); + D2Ajm1 = LaplaceCoeffs(alph12,0.5,j-1,2); + f27k = -k*Ak-0.5*alph23*DAk; + f31j = (j-0.5)*Ajm1+0.5*alph12*DAjm1; + + Df31j = (j-0.5)*DAjm1+0.5*DAjm1+0.5*alph12*D2Ajm1; + + F=sqrt(1); + %calculate and add to the dlambda of the intermediate planet + AlphaFactor = -3/2*f27k*(f31j-(j==2)/(2*alph12^2))*alph23; + MassFactor = mu(j1)*mu(j3)/(1+mu(j2))^2; + SinTerm1 = (j/nk+(1-k)/nj)*sin(Lambda_j+Lambda_k-2*pom(:,j2))/(nj+nk)^2; + SinTerm2 = (j/nk-(1-k)/nj)*sin(Lambda_j-Lambda_k)/(nj-nk)^2; + dLambda(:,j2) = dLambda(:,j2)+F*AlphaFactor*MassFactor*n(j2)^3*(SinTerm1+SinTerm2); + + %calculate and add to the dlambda of the innermost planet + AlphaFactor = 3/2*f27k*(f31j-2*alph12*(j==2))*alph12*alph23; + MassFactor = mu(j2)*mu(j3)/(1+mu(j1))/(1+mu(j2)); + SinTerm1 = sin(Lambda_j+Lambda_k-2*pom(:,j2))/(nj+nk)^2; + SinTerm2 = sin(Lambda_j-Lambda_k)/(nj-nk)^2; + dLambda(:,j1) = dLambda(:,j1)+F*AlphaFactor*MassFactor*n(j1)^2*n(j2)*(j-1)/nk*(SinTerm1+SinTerm2); + + %calculate and add to the dlambda of the outermost planet + AlphaFactor = -3/2*f27k*(f31j-(j==2)/(2*alph12^2)); + MassFactor = mu(j1)*mu(j2)/(1+mu(j2))/(1+mu(j3)); + SinTerm1 = sin(Lambda_j+Lambda_k-2*pom(:,j2))/(nj+nk)^2; + SinTerm2 = -sin(Lambda_j-Lambda_k)/(nj-nk)^2; + dLambda(:,j3) = dLambda(:,j3)+F*AlphaFactor*MassFactor*n(j3)^2*n(j2)*k/nj*(SinTerm1+SinTerm2); + + + end + end +end diff --git a/AnalyticLC/LaplaceCoeffs.m b/AnalyticLC/LaplaceCoeffs.m new file mode 100644 index 0000000..28a11fd --- /dev/null +++ b/AnalyticLC/LaplaceCoeffs.m @@ -0,0 +1,166 @@ +function Output = LaplaceCoeffs(varargin) +% LaplaceCoeff: get a Laplace coefficient (definition: Solar System +% Dynamics, eq. 6.67) + +%Inputs: 'init' - calculate a table of the laplace coefficients and save it +% in ~/Matlab/Data/LaplaceCoeffs.mat +% 'load' = return all the table and the row/col values in a +% structure +% alph,s,j = return the laplace coefficient for the given alph,s,j. +% alph can be a vector +% alph,s,j,D = return the D order derivative of the laplace +% coefficient for the given alph, s, j. + +% Yair Judkovsky, 6.4.2020 + +persistent jVec sVec DVec alphaVec bTable Nj Ns ND + + +%return a specific laplace coeff +if isnumeric(varargin{1}) + + alph = varargin{1}; + s = varargin{2}; + j = abs(varargin{3}); %b_s_j = b_s__-j, Solar System Dynamics eq. 6.69 + + %a recursive path of the function for the case there's a 4th input - + %the laplace coeff derivative. + if length(varargin)>3 + D = varargin{4}; %derivative order + else + D = 0; + end + + jInd = j+1; + sInd = s+0.5; + DInd = D+1; + + %if the requested derivative is too high, calculate by recursive method + if DInd>ND + %2nd and higher derivatives of the coeff, Solar System Dynamics eq. 6.71 + Output = s*(LaplaceCoeffs(alph,s+1,j-1,D-1)-2*alph.*LaplaceCoeffs(alph,s+1,j,D-1)+LaplaceCoeffs(alph,s+1,j+1,D-1)-2*(D-1)*LaplaceCoeffs(alph,s+1,j,D-2)); + return + end + + row = sub2ind([Nj,Ns,ND],jInd,sInd,DInd); + + %find the column corresponding to the alpha given, currently by finding the + %closest value in the table + F = (length(alphaVec)); + col = round(alph*F); + if col==0 + col = 1; + end + + %return the value from the table + Output = bTable(row,col); + return +end + +%perform the calculation. Table will be made such that each row +%corresponds to some j,s combination, and each column corresponds to an +%alpha value. +if isequal(varargin{1},'init') + + %define the integrand function + laplace_integrand=@(j1,s1,alph1) (@(phi) (cos(j1.*phi)./(1-2.*alph1.*cos(phi)+alph1.^2).^s1)); + + warning off + + %set range of j, s, D, alpha. Can also be read from varargin if + %exists, where the 2,3,4,5 arguments will be Ns, Nj,ND and dalpha + Ns =3; + Nj = 5; + ND = 4; + dalpha = 1e-2; + %set integration tolerance + AbsTol = 1e-5; + Path = '~/Matlab/Data/LaplaceCoeffs.mat'; + + %read inputs from varargin, if they exist + if length(varargin)==5 + Ns = varargin{2}; + Nj = varargin{3}; + ND = varargin{4}; + dalpha = varargin{5}; + end + + if length(varargin)>5 + AbsTol = varargin{6}; + end + + if length(varargin)>6 + Path = varargin{7}; + end + + %set a vector of s and j values + sVec = (1:Ns)-0.5; + jVec = (1:Nj)-1; + DVec = (1:ND)-1; + + %define the alpha vector from dalpha + alphaVec = dalpha:dalpha:1; + + %pre-allocate + bTable = zeros(Nj*Ns*ND,length(alphaVec)); + + for jInd = 1:length(jVec) + for sInd = 1:length(sVec) + for DInd = 1:length(DVec) + for alphaInd = 1:length(alphaVec) + + row = sub2ind([Nj,Ns,ND],jInd,sInd,DInd); + + val = CalculateSingleCoeff(jVec(jInd),sVec(sInd),alphaVec(alphaInd),DVec(DInd),AbsTol,laplace_integrand); + + bTable(row,alphaInd)=val; + + end + end + end + end + save(Path,'jVec','sVec','alphaVec','DVec','bTable','Nj','Ns','ND'); + disp('calculated laplace coeffs'); + Output = []; + return +end + +%load the table including the persistent variables for future use. Also +%return the persistent variables in a structure as an output. +if isequal(varargin{1},'load') + Path = '~/Matlab/Data/LaplaceCoeffs.mat'; + if nargin>1 + Path = varargin{2}; + end + load(Path); + Output = load(Path); + return +end + + +% if D==0 +% Output = LaplaceCoeffs(alph,s,j); +% elseif D==1 %first derivative of the coeff, Solar System Dynamics eq. 6.70 +% Output = s*(LaplaceCoeffs(alph,s+1,j-1)-2*alph.*LaplaceCoeffs(alph,s+1,j)+LaplaceCoeffs(alph,s+1,j+1)); +% else %2nd and higher derivatives of the coeff, Solar System Dynamics eq. 6.71 +% Output = s*(LaplaceCoeffs(alph,s+1,j-1,D-1)-2*alph.*LaplaceCoeffs(alph,s+1,j,D-1)+LaplaceCoeffs(alph,s+1,j+1,D-1)-2*(D-1)*LaplaceCoeffs(alph,s+1,j,D-2)); +% end +% return + +function val = CalculateSingleCoeff(j,s,alph,D,AbsTol,laplace_integrand) + +% %define the integrand function +% laplace_integrand=@(j1,s1,alph1) (@(phi) (cos(j1.*phi)./(1-2.*alph1.*cos(phi)+alph1.^2).^s1)); + +if D==0 + val = 1/pi*integral(laplace_integrand(j,s,alph),0,2*pi,'AbsTol',AbsTol); +elseif D==1 %first derivative of the coeff, Solar System Dynamics eq. 6.70 + val = s*(CalculateSingleCoeff(j-1,s+1,alph,D-1,AbsTol,laplace_integrand)... + -2*alph.*CalculateSingleCoeff(j,s+1,alph,D-1,AbsTol,laplace_integrand)... + +CalculateSingleCoeff(j+1,s+1,alph,D-1,AbsTol,laplace_integrand)); +else %2nd and higher derivatives of the coeff, Solar System Dynamics eq. 6.71 + val = s*(CalculateSingleCoeff(j-1,s+1,alph,D-1,AbsTol,laplace_integrand)... + -2*alph.*CalculateSingleCoeff(j,s+1,alph,D-1,AbsTol,laplace_integrand)... + +CalculateSingleCoeff(j+1,s+1,alph,D-1,AbsTol,laplace_integrand)... + -2*(D-1)*CalculateSingleCoeff(j,s+1,alph,D-2,AbsTol,laplace_integrand)); +end diff --git a/AnalyticLC/MALookupTable.m b/AnalyticLC/MALookupTable.m new file mode 100644 index 0000000..92a88cf --- /dev/null +++ b/AnalyticLC/MALookupTable.m @@ -0,0 +1,105 @@ +function Output = MALookupTable(varargin) +% MALookupTable: generate a look-up table for Mandel-Agol model values, or +% use the table and read from it + +%Inputs: 'init' - calculate a table of the MA values for the given rVec, +%zVec, u1, u2 +% 'load' = return all the table and the row/col values in a +% structure +% r,z = return the Mandel-Agol model values for the given r,z. +% u1 and u2 which will be used are the ones that were used for +% generating the table. z can be a vector, r must be a scalar. + +% Yair Judkovsky, 25.10.2020 + +persistent rVec zVec u1 u2 R Z MATable + +%return a specific Mandel-Agol model value +if isnumeric(varargin{1}) + + %get r, z values from the input + r = varargin{1}; + z = varargin{2}; + Requested_u1 = varargin{3}; + Requested_u2 = varargin{4}; + + %if no table exists, or if the requested u1,u2 do not match the table, just calculate the values + if isempty(MATable) || Requested_u1~=u1 || Requested_u2~=u2 + Output = occultquadVec(r,z,Requested_u1,Requested_u2); + return + end + + Output = ones(size(z)); + + %old method - plug in "1" wherever outside the planet radius range +% Output(z<1+r) = interp2(R,Z,MATable,z(z<1+r),r,'linear',1); + + %new method - outside the radius range calculate using occultquadVec. This + %is useful because now the radius range can be reduce, yielding a fast + %calculation for most of the points + ind = r>min(rVec) & r3 +asin0 = asin(sqrt(((1+ex).^2./(aor.^2.*(1-ex.^2-ey.^2).^2)-sin(I).^2.*sin(Omega).^2)./(1-sin(I).^2.*sin(Omega).^2))); +T=1./n*2.*((1-ex.^2-ey.^2).^(3/2))./((1+ex).^2).*asin0; +asin2 = asin(sqrt(((1+r).^2.*(1+ex).^2./(aor.^2.*(1-ex.^2-ey.^2).^2)-sin(I).^2.*sin(Omega).^2)./(1-sin(I).^2.*sin(Omega).^2))); +asin1 = asin(sqrt(((1-r).^2.*(1+ex).^2./(aor.^2.*(1-ex.^2-ey.^2).^2)-sin(I).^2.*sin(Omega).^2)./(1-sin(I).^2.*sin(Omega).^2))); +Tau=1./n.*((1-ex.^2-ey.^2).^(3/2))./((1+ex).^2).*(asin2-asin1); +end \ No newline at end of file diff --git a/AnalyticLC/fooo.m b/AnalyticLC/fooo.m new file mode 100644 index 0000000..3287ad6 --- /dev/null +++ b/AnalyticLC/fooo.m @@ -0,0 +1 @@ +t = 0:100:0.02 diff --git a/AnalyticLC/fooo.octave b/AnalyticLC/fooo.octave new file mode 100644 index 0000000..d58c058 --- /dev/null +++ b/AnalyticLC/fooo.octave @@ -0,0 +1,16 @@ +P = [9, 19, 32] +Tmid0 = [0, 10, 15] +ror = [0.03, 0.12, 0.05] +aor = 25 +mu = [0.02, 0.05, 0.01] +ex0 = [0, 0, 0] +ey0 = [0, 0, 0] +I0 = [90, 90, 90] +Omega0 = [0, 0, 0] +t = [0, 9, 19, 22, 32] +u1 = 0.27 +u2 = 0.13 + +LC = AnalyticLC(P, Tmid0, ror, aor, mu, ex0, ey0, I0, Omega0, t, u1, u2) + +print LC diff --git a/AnalyticLC/test.octave b/AnalyticLC/test.octave new file mode 100644 index 0000000..d58c058 --- /dev/null +++ b/AnalyticLC/test.octave @@ -0,0 +1,16 @@ +P = [9, 19, 32] +Tmid0 = [0, 10, 15] +ror = [0.03, 0.12, 0.05] +aor = 25 +mu = [0.02, 0.05, 0.01] +ex0 = [0, 0, 0] +ey0 = [0, 0, 0] +I0 = [90, 90, 90] +Omega0 = [0, 0, 0] +t = [0, 9, 19, 22, 32] +u1 = 0.27 +u2 = 0.13 + +LC = AnalyticLC(P, Tmid0, ror, aor, mu, ex0, ey0, I0, Omega0, t, u1, u2) + +print LC diff --git a/FittingAndMath/CatStructFields.m b/FittingAndMath/CatStructFields.m new file mode 100644 index 0000000..afb5c65 --- /dev/null +++ b/FittingAndMath/CatStructFields.m @@ -0,0 +1,6 @@ +function S = CatStructFields(S, T, dim) +fields = fieldnames(S); +for k = 1:numel(fields) + aField = fields{k}; % EDIT: changed to {} + S.(aField) = cat(dim, S.(aField), T.(aField)); +end \ No newline at end of file diff --git a/FittingAndMath/Vector1stOrderDE.m b/FittingAndMath/Vector1stOrderDE.m new file mode 100644 index 0000000..cc8230b --- /dev/null +++ b/FittingAndMath/Vector1stOrderDE.m @@ -0,0 +1,37 @@ +function xt = Vector1stOrderDE(A,x0,t) +%Vector1stOrderDE: - solves the vector equation dx/dt=Ax + +%inputs: A - matrix of the equation, x0 - x for t=0, t - time axis to +%evaluate values at. derivation is given at: +%https://www.unf.edu/~mzhan/chapter4.pdf + +%Yair Judkovsky, 3.10.2020 + +if nargin==0 +A=[.1,.01;-.002,-.002].*i; +x0=[.3,0.0001]; +t=linspace(0,100,100000); +end +nt=length(t); +if any(~isfinite(A(:))), xt = 123.123; return, end +[ev,lam]=eig(A); +nev=length(diag(lam)); +phi0i=inv(ev); +% add singleton dimension as the first (time) index and replicate +% eigenvectors along that dimension. +evp=repmat(reshape(ev,[1,nev,nev]),[length(t),1,1]); + +lamt=reshape(exp(diag(lam)*t).',length(t),1,nev); +lamt=repmat(lamt,1,nev,1); +phit=lamt.*evp; + +xt=reshape(reshape(phit,nev*nt,nev)*(phi0i*x0),nt,nev); + +if nargin==0 +figure; +plot(xt); +hold on +plot(real(x0),imag(x0),'o'); +axis equal +end + diff --git a/FittingAndMath/fit_linear.m b/FittingAndMath/fit_linear.m new file mode 100644 index 0000000..c5e195e --- /dev/null +++ b/FittingAndMath/fit_linear.m @@ -0,0 +1,45 @@ +function [a,b,da,db]=fit_linear(x,y,varargin) +%fit_linear: 1st order linear regression. + +%This function fits a dataset of x and y values to the linear equation +%y=ax+b. it is optional for the user to insert also y error values for the +%estimation of a and b errors. +% + +%Yair Judkovsky, 21.8.2013 + +if nargin==2 %case of no error input - assuming that all errors are equal + dy = ones(size(y)); +else + dy = varargin{1}; +end + +dy2=dy.^2; +x2=x.^2; +xy=x.*y; + +S=sum(1./dy2); +Sx=sum(x./dy2); +Sy=sum(y./dy2); +Sxx=sum(x2./dy2); +Sxy=sum(xy./dy2); +D=S*Sxx-Sx^2; + +a=(S*Sxy-Sx*Sy)/D; +b=(Sxx*Sy-Sxy*Sx)/D; + +if nargin==2 %no error given - need to estimate it + Ytheo=a*x+b; + dy2=mean((Ytheo-y).^2)*ones(size(y)); + S=sum(1./dy2); + Sx=sum(x./dy2); + Sxx=sum(x2./dy2); + D=S*Sxx-Sx^2; + da=sqrt(S/D); + db=sqrt(Sxx/D); +else + da=sqrt(S/D); + db=sqrt(Sxx/D); +end + +return diff --git a/FittingAndMath/horizontal.m b/FittingAndMath/horizontal.m new file mode 100644 index 0000000..47e812a --- /dev/null +++ b/FittingAndMath/horizontal.m @@ -0,0 +1,9 @@ +function horzvec=horizontal(vec) +%reshape input vector into row vector + +if min(size(vec))>1 + horzvec=[]; + return +end + +horzvec=reshape(vec,1,length(vec)); \ No newline at end of file diff --git a/FittingAndMath/range.m b/FittingAndMath/range.m new file mode 100644 index 0000000..40b155a --- /dev/null +++ b/FittingAndMath/range.m @@ -0,0 +1,3 @@ +function res = range(X) +res=max(X(:))-min(X(:)); +return \ No newline at end of file diff --git a/FittingAndMath/vertical.m b/FittingAndMath/vertical.m new file mode 100644 index 0000000..fd2763a --- /dev/null +++ b/FittingAndMath/vertical.m @@ -0,0 +1,9 @@ +function vertvec=vertical(vec) +%reshape input vector into row vector + +if min(size(vec))>1 + vertvec=[]; + return +end + +vertvec=reshape(vec,length(vec),1); \ No newline at end of file diff --git a/FluxCalculation/ellecVec.m b/FluxCalculation/ellecVec.m new file mode 100644 index 0000000..90009a0 --- /dev/null +++ b/FluxCalculation/ellecVec.m @@ -0,0 +1,17 @@ +function ellec=ellecVec(k) +% Computes polynomial approximation for the complete elliptic +% integral of the second kind (Hasting's approximation): +m1=1-k.*k; +m1(m1size(z,1) % z is a column vector + z=z.'; +end; + +%set zeros to lambdad, lambda3, 3tad, kap0, kap1 +lambdad=zeros(size(z)); +lambdae=lambdad; +etad=lambdad; +kap0=lambdad; +kap1=lambdad; + +%if p is close to 0.5, make it 0.5 +if abs(p-0.5)<1e-3 + p=0.5; +end; + +omega=1-u1/3-u2/6; +x1=(p-z).^2; +x2=(p+z).^2; +x3=p.^2-z.^2; + +% the source is unocculted: +% Table 3, I. +i=z>=1+p; +if any(i) + lambdad(i)=0; + etad(i)=0; + lambdae(i)=0; +end; + +% the source is partly occulted and the occulting object crosses the limb: +% Equation (26): +i=z>=abs(1-p) & z<=1+p; +if any(i) + kap1(i)=acos(min((1-p^2+z(i).^2)/2./z(i),1)); + kap0(i)=acos(min((p*p+z(i).^2-1)/2/p./z(i),1)); + lambdae(i)=p^2*kap0(i)+kap1(i); + lambdae(i)=(lambdae(i)-0.5*sqrt(max(4*z(i).^2-(1+z(i).^2-p^2).^2,0)))/pi; +end; + +% the occulting object transits the source star (but doesn't +% completely cover it): +i=z<=1-p; +if any(i) + lambdae(i)=p^2; +end; + +% the edge of the occulting star lies at the origin- special +% expressions in this case: + +i=find(abs(z-p)<1e-4*(z+p)); +if ~isempty(i) +% Table 3, Case V.: + j=z(i)>=0.5; + if any(j) + %lam=0.5*pi; + q=0.5/p; + Kk=ellkVec(q); + Ek=ellecVec(q); +% Equation 34: lambda_3 + lambdad(i(j))=1/3+16*p/9/pi*(2*p^2-1)*Ek-(32*p^4-20*p^2+3)/9/pi/p*Kk; +% Equation 34: eta_1 + etad(i(j))=1/2/pi*(kap1(i(j))+p^2*(p^2+2*z(i(j))).*kap0(i(j))-(1+5*p^2+z(i(j)).^2)./4.*sqrt((1-x1(i(j))).*(x2(i(j))-1))); + if p==0.5 +% Case VIII: p=1/2, z=1/2 + lambdad(i(j))=1/3-4/pi/9; + etad(i(j))=3/32; + end; + end; + j=~j; + if any(j) +% Table 3, Case VI.: + %lam=0.5*pi; + q=2*p; + Kk=ellkVec(q); + Ek=ellecVec(q); +% Equation 34: lambda_4 + lambdad(i(j))=1/3+2/9/pi*(4*(2*p^2-1)*Ek+(1-4*p^2)*Kk); +% Equation 34: eta_2 + etad(i(j))=p^2/2*(p^2+2*z(i(j)).^2); + end +end + + +% the occulting star partly occults the source and crosses the limb: +% Table 3, Case III: +i=find(z>=0.5+abs(p-0.5) & z<1+p | (p>0.5 & z>abs(1-p)*(1+10*eps(1)) & z10*eps); +if any(i) + %lam=0.5*pi; + q=sqrt((x2(i)-x1(i))./(1-x1(i))); + Kk=ellkVec(q); + Ek=ellecVec(q); + n=x2(i)./x1(i)-1; + Pk=Kk-n/3.*rj(0,1-q.^2,1,1+n,TolFactor); + % Equation 34, lambda_2: + lambdad(i)=2/9/pi./sqrt(1-x1(i)).*((1-5*z(i).^2+p^2+x3(i).^2).*Kk+(1-x1(i)).*(z(i).^2+7*p^2-4).*Ek-3*x3(i)./x1(i).*Pk); + j=z(i)=1 + i=z<=p-1; + if any(i) + %lambdad(i)=1; + %etad(i)=1; + %lambdae(i)=1; + muo1(i)=0; + mu0(i)=0; + end; +end + +if ToSmooth>0 + deg=2; + CorrectionRange=2e-4; + FitRange=5*CorrectionRange; +% i=find(z>=1-p-CorrectionRange & z<=1-p+CorrectionRange); % inner than II/III contact +% if ~isempty(i) +% if p<1 +% tmpz=1-p+[linspace(-FitRange,-CorrectionRange,deg*10).'; linspace(CorrectionRange,FitRange,deg*10).']; +% else +% tmpz=p-1+[linspace(-FitRange,-CorrectionRange,deg*10).'; linspace(CorrectionRange,FitRange,deg*10).']; +% end +% tmpz(tmpz<0)=[]; +% tmp_muo1=occultquadVec(p,tmpz,u1,u2,0,'ToSmooth',false); +% if ~isempty(tmpz) +% pol=IterativePolyfit(tmpz,tmp_muo1,deg,3); +% muo1(i)=polyval(pol,z(i)); +% end +% end; + + % inner than II/III contact + i=z>=1-p-CorrectionRange & z<=1-p; + if ~isempty(i) + %if p<1 + % tmpz=1-p+linspace(-FitRange,-CorrectionRange,deg*10).'; + %else + % tmpz=p-1+linspace(-FitRange,-CorrectionRange,deg*10).'; + %end + tmpz=abs(1-p)+linspace(-FitRange,-CorrectionRange,deg*10).'; + tmpz(tmpz<0)=[]; + tmp_muo1=occultquadVec(p,tmpz,u1,u2,0,'ToSmooth',false); + if ~isempty(tmpz) + pol=polyfit(tmpz,tmp_muo1,deg); + muo1(i)=polyval(pol,z(i)); + end + end; + + % outer than II/III contact + j=z>=1-p & z<=1-p+CorrectionRange; + if ~isempty(j) + %if p<1 + % tmpz=1-p+linspace(CorrectionRange,FitRange,deg*10).'; + %else + % tmpz=p-1+linspace(CorrectionRange,FitRange,deg*10).'; + %end + tmpz=abs(1-p)+linspace(CorrectionRange,FitRange,deg*10).'; + tmpz(tmpz<0)=[]; + tmp_muo1=occultquadVec(p,tmpz,u1,u2,0,'ToSmooth',false); + if ~isempty(tmpz) + pol=polyfit(tmpz,tmp_muo1,deg); + muo1(j)=polyval(pol,z(j)); + end + end; + + % close to z=p + if p>0.5 + deg=3; + end; + CorrectionRange=1e-2; + FitRange=5*CorrectionRange; + if p<0.5 + i=find(z-p>-CorrectionRange & z-pmax(-CorrectionRange,1-2*p) & z-p0; +if any(i) + xt(i)=x(i); + yt(i)=y(i); + w(i)=1; +end; +i=~i; +if any(i) + xt(i)=x(i)-y(i); + yt(i)=-y(i); + w(i)=sqrt(x(i))./sqrt(xt(i)); +end +i=true(size(x)); +while any(i) + alamb(i)=2*sqrt(xt(i)).*sqrt(yt(i))+yt(i); + xt(i)=.25*(xt(i)+alamb(i)); + yt(i)=.25*(yt(i)+alamb(i)); + ave(i)=THIRD*(xt(i)+yt(i)+yt(i)); + s(i)=(yt(i)-ave(i))./ave(i); + i(i)=abs(s(i))>ERRTOL; +end; +rc=w.*(1+s.^2.*(C1+s.*(C2+s.*(C3+s.*C4))))./sqrt(ave); +% (C) Copr. 1986-92 Numerical Recipes Software + +function rj=rj(x,y,z,p,ERRTOL) +% built for x,z scalars and y,p vectors +fac=ones(size(y)); +x=x*fac; z=z*fac; rj=0*fac; +xt=rj; yt=rj; zt=rj; pt=rj; a=rj; b=rj; rho=rj; tau=rj; rcx=rj; +sqrtx=rj; sqrty=rj; sqrtz=rj; alamb=rj; alpha=rj; beta=rj; Sum=rj; ave=rj; +ave4=ave*ones(1,4); + +ERRTOL=.05*ERRTOL; +C1=3/14; +C2=1/3; +C3=3/22; +C4=3/26; +C5=.75*C3; +C6=1.5*C4; +C7=.5*C2; +C8=C3+C3; +i=p>0; +if any(i) + xt(i)=x(i); + yt(i)=y(i); + zt(i)=z(i); + pt(i)=p(i); +end; +i=~i; +if any(i) + xt(i)=min([x(i),y(i),z(i)],[],2); + zt(i)=max([x(i),y(i),z(i)],[],2); + yt(i)=x(i)+y(i)+z(i)-xt(i)-zt(i); + a(i)=1./(yt(i)-p(i)); + b(i)=a(i).*(zt(i)-yt(i)).*(yt(i)-xt(i)); + pt(i)=yt(i)+b(i); + rho(i)=xt(i).*zt(i)./yt(i); + tau(i)=p(i).*pt(i)./yt(i); + rcx(i)=rc(rho(i),tau(i),ERRTOL); +end +%continue % ************** 1 *************** +i=true(size(y)); +while any(i) + sqrtx(i)=sqrt(xt(i)); + sqrty(i)=sqrt(yt(i)); + sqrtz(i)=sqrt(zt(i)); + alamb(i)=sqrtx(i).*(sqrty(i)+sqrtz(i))+sqrty(i).*sqrtz(i); + alpha(i)=(pt(i).*(sqrtx(i)+sqrty(i)+sqrtz(i))+sqrtx(i).*sqrty(i).*sqrtz(i)).^2; + beta(i)=pt(i).*(pt(i)+alamb(i)).^2; + Sum(i)=Sum(i)+fac(i).*rc(alpha(i),beta(i),ERRTOL); + fac(i)=.25*fac(i); + xt(i)=.25*(xt(i)+alamb(i)); + yt(i)=.25*(yt(i)+alamb(i)); + zt(i)=.25*(zt(i)+alamb(i)); + pt(i)=.25*(pt(i)+alamb(i)); + ave(i)=.2*(xt(i)+yt(i)+zt(i)+pt(i)+pt(i)); + ave4(i,:)=ave(i)*ones(1,4); + %delx=(ave-xt)/ave; + %dely=(ave-yt)/ave; + %delz=(ave-zt)/ave; + %delp=(ave-pt)/ave; + del4=(ave4-[xt yt zt pt])./ave4; + i=max(abs(del4),[],2)>ERRTOL; +end; +ea=del4(:,1).*(del4(:,2)+del4(:,3))+del4(:,2).*del4(:,3); +eb=del4(:,1).*del4(:,2).*del4(:,3); +ec=del4(:,4).^2; +ed=ea-3*ec; +ee=eb+2*del4(:,4).*(ea-ec); +rj=3*Sum+fac.*(1+ed.*(-C1+C5*ed-C6*ee)+eb.*(C7+del4(:,4).*(-C8+del4(:,4).*C4))+del4(:,4).*ea.*(C2-del4(:,4)*C3)-C2*del4(:,4).*ec)./(ave.*sqrt(ave)); +i=p<=0; +if any(i) + rj(i)=a.*(b.*rj(i)+3*(rcx(i)-rf(xt(i),yt(i),zt(i),ERRTOL))); +end; +% (C) Copr. 1986-92 Numerical Recipes Software + +function rf=rf(x,y,z,ERRTOL) +xyz=[x y z]; +sqrtxyz=zeros(size(xyz)); +alamb=sqrtxyz; +ave=sqrtxyz; +ave3=sqrtxyz; +del3=sqrtxyz; + +ERRTOL=.08*ERRTOL; +THIRD=1/3; +C1=1/24; +C2=.1; +C3=3/44; +C4=1/14; +xyzt=xyz; +%continue % **************** 1 ***************** +i=true(size(x)); +while any(i) + sqrtxyz(i,:)=sqrt(xyz(i,:)); + %alamb(i,:)=sqrtxyz(i,1).*(sqrtxyz(i,2)+sqrtxy(i,3))+sqrtxyz(i,2).*sqrtxyz(i,3); + alamb(i,:)=sqrtxyz(i,1).*(sqrtxyz(i,2)+sqrtxyz(i,3))+sqrtxyz(i,2).*sqrtxyz(i,3); + xyzt(i,:)=.25*(xyzt(i,:)+alamb(i,:)*ones(1,3)); + ave(i)=THIRD*(sum(xyzt(i,:),2)); + ave3(i,:)=ave(i)*ones(1,3); + del3(i,:)=(ave3(i,:)-xyzt(i,:))./ave3(i,:); + i(i)=max(abs(del3(i,:)),[],2)>ERRTOL; +end; +e2=del3(:,1).*del3(:,2)-del3(:,3)^2; +e3=del3(:,1).*del3(:,2).*del3(:,3); +rf=(1+(C1*e2-C2-C3*e3).*e2+C4*e3)./sqrt(ave); +% (C) Copr. 1986-92 Numerical Recipes Software 0NL&WR2. diff --git a/LaplaceCoeffs.mat b/LaplaceCoeffs.mat new file mode 100644 index 0000000..a276742 Binary files /dev/null and b/LaplaceCoeffs.mat differ diff --git a/anac_2_2p.py b/anac_2_2p.py new file mode 100644 index 0000000..cf35f30 --- /dev/null +++ b/anac_2_2p.py @@ -0,0 +1,287 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +abgewandelt von: +Created on Fri Dec 3 08:49:07 2021 + +@author: dreizler +""" +import pdb + +from scipy.stats import gamma,norm,beta,truncnorm +import numpy as np + +from types import SimpleNamespace + +import dynesty +from dynesty.utils import resample_equal + +# system functions that are always useful to have +import time, sys, os + +# basic numeric setup +import numpy as np + +from dynesty import plotting as dyplot + +from oct2py import Oct2Py + +from multiprocessing import Pool + +# plotting +import matplotlib +from matplotlib import pyplot as plt +# inline plotting +#%matplotlib inline + +import csv +import json +import pickle + +def parse_parameters(settings): + p = [] + fix = {} + for planet in settings["planets"]: + for key, value in planet.items(): + fix[key] = [(None, None)] * len(settings["planets"]) + + i = 0 + for planet in settings["planets"]: + for key, value in planet.items(): + if value["type"] != "fixed": + p += [(i, key, value["type"], value["values"])] + else: + fix[key][i] = value["values"] + i += 1 + + for key, value in settings["parameters"].items(): + if value["type"] != "fixed": + p += [(-1, key, value["type"], value["values"])] + else: + fix[key] = value["values"] + + return p, fix + +def parse_settings(): + settings = json.load(open("parameter.json")) + + p, fix = parse_parameters(settings) + + return p, fix, settings["nkern"], settings["nlivepoints"] + + +def generate_mockdata(): + np.random.seed(56101) + + P_true=[12.3456, 29.3] + Tmid0_true=[9.8765432, 15.4029] + ror_true=[0.021, 0.027] + aor_true=30 + mu_true=[1.24e-5, 3e-5] + ex0_true=[0.01234, 0.00998] + ey0_true=[0.02345, 0.0189] + I0_true=[0.01745329, 0.0203] + Omega0_true=[0.4014257, 0.369] + sig_true = np.log(0.000084) + + #t_mock = np.linspace(0,100,500) + u1=0.5; + u2=0.3; + t=np.linspace(0,100,2) #0:0.002:100; + tRV=np.linspace(0,100,1000) + + init_octpy() + LC,RV_o_r,Y_o_r,Z_o_r,O = octave.feval("AnalyticLC",P_true, Tmid0_true, ror_true, aor_true, mu_true, ex0_true, ey0_true, I0_true, Omega0_true, t, u1, u2,'tRV',tRV, nout=5) + + #y_old = octave.AnalyticLC(P_true, Tmid0_true, ror_true, aor_true, mu_true, ex0_true, ey0_true, I0_true, Omega0_true, t, u1, u2) + y_true2 = [] + for i in RV_o_r: + y_true2.append(i[0]) + y_true2 = np.array(y_true2) + y_true = y_true2[(np.isfinite(y_true2))] + + y = y_true + 0.0001*np.random.randn(len(y_true)) #0.0001 + yerr = 0.00005 + np.zeros(len(y)) #0.00005 + + # plot results + plt.figure(figsize=(10, 5)) + plt.errorbar(tRV, y, yerr=yerr, fmt='ko', ecolor='red') + plt.plot(tRV, y_true, color='blue', lw=3) + plt.xlabel(r'$t$') + plt.ylabel(r'$y_true$') + plt.tight_layout() + plt.show() + plt.close() + + return y, yerr, y_true + +def calculate(nkern, nlivepoints, parameters, fixed_parameters): + ndim = len(parameters) + p = Pool(nkern, init, [parameters, fixed_parameters]) + + #P, Tmid0, ror, aor, mu, ex0, ey0, I0, Omega0, sig + dsampler = dynesty.NestedSampler(loglike, prior_transform, ndim=ndim, periodic=None, + bound='multi', sample='auto',nlive=nlivepoints,pool=p, queue_size=nkern) #n*n/2 live points mind. + dsampler.run_nested(dlogz=0.001) + dres = dsampler.results + with open('data_4dim.pickle', 'wb') as f: + pickle.dump(dres, f, pickle.HIGHEST_PROTOCOL) + + return dres + +def create_plot(dres): + P_true=[12.3456, 29.3] + Tmid0_true=[9.8765432, 15.4029] + ror_true=[0.021, 0.027] + aor_true=30 + mu_true=[1.24e-5, 3e-5] + ex0_true=[0.01234, 0.00998] + ey0_true=[0.02345, 0.0189] + I0_true=[0.01745329, 0.0203] + Omega0_true=[0.4014257, 0.369] + sig_true = np.log(0.000084) + + truths = [P_true[0], P_true[1], sig_true] + labels = [r'P $d$',r'P $d$', r'T $d$',r'T $d$', r'ex0', r'ex0', r'ey0', r'ey0',r'$\sigma \log(m/s)$'] + fig, axes = dyplot.traceplot(dres, truths=truths, truth_color='black', labels=labels, connect=True, + fig=plt.subplots(ndim, 2, constrained_layout=True)) + fig.tight_layout() + plt.savefig("tmp_2p/plot_test1.png") + fig, axes = dyplot.cornerplot(dres, truths=truths, show_titles=True, + title_kwargs={'y': 1.04}, labels=labels, + fig=plt.subplots(ndim, ndim,constrained_layout=True)) + + fig.tight_layout() + plt.savefig("tmp_2p/plot_test2.png") + #lt.close() + +def init(parameters, fixed_parameters): + global p + p = parameters + global fix + fix = fixed_parameters + init_octpy() + +def init_octpy(): + global octave + octave = Oct2Py() + BaseDir='/home/end/anal/endsversion/' + octave.addpath(BaseDir) + octave.addpath(BaseDir + 'AnalyticLC/') + octave.addpath(BaseDir + 'FittingAndMath/') + octave.addpath(BaseDir + 'FluxCalculation/') + octave.LaplaceCoeffs('load',BaseDir + 'LaplaceCoeffs.mat') + +# Prior transforms for nested samplers: +def transform_uniform(x, hyperparameters): + a, b = hyperparameters + return a + (b-a)*x + +def transform_loguniform(x, hyperparameters): + a, b = hyperparameters + la = np.log(a) + lb = np.log(b) + return np.exp(la + x * (lb - la)) + +def transform_normal(x, hyperparameters): + mu, sigma = hyperparameters + return norm.ppf(x, loc=mu, scale=sigma) + +def transform_beta(x, hyperparameters): + a, b = hyperparameters + return beta.ppf(x, a, b) + +def transform_exponential(x, hyperparameters): + a = hyperparameters + return gamma.ppf(x, a) + +def transform_truncated_normal(x, hyperparameters): + mu, sigma, a, b = hyperparameters + ar, br = (a - mu) / sigma, (b - mu) / sigma + return truncnorm.ppf(x, ar, br, loc=mu, scale=sigma) + +def transform_modifiedjeffreys(x, hyperparameters): + turn, hi = hyperparameters + return turn * (np.exp( (x + 1e-10) * np.log(hi/turn + 1)) - 1) + +# prior transform +def prior_transform(utheta): + global p; + global fix; + + transforms = [x[2] for x in p] + hyperparams = [x[3] for x in p] + transformed_priors = np.zeros(len(utheta)) + for k,param in enumerate(utheta): + if transforms[k] == 'uniform': + transformed_priors[k] = transform_uniform(param,hyperparams[k]) + elif transforms[k] == 'loguniform': + transformed_priors[k] = transform_loguniform(param,hyperparams[k]) + elif transforms[k] == 'normal': + transformed_priors[k] = transform_normal(param,hyperparams[k]) + else: + raise ValueError("Transform not known") + + return transformed_priors + +# log-likelihood +def loglike(theta): + global octave; + global p; + global fix; + parameters = fix.copy() + + for i in range(len(theta)): + planetid = p[i][0] + keyname = p[i][1] + if planetid == -1: + parameters[keyname] = theta[i] + else: + parameters[keyname][planetid] = theta[i] + + n = SimpleNamespace(**parameters) + + u1=0.5; + u2=0.3; + t=np.linspace(0,100,2) #0:0.002:100; + tRV=np.linspace(0,100,1000) + + try: + LC,RV_o_r,Y_o_r,Z_o_r,O = octave.feval("AnalyticLC",n.P, n.Tmid0, n.ror, n.aor, n.mu, n.ex0, n.ey0, n.I0, n.Omega0, t, u1, u2,'tRV',tRV, nout=5) + + model =[] + for i in RV_o_r: + model.append(i[0]) + model = np.array(model) + model23 = model[~(np.isnan(model))] + model2 = model23[~(np.iscomplex(model23))] + #print("Model:") + #print(model2) + #print("y:") + #print(y) + if (model2.size == 0): + #print("error") + return -np.inf + inv_sigma2 = 1.0 / (n.yerr**2 + np.exp(2 * n.sig)) + ret = -0.5 * (np.sum((n.y-model2)**2 * inv_sigma2 - np.log(inv_sigma2))) + #if np.imag(ret)!=0: + # return -np.inf + #print(ret) + return ret + except Exception as e: + print(e) + return -np.inf + +def main(): + parameters, fixed, nkern, nlivepoints = parse_settings() + + y, yerr, y_true = generate_mockdata() + + fixed["y"] = y + fixed["yerr"] = yerr + fixed["y_true"] = y_true + + results = calculate(nkern, nlivepoints, parameters, fixed) + create_plot(results) + +main() diff --git a/octave-workspace b/octave-workspace new file mode 100644 index 0000000..2db5e86 Binary files /dev/null and b/octave-workspace differ diff --git a/parameter.json b/parameter.json new file mode 100644 index 0000000..ad93ebb --- /dev/null +++ b/parameter.json @@ -0,0 +1,39 @@ +{ + "planets": [ + { + "P": {"type": "uniform", "values": [11,13]}, + "Tmid0": {"type": "fixed", "values": 9.876}, + "ex0": {"type": "fixed", "values": 0.01234}, + "ey0": {"type": "fixed", "values": 0.02345}, + "ror": {"type": "fixed", "values": 0.021}, + "mu": {"type": "fixed", "values": 1.24e-5}, + "I0": {"type": "fixed", "values": 0.01745}, + "Omega0": {"type": "fixed", "values": 0.4014257} + }, + { + "P": {"type": "uniform", "values": [29,31]}, + "Tmid0": {"type": "fixed", "values": 15.4029}, + "ex0": {"type": "fixed", "values": 0.00998}, + "ey0": {"type": "fixed", "values": 0.0189}, + "ror": {"type": "fixed", "values": 0.027}, + "mu": {"type": "fixed", "values": 3e-5}, + "I0": {"type": "fixed", "values": 0.0203}, + "Omega0": {"type": "fixed", "values": 0.369} + } + ], + "parameters": { + "sig": {"type": "uniform", "values": [-10000,100]}, + "aor": {"type": "fixed", "values": 30}, + "u1" : {"type": "fixed", "values": 0.5}, + "u2" : {"type": "fixed", "values": 0.3} + }, + + "LC_datapath": "/dev/null", + "RV_datapath": "/dev/null", + "AnalyticLC_datapath": "/dev/null", + "atrometry_datapath": "/dev/null", + + "nlivepoints": 1000, + "nkern": 16 +} + diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..36c8342 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,3 @@ +oct2py +dynesty +matplotlib