1230 lines
68 KiB
Matlab
1230 lines
68 KiB
Matlab
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
|