From fe302e0a3dc01b770181f500c425337fcd4036e5 Mon Sep 17 00:00:00 2001
From: Marcodagostino2 <mcdago25@gmail.com>
Date: Sat, 16 Nov 2024 14:30:38 +0100
Subject: [PATCH] [design-updates] Fixed model gravity. Gravity model use a
 point-mass gravitation model with an Earth radius defined in the Environment
 class.

---
 classes/components/Environment.m   | 62 +++++++++---------------------
 functions/odeFunctions/ballistic.m |  6 +--
 2 files changed, 20 insertions(+), 48 deletions(-)

diff --git a/classes/components/Environment.m b/classes/components/Environment.m
index c0c5cd9..52ba0b1 100644
--- a/classes/components/Environment.m
+++ b/classes/components/Environment.m
@@ -12,29 +12,28 @@ classdef Environment < Component
 %               file
 
     properties
-        lat0            double                % [deg] Launchpad latitude
-        lon0            double                % [deg] Launchpad longitude
-        z0              double                % [m] Launchpad Altitude
-        omega           double                % [deg] Launchpad Elevation
-        phi             double                % [deg] Launchpad Azimuth
-        pin1Length      double                % [m] Distance from the upper pin to the upper tank cap
-        pin2Length      double                % [m] Distance from the lower pin to the lower tank cap
-        rampLength      double                % [m] Total launchpad length
+        lat0            double                  % [deg] Launchpad latitude
+        lon0            double                  % [deg] Launchpad longitude
+        z0              double                  % [m] Launchpad Altitude
+        omega           double                  % [deg] Launchpad Elevation
+        phi             double                  % [deg] Launchpad Azimuth
+        pin1Length      double                  % [m] Distance from the upper pin to the upper tank cap
+        pin2Length      double                  % [m] Distance from the lower pin to the lower tank cap
+        rampLength      double                  % [m] Total launchpad length
 
-        temperature     double = 288.15       % [K] Ground temperature
-        pressure        double                % [Pa] Ground pressure
-        rho             double                % [Kg/m^3] Ground air density
-        gamma           double = 1.4          % [-] Gas constant
+        temperature     double = 288.15         % [K] Ground temperature
+        pressure        double                  % [Pa] Ground pressure
+        rho             double                  % [Kg/m^3] Ground air density
+        gamma           double = 1.4            % [-] Gas constant
     end
 
     properties(SetAccess = private)
-        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
-        effectiveRampAltitude   double        % [m] Projection of effectiveRampLength on z axis
+        g0                      double              % [-] Gravity costant at launch latitude and altitude
+        earthRadius             double = 6371007    % [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
+        effectiveRampAltitude   double              % [m] Projection of effectiveRampLength on z axis
     end
 
     properties(Access = protected)
@@ -66,8 +65,6 @@ classdef Environment < Component
     methods
         function updateAll(obj)
             obj.updateG0;
-            obj.updateG0Centrifugal;
-            obj.updateEarthRadius;
             obj.updatePinDistance;
             obj.updateRamp;
             obj.updateLocal;
@@ -81,28 +78,7 @@ classdef Environment < Component
 
         function obj = updateG0(obj)
             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);
+                obj.g0 = gravitywgs84(obj.z0, obj.lat0);
             end
         end
 
diff --git a/functions/odeFunctions/ballistic.m b/functions/odeFunctions/ballistic.m
index 391f3c3..cce9648 100644
--- a/functions/odeFunctions/ballistic.m
+++ b/functions/odeFunctions/ballistic.m
@@ -58,11 +58,7 @@ dY = zeros(14, 1);
 %% GETTING DATA
 S = rocket.crossSection;                                                    % [m^2]   cross surface
 C = rocket.diameter;                                                        % [m]     caliber
-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
+g = environment.g0/(1 + (altitude/environment.earthRadius))^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