-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathMission.m
162 lines (132 loc) · 7.64 KB
/
Mission.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
153
154
155
156
157
158
159
160
161
%# Mission is a class for determining Interplanetary Trajectories based on given Bodies and given time inputs
classdef Mission
%#
%# PROPERTIES:
%#
%#
%# Trajectory; : Trajectory Of Mission
%# Mission_Times; : Array of relative Encounter Times (first element is Absolute Time)
%# Absolute_Times; : Array of Absolute Encounter Times
%# Spice_Min_Times; : Array of Minimum Times for SPICE Calculation
%# Spice_Max_Times; : Array of Maximum Times for SPICE Calculation
%# Out_Of_Spice_Bounds =0; : Flag to indicate SPICE was unable to be used for a body i
%# TotaldV; : Total DeltaV of Mission
%# FlybyRendez=0; : Flyby or Rendezvous Flag =0 (FLYBY) = 1 (RENDEZVOUS)
%# wayflag=1; : Default to Prograde Only
%#
%# METHODS:
%#
%#
%# Mission : Class Constructor
%# Set_Mission_Times : Function to Calculate Mission Times from input array of Absolute Times
%# Set_Absolute_Times : Function to Calculate Absolute Times from input array of Mission Times
%# Compute_DeltaV : Function to Update Trajectory Based on input Mission Times and output DeltaV
%#
properties
Trajectory; % Trajectory Of Mission
Mission_Times; % Array of relative Encounter Times
Absolute_Times; % Array of Absolute Encounter Times
Spice_Min_Times; % Array of Minimum Times for SPICE Calculation
Spice_Max_Times; % Array of Maximum Times for SPICE Calculation
Out_Of_Spice_Bounds =0; % Flag to indicate SPICE was unable to be used for a body i
TotaldV; % Total DeltaV of Mission
FlybyRendez=0; % Flyby or Rendezvous Flag =0 (FLYBY) = 1 (RENDEZVOUS)
wayflag=1; % Default to Prograde Only
home_periapsis=0; % Periapsis radial distance from Home's Centre
target_periapsis=0; % Periapsis radial distance from Target's Centre
end
methods
function obj = Mission( bodies , times , MinSpice, MaxSpice)
%# Mission : Class Constructor
% Initialise SPICE Min & MAx Times
obj.Spice_Min_Times = MinSpice
obj.Spice_Max_Times = MaxSpice
% Initialise Positions and Osculating Orbits of Various Bodies
for i=1:size(bodies,2)
if (contains(bodies(i).ID,'CUSTOM BODY','IgnoreCase',true))
mode(i)=0;
bodies(i)=bodies(i).compute_ephem_at_t(times(i),0,1e-4);
bodies(i)=bodies(i).calculate_orbit_from_ephem(times(i));
continue;
end
if times(i)<obj.Spice_Min_Times(i)
mode(i)=1;
obj.Out_Of_Spice_Bounds=bitor(obj.Out_Of_Spice_Bounds,2^i,'uint32');
bodies(i)=bodies(i).compute_ephem_at_t(obj.Spice_Min_Times(i),2,1e-4);
bodies(i)=bodies(i).calculate_orbit_from_ephem(obj.Spice_Min_Times(i));
bodies(i).ephem0=bodies(i).ephemt;
elseif times(i)>obj.Spice_Max_Times(i)
mode(i)=1;
obj.Out_Of_Spice_Bounds=bitor(obj.Out_Of_Spice_Bounds,2^i,'uint32');
bodies(i)=bodies(i).compute_ephem_at_t(obj.Spice_Max_Times(i),2,1e-4);
bodies(i)=bodies(i).calculate_orbit_from_ephem(obj.Spice_Max_Times(i));
bodies(i).ephem0=bodies(i).ephemt;
else
if bitand(obj.Out_Of_Spice_Bounds,2^i)
obj.Out_Of_Spice_Bounds=obj.Out_Of_Spice_Bounds-2^i;
end
mode(i)=2;
bodies(i)=bodies(i).compute_ephem_at_t(times(i),mode(i),1e-4);
bodies(i)=bodies(i).calculate_orbit_from_ephem(times(i));
end
end
% Create an object of Type Nbody_Trajectory_With_Encounters
obj.Trajectory = Nbody_Trajectory_With_Encounters(bodies);
obj.Mission_Times = zeros(1,obj.Trajectory.Nbody);
obj.Absolute_Times = zeros(1,obj.Trajectory.Nbody);
obj = obj.Set_Mission_Times( times );
[obj DeltaV gradient] = obj.Compute_DeltaV( obj.Mission_Times );
end
function obj = Set_Mission_Times( obj, times )
%# Set_Mission_Times : Function to Calculate Mission Times from input array of Absolute Times
% Function to Initialise Mission_Times
% Input times is array of absolute times
% Output has first item as launch time and
% remaining times as cruise durations
obj.Mission_Times(1) = times(1);
for i=2:obj.Trajectory.Nbody
obj.Mission_Times(i)=times(i)-times(i-1);
end
end
function obj = Set_Absolute_Times( obj, times )
%# Set_Absolute_Times : Function to Calculate Absolute Times from input array of Mission Times
% Reverse of Set_Mission_Times
obj.Absolute_Times(1)=times(1);
for i=2:obj.Trajectory.Nbody
obj.Absolute_Times(i)=times(i)+obj.Absolute_Times(i-1);
end
end
function [obj, DeltaV, gradient] = Compute_DeltaV( obj, times )
%# Compute_DeltaV : Function to Update Trajectory Based on input Mission Times and output DeltaV
obj = obj.Set_Absolute_Times(times);
for i=1:obj.Trajectory.Nbody
if (obj.Trajectory.Body_Set(i).Fixed_Point==-1)&&(contains(obj.Trajectory.Body_Set(i).ID,'CUSTOM BODY','IgnoreCase',true))
mode(i)=0;
obj.Trajectory.Body_Set(i)=obj.Trajectory.Body_Set(i).compute_ephem_at_theta(obj.Trajectory.Body_Set(i).true_anomaly);
% obj.Trajectory.Body_Set(i)=obj.Trajectory.Body_Set(i).calculate_orbit_from_ephem(obj.Absolute_Times(i));
continue;
end
if obj.Absolute_Times(i)<obj.Spice_Min_Times(i)
mode(i)=1;
obj.Out_Of_Spice_Bounds=bitor(obj.Out_Of_Spice_Bounds,2^i,'uint32');
obj.Trajectory.Body_Set(i)=obj.Trajectory.Body_Set(i).compute_ephem_at_t(obj.Spice_Min_Times(i),2,1e-4);
obj.Trajectory.Body_Set(i)=obj.Trajectory.Body_Set(i).calculate_orbit_from_ephem(obj.Spice_Min_Times(i));
obj.Trajectory.Body_Set(i).ephem0=obj.Trajectory.Body_Set(i).ephemt;
elseif obj.Absolute_Times(i)>obj.Spice_Max_Times(i)
mode(i)=1;
obj.Out_Of_Spice_Bounds=bitor(obj.Out_Of_Spice_Bounds,2^i,'uint32');
obj.Trajectory.Body_Set(i)=obj.Trajectory.Body_Set(i).compute_ephem_at_t(obj.Spice_Max_Times(i),2,1e-4);
obj.Trajectory.Body_Set(i)=obj.Trajectory.Body_Set(i).calculate_orbit_from_ephem(obj.Spice_Max_Times(i));
obj.Trajectory.Body_Set(i).ephem0=obj.Trajectory.Body_Set(i).ephemt;
else
obj.Out_Of_Spice_Bounds=bitset(obj.Out_Of_Spice_Bounds,i+1,0,'uint32');
mode(i)=2;
end
end
obj.Trajectory = obj.Trajectory.Compute_Total_Deltav( obj.Absolute_Times, mode, 1e-4, 1000,obj.wayflag,obj.FlybyRendez,obj.home_periapsis,obj.target_periapsis);
DeltaV = obj.Trajectory.BestDeltaV;
gradient=[];
obj.TotaldV = DeltaV;
end
end
end