-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathLyapOrbitParameters.m
153 lines (119 loc) · 5.66 KB
/
LyapOrbitParameters.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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
%{
...
Created on 9/3/2020 18:47
This File does the continuation and gets all the Lyapunov orbit parameters.
Inputs
------
1) UserDat - Supplied UserData in main file.
2) G_var - Global data.
Outputs
--------
1) LyapOrb - A structure (1 x NoofEqPoints) containing
- time - full time period for each lyapunov orbits.
- IC - Initial Conditions for all Lyapunov Orbits.
- Energy - Energy Value for all Lyapunov Orbits.
- Monodromy - Monodromy matrix of all Lyapunov Orbits.
- Eigens - Eigen Values and Eigen Vectors for all Lyapunov Orbits.
- Stable,Unstable and Center Eigen Values
- Stable,Unstable and Center Eigen Vectors
Note that the size of eigen values and eigen vectors might change
Dependencies
------------
1) InitialGuess(PointLoc,G_var)
2) DiffCorrec(X_Guess,Plot,G_var)
Reference for continuation
--------------------------
1) Wang Sang Koon, Martin W Lo,JE Marsden, Shane D Ross - "Dynamical
Systems,the Three Body Problem and Space Mission Design", 2011
...
%}
function [LyapOrb] = LyapOrbitParameters(UserDat,G_var)
% Extract the parameters
NoofFam = UserDat.NoofFam;
CorrecPlot = UserDat.CorrectionPlot;
NoofEqPoints = UserDat.PointLoc;
Dimension = UserDat.Dimension;
mu = G_var.Constants.mu;
funVarEq = G_var.IntFunc.VarEqAndSTMdot;
% Initialize
tCorrec = zeros(NoofFam,1);
xCorrec = zeros(NoofFam,2*Dimension);
Energy = zeros(NoofFam,1);
Monodromy = zeros(2*Dimension,2*Dimension,NoofFam);
for PointLoc = 1:NoofEqPoints
fprintf('*******************************************\n')
fprintf('Starting the process for L%d Eq.Point\n',PointLoc)
fprintf('*******************************************\n')
[XGuessL] = InitialGuess(PointLoc,G_var);
if Dimension == 3
XGuessL = XGuessL(1);
elseif Dimension == 2
XGuessL = XGuessL(2);
end
family = 1;
fprintf('===============================================\n')
fprintf('Obtaining the Corrected Values for family:- %d\n',family)
fprintf('===============================================\n')
fprintf('\nThe Guess X value is : [%f %f %f %f %f %f]\n',XGuessL.one)
[tCorrec(family,1),xCorrec(family,:),DF(:,:,family)] = DiffCorrec(XGuessL.one,CorrecPlot,G_var);
fprintf('\n')
fprintf('Calculating the orbits Parameters ..\n')
fprintf('\n')
Energy(family,1) = jacobiValue3D(xCorrec(family,:),mu);
[~,Monodromy(:,:,family),~,~] = StateTransAndX(G_var,xCorrec(family,:),funVarEq,tCorrec(family,1));
[Eigens(family).S_EigVal,Eigens(family).US_EigVal,Eigens(family).C_Val,Eigens(family).S_EigVec,...
Eigens(family).US_EigVec,Eigens(family).C_EigVec] = CalcEigenValVec(Monodromy(:,:,family),1) ;
family = 2;
fprintf('\n===============================================\n')
fprintf('Obtaining the Corrected Values for family:- %d\n',family)
fprintf('===============================================\n')
fprintf('\nThe Guess X value is : [%f %f %f %f %f %f]\n',XGuessL.two)
[tCorrec(family,1),xCorrec(family,:),DF(:,:,family)] = DiffCorrec(XGuessL.two,CorrecPlot,G_var);
fprintf('\n')
fprintf('Calculating the orbits Parameters ..\n')
fprintf('\n')
Energy(family,1) = jacobiValue3D(xCorrec(family,:),mu);
[~,Monodromy(:,:,family),~,~] = StateTransAndX(G_var,xCorrec(family,:),funVarEq,tCorrec(family,1));
[Eigens(family).S_EigVal,Eigens(family).US_EigVal,Eigens(family).C_Val,Eigens(family).S_EigVec,...
Eigens(family).US_EigVec,Eigens(family).C_EigVec] = CalcEigenValVec(Monodromy(:,:,family),1) ;
% fprintf('Energy of the orbit:- %d\n',Energy(family,1))
% fprintf('Monodromy Matrix:- %d\n',Monodromy(:,:,family))
% fprintf('Stable Eigen Value:- %d\n',Eigens(family).S_EigVal)
% fprintf('Unstable Eigen Value:- %d\n',Eigens(family).US_EigVal)
% fprintf('Center Eigen Value:- %d\n',Eigens(family).C_Val)
% fprintf('Stable Eigen Vector:- %d\n',Eigens(family).S_EigVec)
% fprintf('Unstable Eigen Vector:- %d\n',Eigens(family).US_EigVec)
% fprintf('Center Eigen Vector:- %d\n',Eigens(family).C_EigVec)
%% Start Continuation
for family = 3:NoofFam
fprintf('\n===============================================\n')
fprintf('Obtaining the Corrected Values for family:- %d\n',family)
fprintf('===============================================\n')
delta = xCorrec(family-1,:) - xCorrec(family-2,:);
GuessX = xCorrec(family-1,:) + delta;
fprintf('\nThe Guess X value is : [%f %f %f %f %f %f]\n',GuessX)
[tCorrec(family,1),xCorrec(family,:),DF(:,:,family)] = DiffCorrec(GuessX,CorrecPlot,G_var);
fprintf('\n')
fprintf('Calculating the orbits Parameters ..\n')
fprintf('\n')
Energy(family,1) = jacobiValue3D(xCorrec(family,:),mu);
[~,Monodromy(:,:,family),~,~] = StateTransAndX(G_var,xCorrec(family,:),funVarEq,tCorrec(family,1));
[Eigens(family).S_EigVal,Eigens(family).US_EigVal,Eigens(family).C_Val,Eigens(family).S_EigVec,...
Eigens(family).US_EigVec,Eigens(family).C_EigVec] = CalcEigenValVec(Monodromy(:,:,family),1) ;
% fprintf('Energy of the orbit:- %d\n',Energy(family,1))
% fprintf('Monodromy Matrix:- %d\n',Monodromy(:,:,family))
% fprintf('Stable Eigen Value:- %d\n',Eigens(family).S_EigVal)
% fprintf('Unstable Eigen Value:- %d\n',Eigens(family).US_EigVal)
% fprintf('Center Eigen Value:- %d\n',Eigens(family).C_Val)
% fprintf('Stable Eigen Vector:- %d\n',Eigens(family).S_EigVec)
% fprintf('Unstable Eigen Vector:- %d\n',Eigens(family).US_EigVec)
% fprintf('Center Eigen Vector:- %d\n',Eigens(family).C_EigVec)
%% Store all the values
LyapOrb(PointLoc).time = tCorrec; %(NoofFam x 1) - Full Orbit Time
LyapOrb(PointLoc).IC = xCorrec; %(NoofFam x UserDat.Dimension)
LyapOrb(PointLoc).Energy = Energy;
LyapOrb(PointLoc).Monodromy = Monodromy;
LyapOrb(PointLoc).Eigens = Eigens;
end
end
end