1 % Started 12/16/25
2 % Truman Abbe | Utah State University | truman.abbe23@[Link]
3 % quaternion8.1.m
4
5
6 %%%%%%%%%%%%% RUN SIMULATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7
8 % initial state y0 is given with scalar first quaternion:
9 % [w, x, y, z, w_x, w_y, w_z]
10 y0 = [1 0 0 0 deg2rad(5) deg2rad(7) deg2rad(9)]';
11
12 fprintf('Starting 30-day simulation... this may take a few minutes.\n');
13 tic; % starts timer
14
15 % tolerances are relaxed to speed up integration
16 options = odeset('RelTol', 1e-4, 'AbsTol', 1e-4, 'MaxStep', 5);
17
18 % define 30 days in seconds
19 t_end = 30 * 24 * 3600;
20
21 % use a specified output time vector to save memory
22 % instead of saving every tiny step, save data every 120 seconds
23 t_span = 0:120:t_end;
24
25 % stiff solver ode15s is used to accommodate magnetic torques
26 [T, Y] = ode15s(@myODE, t_span, y0, options);
27
28
29 elapsedTime = toc; % ends timer
30 fprintf('Simulation complete in %.2f seconds.\n', elapsedTime);
31 final_state = Y(end, :);
32 disp(T(end, :)); % time elapsed is displayed for debugging
33 disp(final_state); % final state is displayed for debugging
34
35
36 %%%%%% FINAL MISSION DASHBOARD: 14/30-DAY ANALYTICS %%%%%%%%%%%%%%%%%%%%%%%%%%
37 T_days = T(:) / 86400; % column vector
38 Y = Y(:,:); % ensures matrix is in correct shape
39
40 % downsampling factor for visual clarity (plots every 10th point)
41 % this prevents "solid block" artifacts in long-duration simulations
42 idx = 1:20:size(Y,1);
43
44 %disp(size(T_days))
45 %disp(size(Y))
46 %disp(idx(1:10))
47
48 % --- Minimal test to check plotting ---
49 figure;
50 plot(T_days(idx), Y(idx,1), 'r', 'LineWidth', 1);
51 title('Minimal Test: First Quaternion');
52 xlabel('Time (Days)'); ylabel('q_w');
53 grid on;
54
55 fig = figure('Name', 'Satellite Mission Dashboard: Long-Term Analysis', 'NumberTitle', 'off',
'Color', 'w');
56
57 % --- 1. Top Left: Raw Quaternions (Cleaned) ---
58 subplot(2,2,1);
59 plot(T_days(idx), Y(idx, 1:4), 'LineWidth', 1);
60 title('State Vector: Raw Quaternions');
61 xlabel('Time (Days)'); ylabel('Component Value');
62 legend('q_w', 'q_x', 'q_y', 'q_z', 'Location', 'northeast');
63 grid on; xlim([0 T_days(end)]);
64 ylim([-1.1 1.1]);
65
66 % --- 2. Top Right: High-Precision Norm Error ---
67 subplot(2,2,2);
68 q_norm = sqrt(sum(Y(:,1:4).^2, 2));
69 error_norm = 1 - q_norm;
70 plot(T_days(idx), error_norm(idx), 'Color', [0, 0.447, 0.741]);
71 title('High-Precision Norm Error (1 - ||q||)');
72 xlabel('Time (Days)'); ylabel('Error Magnitude');
73 grid on; xlim([0 T_days(end)]);
74
75 % --- 3. Bottom Left: Angular Velocity (The Physics) ---
76 subplot(2,2,3);
77 Y_deg = Y(:, 5:7) * (180/pi);
78 w_total_deg = sqrt(sum(Y_deg.^2, 2));
79
80 plot(T_days, Y_deg, 'LineWidth', 1); hold on;
81 plot(T_days, w_total_deg, 'k', 'LineWidth', 2);
82
83 % Detumble Threshold Line
84 yline(1, '--r', 'Detumble Threshold', 'LabelVerticalAlignment', 'bottom');
85
86 title('Physics: Angular Velocity Decay');
87 xlabel('Time (Days)'); ylabel('deg/sec');
88 legend('\omega_x', '\omega_y', '\omega_z', '|\omega|_{total}');
89 grid on; xlim([0 T_days(end)]);
90 ylim([0 max(w_total_deg)*1.1]); % Ensures magnitude is always visible
91
92 % --- 4. Bottom Right: Unit Quaternion Health (Visual Fix) ---
93 subplot(2,2,4);
94 % We relax the ylim slightly to avoid the "Solid Red Block" effect
95 plot(T_days(idx), q_norm(idx), 'r', 'LineWidth', 1.5);
96 title('Unit Quaternion Health (Stability)');
97 xlabel('Time (Days)'); ylabel('Magnitude');
98 grid on; xlim([0 T_days(end)]);
99
100 % adjusting ylim to show a stable line rather than a thick band
101 % if error is 1e-4, this range shows the stability clearly.
102 ylim([0.999 1.001]);
103
104 % Auto-save the figure as a high-res PNG for your report
105 drawnow;
106 saveas(fig, 'Final_30Day_Mission_Dashboard8.[Link]');
107
108
109 %%%%%%%%%%%% ODE FUNCTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
110 function dydt = myODE(t, y)
111 % y is [qw, qx, qy, qz, wx, wy, wz] column vector
112
113 % calculates earth's magnetic flux density
114 pos = circularPropagator(t); % orbital position
115 B_eci = getEarthField(pos); % magnetic field components given position
116 R_bi = quat2rotm_manual(y(1:4)');
117 B_body = R_bi * B_eci; % conversion to body frame
118
119 % calculates magnetic torques of hysteresis rods
120 V_hyst = pi * (0.0005)^2 * (0.095) * 2; % rod volume
121 B_s = 0.027; H_c = 12; B_r = 0.004; % apparent values from paper
122 m_hyst = calculate_m_hyst(y, B_body, V_hyst, B_s, H_c, B_r);
123 tau_hyst = cross(m_hyst, B_body);
124
125 % permanent magnet torques
126 m_pm_val = 0.003; % dipole strength in Ampere-meters^2
127 % assume the magnet is aligned with the satellite's Z-axis (long axis)
128 m_pm_body = [0; 0; m_pm_val];
129 tau_pm = cross(m_pm_body, B_body);
130
131 %if norm(y(5:7)) < deg2rad(1)
132 % tau_total = tau_hyst + tau_pm;
133 %else
134 % tau_total = tau_hyst;
135 %end
136
137 % total external torques
138 tau_total = tau_hyst; % tau_pm;
139
140 % unpack quaternion and angular velocity into row vectors
141 q = y(1:4)';
142 w = y(5:7)';
143
144 % get pure quaternion row
145 wq = [0 w];
146
147 % compute qdot
148 qdot = 0.5 * quatmultiply(q, wq);
149
150 % compute wdot
151 Ix = 0.01; Iy = 0.01; Iz = 0.01;
152 I = diag([Ix Iy Iz]);
153 wdot = (I \ (tau_total' - cross(w', I*w')))';
154
155 % baumgarte stabilization
156 % this forces the norm back to 1.0 if it drifts
157 k = 0.2; % gain factor
158 norm_error = 1 - (q(1)^2 + q(2)^2 + q(3)^2 + q(4)^2);
159 q_dot_corrected = qdot + k * norm_error * q;
160
161 final_q_dot = q_dot_corrected(1:4);
162 final_w_dot = wdot(1:3);
163
164 % Pack into a single 7x1 COLUMN vector
165 dydt = [final_q_dot(:); final_w_dot(:)];
166 end
167
168
169 %%%%%%%%%%%%%%%%HAMILTON PRODUCT%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
170 function r = quatmultiply(p, q)
171 % r = quatmultiply(p, q)
172 % Performs the Hamilton product (p * q) for quaternions.
173 % *** ASSUMES SCALAR-FIRST CONVENTION: q = [qw, qx, qy, qz] ***
174 %
175 % Inputs p and q MUST be 1x4 row vectors.
176 % Output r is a 1x4 row vector.
177
178 % --- Unpack components of the first quaternion (p) ---
179 % In your case, p = Omega = [0, wx, wy, wz]
180 ps = p(1); % Scalar part (p_s)
181 pv = p(2:4); % Vector part (p_v)
182
183 % --- Unpack components of the second quaternion (q) ---
184 % In your case, q = Orientation = [qw, qx, qy, qz]
185 qs = q(1); % Scalar part (q_s)
186 qv = q(2:4); % Vector part (q_v)
187
188 % --- Calculate the Product Components ---
189
190 % 1. Calculate the NEW Scalar Part (r_s):
191 % r_s = p_s*q_s - dot(p_v, q_v)
192 r_s = ps * qs - dot(pv, qv);
193
194 % 2. Calculate the NEW Vector Part (r_v):
195 % r_v = p_s*q_v + q_s*p_v + cross(p_v, q_v)
196 r_v = ps * qv + qs * pv + cross(pv, qv);
197
198 % Combine product: r = [r_s, r_v] (SCALAR-FIRST order)
199 r = [r_s, r_v];
200 end
201
202
203 %%%%%%%%%%%%%%%%%%%%%QUAT2ROTM FUNCTION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
204 function R_bi = quat2rotm_manual(q)
205 % QUAT2ROTM_MANUAL Converts a quaternion to a 3x3 rotation matrix.
206 % This replaces the need for the Navigation or Robotics toolboxes.
207 %
208 % Input: q - A 4x1 or 1x4 vector [w, x, y, z] (Scalar first)
209 % Output: R_bi - A 3x3 Rotation Matrix
210
211 % 1. Ensure q is a unit quaternion (essential for 30-day stability)
212 %q = q / norm(q); this is redundant and taken care of my Baumgarte
213
214 % 2. Extract components
215 w = q(1);
216 x = q(2);
217 y = q(3);
218 z = q(4);
219
220 % 3. Calculate the rotation matrix (Hamiltonian convention)
221 R_bi = [1 - 2*y^2 - 2*z^2, 2*x*y - 2*w*z, 2*x*z + 2*w*y;
222 2*x*y + 2*w*z, 1 - 2*x^2 - 2*z^2, 2*y*z - 2*w*x;
223 2*x*z - 2*w*y, 2*y*z + 2*w*x, 1 - 2*x^2 - 2*y^2];
224 end
225
226
227 %%%%%%%%%%%%%%%%%MAGNETIC FIELD FUNCTION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
228 function B_eci = getEarthField(r_eci)
229 % --- Constants ---
230 M_earth = 7.94e22; % Earth's magnetic moment (Amp-m^2)
231 mu0_4pi = 1e-7; % Constant (mu0 / 4*pi)
232
233 % --- Distance Calculation ---
234 r_mag = norm(r_eci); % Distance from Earth center (meters)
235 r_hat = r_eci / r_mag; % Unit vector pointing to the satellite
236
237 % --- Magnetic North Pole Orientation ---
238 % For a standard model, the dipole is aligned near the Z-axis
239 % You can refine this later with an 11.5 degree tilt if needed.
240 m_hat = [0; 0; 1];
241
242 % --- The Dipole Equation ---
243 % B = (mu0/4pi) * (3*r_hat*(m_dot_r) - m) / r^3
244 m_dot_r = dot(m_hat, r_hat);
245 B_eci = (mu0_4pi * M_earth / r_mag^3) * (3 * r_hat * m_dot_r - m_hat);
246 end
247
248
249 %%%%%%%%%%%%%%%%%%%ORBITAL PROPAGATOR FUNCTION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
250 function position = circularPropagator(t)
251
252 R = 6371 * 10^3 + 400 * 10^3;
253 i = deg2rad(51.6);
254 mu = 3.985892 * 10^14;
255 n = sqrt(mu / R^3);
256 theta0 = 0;
257 theta = n * t + theta0;
258 x = R * cos(theta);
259 y = R * sin(theta) * cos(i);
260 z = R * sin(theta) * sin(i);
261
262 position = [x; y; z];
263 end
264
265
266 %%%%%%%%%%%%%%%%%%CALCULATE_B_HYST FUNCTION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
267 % function m_hyst = calculate_m_hyst(y, B_body, V_hyst, B_s, H_c, B_r)
268 % y has to be a vertical matrix
269
270 % mu_0 = 4 * pi * 10^(-7);
271 % H = B_body / mu_0;
272 % p = (1 / H_c) * tan((pi * B_r) / (2 * B_s));
273
274 % H_dot = -cross(y(5:7), H); % assuming y(5:7) is your [wx; wy; wz]
275 % find the sign (s) for each axis (x, y, z)
276 % s will be a 3x1 vector of 1, -1, or 0
277
278 % first sign flip method
279 % s = sign(H_dot);
280 % fallback: if sign is 0 (not moving), assume increasing to avoid errors
281 % s(s == 0) = 1;
282
283 % Hdot_scale = 0.05; % A/m/s (tuned to LEO dynamics)
284 % s = tanh(H_dot / Hdot_scale);
285
286 % continuous approximation - better sign flip method
287 % eps = 1e-6; % smoothing parameter
288 % s = H_dot ./ sqrt(H_dot.^2 + eps^2);
289
290 % Hdot_ref = 5e-4; % A/m/s, consistent with Gerhardt test rates
291 % s = tanh(H_dot / Hdot_ref);
292
293 % B_hyst = (2 / pi) * B_s * atan(p * (H - s * H_c));
294
295 % m_hyst = (B_hyst * V_hyst) / mu_0;
296
297 % end
298
299 function m_hyst = calculate_m_hyst(y, B_body, V_hyst, B_s, H_c, B_r)
300 mu_0 = 4 * pi * 10^(-7);
301 H = B_body / mu_0;
302
303 % Calculation of the 'p' parameter for the Henretty model
304 p = (1 / H_c) * tan((pi * B_r) / (2 * B_s));
305
306 % H_dot represents the change in magnetic field relative to the body axes
307 H_dot = -cross(y(5:7), H);
308
309 % IMPROVEMENT: Dynamic Scaling
310 % Using a slightly tighter scale (0.01) ensures the hysteresis loop
311 % closes properly at very low angular velocities.
312 Hdot_scale = 0.01;
313 s = tanh(H_dot / Hdot_scale);
314
315 % ENERGY LOCK: Prevent "Negative Damping"
316 % If the angular velocity is nearly zero, force s to zero to ensure
317 % the rods don't accidentally add energy due to numerical noise.
318 if norm(y(5:7)) < 1e-6
319 s = [0; 0; 0];
320 end
321
322 % The Henretty sigmoid curve equation
323 B_hyst = (2 / pi) * B_s * atan(p * (H - s * H_c));
324 m_hyst = (B_hyst * V_hyst) / mu_0;
325 end