-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathDfmatrix3D.m
127 lines (93 loc) · 3.15 KB
/
Dfmatrix3D.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
%{
...
Gives matrix Df(x) (i.e., the matrix of derivatives of f, where xdot=f(x) )
for a 4-dimensional point x in the planar CR3BP's phase space
-----------------------------------------------------------------------
CR3BP (Circular Restricted Three-Body [Gravitational] Problem)
with the LARGER MASS, M1 to the left of the origin at (-mu,0)
and the smaller mass, M2, or the planet (ie. Earth), is at (1 - mu, 0)
(rotating coords)
L4
-L3------M1--+-----L1--M2--L2-
L5
The following is the Jacobian matrix
Df =[ 0 0 1 0 ;
0 0 0 1 ;
-Uxx -Uxy 0 2 ;
-Uxy -Uyy -2 0 ];
====================================================
Source File : Prof.Shane Ross (revised 2.19.04)
WebSite(Codes) : http://www.dept.aoe.vt.edu/~sdross/
Modification
1) Added a switch-case to calculate for planar and 3-D verrsions on
Feb 20,2020 by 15:30
====================================================
Inputs:
-------
1) mu - Mass parameter
2) x - X vector
Outputs
-------
1) Df - Jacobian Matrix
...
%}
function Df = Dfmatrix3D(x,mu)
if length(x)>4
type = 'ThreeDim';
else
type = 'Planar';
end
% mu = mass paramater
mu1 = 1-mu ;
mu2 = mu ;
switch type
case 'Planar'
r2= (x(1)+mu2)^2 + x(2)^2; % r: distance to m1, LARGER MASS
R2= (x(1)-mu1)^2 + x(2)^2; % R: distance to m2, smaller mass
r3= r2^1.5;
r5= r2^2.5;
R3= R2^1.5;
R5= R2^2.5;
%The following are three double partial derivatives of the
% effective potential U(x,y)
% First aprtial currently unused
Ux = x(1) - mu1*(x(1)+mu2)/r3 - mu2*(x(1)-mu1)/R3 ;
Uy = x(2) - mu1* x(2) /r3 - mu2* x(2) /R3 ;
Uxx = 1 - mu1/r3 + (3*mu1*(x(1)+mu)^2)/r5 - mu2/R3 + (3*mu*(x(1)-1+mu)^2)/R5 ;
Uyy = 1 - mu1/r3 + (3*mu1*x(2)^2)/r5 - mu2/R3 + (3*mu2*x(2)^2)/R5 ;
Uxy = (3*mu1*(x(1)+mu2)*x(2))/r5 + (3*mu2*(x(1)-1+mu2)*x(2))/R5 ;
Uyx = Uxy;
I = eye(2);
zer = zeros(2,2);
UXX = [Uxx Uxy; Uyx Uyy];
sig = [0 2;-2 0];
Df = [zer I;UXX sig];
%==================================================================================
case 'ThreeDim'
r2= (x(1)+mu2)^2 + x(2)^2 + x(3)^2; % r: distance to m1, LARGER MASS
R2= (x(1)-mu1)^2 + x(2)^2 + x(3)^2; % R: distance to m2, smaller mass
r3= r2^1.5;
r5= r2^2.5;
R3= R2^1.5;
R5= R2^2.5;
% The following are the two partial derivatives of the
% effective potential U(x,y)
% This is first partial currently unused
Ux = x(1) - mu1*(x(1)+mu2)/r3 - mu2*(x(1)-mu1)/R3 ;
Uy = x(2) - mu1* x(2) /r3 - mu2* x(2) /R3 ;
Uz = -(mu1*x(3))/r3 - (mu*x(3))/R3;
Uxx = 1 - mu1/r3 + (3*mu1*(x(1)+mu)^2)/r5 - mu2/R3 + (3*mu*(x(1)-1+mu)^2)/R5 ;
Uyy = 1 - mu1/r3 + (3*mu1*x(2)^2)/r5 - mu2/R3 + (3*mu2*x(2)^2)/R5 ;
Uzz = -mu1/r3 + (3*mu1*x(3)^2)/r5 - mu2/R3 + (3*mu2*x(3)^2)/R5 ;
Uxy = (3*mu1*(x(1)+mu2)*x(2))/r5 + (3*mu2*(x(1)-1+mu2)*x(2))/R5 ;
Uyx = Uxy;
Uxz = (3*mu1*(x(1)+mu2)*x(3))/r5 + (3*mu2*(x(1)-1+mu2)*x(3))/R5 ;
Uzx = Uxz ;
Uyz = (3*mu1*x(2)*x(3))/r5 + (3*mu2*x(2)*x(3))/R5 ;
Uzy = Uyz ;
I = eye(3);
zer = zeros(3,3);
UXX = [Uxx Uxy Uxz; Uyx Uyy Uyz;Uzx Uzy Uzz];
sig = [0 2 0;-2 0 0;0 0 0];
Df = [zer I;UXX sig];
end