From 5a9ca0774746062801a3c887d97a22b3970e5005 Mon Sep 17 00:00:00 2001
From: Marcodagostino2 <mcdago25@gmail.com>
Date: Fri, 15 Nov 2024 21:33:00 +0100
Subject: [PATCH] [design-updates] Fixed gravity model Added variable earth's
radius in accordance with the position on the earth's ellipsoid. Updated the
model in g computation.
---
classes/components/Environment.m | 31 ++++++++++++++++++++++++++++--
functions/odeFunctions/ballistic.m | 6 +++++-
2 files changed, 34 insertions(+), 3 deletions(-)
diff --git a/classes/components/Environment.m b/classes/components/Environment.m
index 7ee1a91..c0c5cd9 100644
--- a/classes/components/Environment.m
+++ b/classes/components/Environment.m
@@ -28,7 +28,9 @@ classdef Environment < Component
end
properties(SetAccess = private)
- g0 double % [-] Gravity costant at launch latitude and altitude
+ g0 double % [-] Gravity costant at launch latitude and z = 0m
+ g0Centrifugal double % [m/s^2] Gravity due to earth rotation at launch latitude and z = 0m
+ earthRadius double % [m] Radius used for g computation
local double % [-] Vector conatining inputs for atmosphereData
pinDistance double % [m] Distance of the upper pin from the rail base (upper pin-boat + boat-rail base)
effectiveRampLength double % [m] Total launchpad length
@@ -64,6 +66,8 @@ classdef Environment < Component
methods
function updateAll(obj)
obj.updateG0;
+ obj.updateG0Centrifugal;
+ obj.updateEarthRadius;
obj.updatePinDistance;
obj.updateRamp;
obj.updateLocal;
@@ -76,7 +80,30 @@ classdef Environment < Component
end
function obj = updateG0(obj)
- obj.g0 = gravitywgs84(obj.z0, obj.lat0);
+ if(~isempty(obj.lat0))
+ obj.g0 = gravitywgs84(0, obj.lat0);
+ end
+ end
+
+ function obj = updateG0Centrifugal(obj)
+ if(~isempty(obj.lat0))
+ earthOmega = 2*pi/86164; % [rad/s] Earth angular velocity
+ aEarth = 6378122; % [m] Earth semi major axis (radius at the equator)
+ eEarth = 0.08182; % [/] Earth's ellipsoid eccentricity
+ lat0rad = deg2rad(obj.lat0);
+ obj.g0Centrifugal = - earthOmega^2 * aEarth * ((cos(lat0rad)^2)/(sqrt(1 - eEarth^2*(sin(lat0rad)^2))));
+ end
+ end
+
+ function obj = updateEarthRadius(obj)
+ if(~isempty(obj.lat0))
+ earthOmega = 2*pi/86164;
+ R0 = 6370000;
+ lat0rad = deg2rad(obj.lat0);
+ mu = 3.986e14; % [m^3/s^2] Earth's gravitational constant (GM)
+ f = @(R) obj.g0 - (mu/(R^2)) + earthOmega^2*R*(cos(lat0rad)^2);
+ obj.earthRadius = fzero(f, R0);
+ end
end
function obj = updatePinDistance(obj)
diff --git a/functions/odeFunctions/ballistic.m b/functions/odeFunctions/ballistic.m
index 0194eff..391f3c3 100644
--- a/functions/odeFunctions/ballistic.m
+++ b/functions/odeFunctions/ballistic.m
@@ -58,7 +58,11 @@ dY = zeros(14, 1);
%% GETTING DATA
S = rocket.crossSection; % [m^2] cross surface
C = rocket.diameter; % [m] caliber
-g = environment.g0/(1 + (altitude*1e-3/6371))^2; % [N/kg] module of gravitational field
+earthOmega = 2*pi/86164; % [rad/s] Earth's angular velocity
+g = (environment.g0 - environment.g0Centrifugal)/...
+ ((1 + ((altitude + environment.z0)/environment.earthRadius))^2) + ...
+ environment.g0Centrifugal - earthOmega^2*...
+ (altitude + environment.z0)*(cos(environment.lat0*pi/180)^2); % [N/kg] module of gravitational field
tb = rocket.motor.cutoffTime; % [s] Burning Time
theta = t-t0; % [s] Time shift (needed for control)
local = environment.local; % vector containing inputs for atmosphereData
--
GitLab