Planetary heatmap creation in matlab

[Naima]

I believe that is the primary functionality of matlab, yes.

I'm not sure what physical quantity "angle_of_incidence_Function" is supposed to be calculating. Can you elaborate on what it is and how you got to those equations? I recognize that they are similar to the ones given on the Wikipedia "eccentric anomaly" page for angle and radius, but they are a bit different as the radius needs to be calculated from the ellipse center rather than one of the foci.

I vaguely recall reading once that Earthlike worlds probably don't occur with significantly ellipsoidal orbits because it messes with the atmospheric processes (or maybe it was excessive temperature swings that would make life unlikely). But that would have been a very long time ago and just for that one paper, so I may be recalling incorrectly.

here are the parameters I have used to calculate the angle of incidence of the heat :
Some fixed params :

rho = 1.225 ; % atmospheric density
P = 1; % atmospheric pressure
To = 20; % average temperature of ocean in °C
Tp = 15; % average temperature of plains in °C
Tm = 5; % average temperature of mountains in °C
Orbital_Period = 365; % Define the number of days in the year
axial_tilt = (23 * pi) / 180; % Define planet's axial tilt
day_of_year = 1; % any value between 1 and 365 
months_in_year = 12;
days_in_month = 365 / months_in_year;
Eccentricity = 0.0167086;
perihelion_day = 3;
Mean_Anomaly  = 2 * pi * (day_of_year - perihelion_day) / Orbital_Period

%some variable naming

"day_of_year" - the day of the year
"axial_tilt" - the axial tilt of the object
"Positional_Matrix" - the positional matrix latitude parameter on the grid map.
"Orbital_Period" - the orbital period of the object
"Eccentricity" - the eccentricity of the object's orbit
"Mean_Anomaly" - the mean anomaly of the object
"perihelion_day" - the perihelion day of the object.

% Define the latitude of each element in the map

Positional_Matrix = zeros(n, m); % latitude in radians
for latitude_index = 1:n
for longitude_index = 1:m
Positional_Matrix(latitude_index, longitude_index) = ((n/2 - latitude_index ) * pi / n);  %from  -π/2 to  π/2 in radiants (-1.57 to 1.57 )
end
end

% Calculate the eccentric anomaly

Eccentric_Anomaly = 2 * atan(sqrt((1 - Eccentricity) / (1 + Eccentricity)) * tan(Mean_Anomaly / 2));

% Calculate the solar declination

delta = asin(sin(axial_tilt) * sin(2 * pi * (Eccentric_Anomaly - day_of_year) / Orbital_Period));

% Calculate the angle of incidence

Theta = acos(sin(Positional_Matrix) * sin(delta) + cos(Positional_Matrix) * cos(delta));

% Calculate the heat balance of the planet

H = S * (1 - Alb) * cos(Theta) / (4 * Emiss) * rho * P; % heat balance in W/m^2;