Added ability to compute atmo density/pressure/temp from interpolant.…

… Call function from main KSPTOT GUI -> Edit menu.

helper_methods/astrodynamics/getPositOfBodyWRTSun.m

@@ -5,11 +5,11 @@
         celBodyData
     end
 
-%     if(numel(time) == 1 && bodyInfo.lastComputedSunTime == time)
-%         rVectB = bodyInfo.lastComputedSunRVect;
-%         vVectB =  bodyInfo.lastComputedSunVVect;
-%         return;
-%     end
+    % if(isscalar(time) && bodyInfo.lastComputedSunTime == time)
+    %     rVectB = bodyInfo.lastComputedSunRVect;
+    %     vVectB =  bodyInfo.lastComputedSunVVect;
+    %     return;
+    % end
 
     try
         if(bodyInfo.propTypeIsTwoBody || (isscalar(time) && time == bodyInfo.epoch))

helper_methods/ksptot_ma/propagation/aerobrake/getAtmoDensityAtAltitude.m

@@ -2,24 +2,42 @@
 %getAtmoDensityAtAltitude Summary of this function goes here
 %   Detailed explanation goes here
 %   See here for math: https://forum.kerbalspaceprogram.com/index.php?/topic/142686-modeling-atmospheres-in-ksp/
+    arguments
+        bodyInfo(1,1) KSPTOT_BodyInfo
+        altitude(1,1) double
+        lat(1,1) double
+        ut(1,1) double
+        long(1,1) double
+    end
 
     if(altitude <= bodyInfo.atmohgt && altitude >= 0)
-%         pressure = getPressureAtAltitude(bodyInfo, altitude);
-        pressure = getBodyAtmoPressure(bodyInfo, altitude);
-
-        if(pressure > 0)
-            temperature = getBodyAtmoTemperature(bodyInfo, ut, lat, long, altitude);
+        if(not(isempty(bodyInfo.densityGI)))
+            lat = angleNegPiToPi_mex(lat);
+            long = AngleZero2Pi(long);
+
+            out = bodyInfo.densityGI(lat, long, altitude);
+            % out = squeeze(out);
 
-            density = getDensityFromIdealGasLaw(pressure, temperature, bodyInfo.atmomolarmass);
-            if(density < 0)
+            density = out(1);
+            pressure = out(2);
+            temperature = out(3);
+        else
+            pressure = getBodyAtmoPressure(bodyInfo, altitude);
+
+            if(pressure > 0)
+                temperature = getBodyAtmoTemperature(bodyInfo, ut, lat, long, altitude);
+
+                density = getDensityFromIdealGasLaw(pressure, temperature, bodyInfo.atmomolarmass);
+                if(density < 0)
+                    density = 0;
+                    pressure = 0;
+                    temperature = 0;
+                end
+            else
                 density = 0;
                 pressure = 0;
                 temperature = 0;
             end
-        else
-            density = 0;
-            pressure = 0;
-            temperature = 0;
         end
     elseif(altitude <= 0)
         density = 0;

helper_methods/zz_classes/@KSPTOT_BodyInfo/KSPTOT_BodyInfo.m

@@ -27,6 +27,8 @@
         axialtempsunmultcurve = griddedInterpolant([0 1]',[0 0]','spline', 'nearest');
         ecctempbiascurve = griddedInterpolant([0 1]',[0 0]','spline', 'nearest');
         atmomolarmass double
+
+        densityGI = []; %atmo density/pressure/temp gridded interp <- Based on the above curves and a UT = 0
 
         %Rotation
         rotperiod double
@@ -621,6 +623,51 @@ function setBaseBodyHeightMap(obj)
             end
         end
     end
+
+    methods
+        function createDensityGriddedInterp(obj)
+            arguments
+                obj(1,1) KSPTOT_BodyInfo
+            end
+
+            lats = linspace(-pi/2,pi/2,100);
+            longs = linspace(0,2*pi,100);
+            alts = linspace(0,obj.atmohgt,100);
+
+            combos = combvec(lats, longs, alts)';
+            numRuns = height(combos);
+
+            density = NaN(numRuns,1);
+            pressureKPA = NaN(numRuns,1);
+            temperature = NaN(numRuns,1);
+            ut = 0;
+
+            pp=gcp('nocreate');
+            if(isempty(pp))
+                M = 0;
+            else
+                M = pp.NumWorkers;
+            end
+
+            obj.densityGI = []; %clears the existing interpolant so that the actual functions get called
+            parfor(i=1:numRuns, M)
+                lat = combos(i,1);
+                long = combos(i,2);
+                altitude = combos(i,3);
+
+                [density(i), pressureKPA(i), temperature(i)] = getAtmoDensityAtAltitude(obj, altitude, lat, ut, long);
+            end
+
+            densityRS = reshape(density, [length(lats), length(longs), length(alts)]);
+            pressureKPA_RS = reshape(pressureKPA, [length(lats), length(longs), length(alts)]);
+            temperatureRS = reshape(temperature, [length(lats), length(longs), length(alts)]);
+            V(:,:,:,1) = densityRS;
+            V(:,:,:,2) = pressureKPA_RS;
+            V(:,:,:,3) = temperatureRS;
+
+            obj.densityGI = griddedInterpolant({lats',longs',alts'}, V, "linear","nearest");
+        end
+    end
 
 	methods(Static)
         function obj = loadobj(obj)