diff --git a/classes/components/Environment.m b/classes/components/Environment.m index 7ee1a91d3532686e58dbb53ce4c369a3f1e907ca..c0c5cd98a47f9a85b479ae6a64526dc10db4d427 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 0194eff2564ba79ba00e30e47181915a47a14891..391f3c366d3b40b5668dae44678659c9e4b7b6ce 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