Planetary heatmap creation in matlab

[Naima]

Hello , I am trying to make a planet heatmap on matlab but I woul dlike to have some outside check for formulas and conditions I gived as assumptions ...

% Load the saved colored heightmap data

clear all
load('colored_heightmap.mat', 'colored_heightmap');

[Orbital_Period, axial_tilt, day_of_year, months_in_year, days_in_month, To, Tp, Tm, S, Alb, Emiss, Eccentricity, Mean_Anomaly, perihelion_day] = heatmapwithlatitudeanddensityair_inputs();

% Calculate the heatmap of the planet
n = size(colored_heightmap, 1); % number of rows in the map
m = size(colored_heightmap, 2); % number of columns in the map

% 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

T = zeros(n, m); % temperature of each element in °C
for latitude_index = 1:n
for longitude_index = 1:m
color = colored_heightmap(latitude_index, longitude_index);
if color == 0 % blue (ocean)
T(latitude_index, longitude_index) = To;
elseif color == 1 % green (plains)
T(latitude_index, longitude_index) = Tp;
elseif color == 2 % white (mountains)
T(latitude_index, longitude_index) = Tm;
end
end
end

% Calculate the heat balance of the planet
Theta = angle_of_incidence_Function(day_of_year, axial_tilt, Positional_Matrix, Orbital_Period, Eccentricity, Mean_Anomaly, perihelion_day); % update the angle of incidence
rho = 1.225 ; % atmospheric density
P = 1; % atmospheric pressure
H = S * (1 - Alb) * cos(Theta) / (4 * Emiss) * rho * P; % heat balance in W/m^2
T = T + H; % update the temperature of each element in °C

% Plot the original map
plot_original_map(colored_heightmap);

% Plot the heatmap
plot_heat_map(T, Positional_Matrix, colored_heightmap, months_in_year, days_in_month, axial_tilt, Orbital_Period, S, Alb, Emiss, rho, P,Eccentricity,Mean_Anomaly, perihelion_day);

% Set the same scale for both maps
axis([1 m 1 n]);

then the functions to load

calculate angle of incidence of sun , hope someone could tell me if calculations math is correct too ...

function Theta = angle_of_incidence_Function(day_of_year, axial_tilt, Positional_Matrix, Orbital_Period, Eccentricity, Mean_Anomaly,perihelion_day)
% Calculate the angle of incidence of the heat from the sun in radians

% 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));
end

function with static parameters

function [Orbital_Period, axial_tilt, day_of_year, months_in_year, days_in_month, To, Tp, Tm, S, Alb, Emiss, Eccentricity, Mean_Anomaly, perihelion_day] = heatmapwithlatitudeanddensityair_inputs()
 

% Define the average temperature of the ocean, the plains, and the mountains
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 ;


% Define the solar radiation constant, the albedo, and the emissivity
S = 1366; % solar radiation constant in W/m^2
Alb = 0.3; % albedo
Emiss = 0.6; % emissivity

end

plot heat map function

function plot_heat_map(T, Positional_Matrix, colored_heightmap, months_in_year, days_in_month, axial_tilt, Orbital_Period, S, Alb, Emiss, rho, P,Eccentricity,Mean_Anomaly, perihelion_day)
    months = {'January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'};
    figure;
    colormap('parula');
    % Match the proportions of the x and y axis to the ones of the colored_heightmap
    [n, m] = size(colored_heightmap);
    imagesc([1 m], [1 n], T);
    colorbar;
    
    % Loop through each month of the year
    for month = 1:months_in_year
        day_of_year = (month - 1) * days_in_month + 1; % update the day of the year
        Theta = angle_of_incidence_Function(day_of_year, axial_tilt, Positional_Matrix, Orbital_Period, Eccentricity, Mean_Anomaly, perihelion_day); % update the angle of incidence
 % update the angle of incidence
        H = S * (1 - Alb) * cos(Theta)  / (4 * Emiss) * rho * P; % update the heat balance
        T = T + H; % update the temperature of each element 
        imagesc(T); % update the heatmap
        title(['Heatmap - ' months{month}]);
        yticks([1 size(colored_heightmap, 1)/6 size(colored_heightmap, 1)/3 size(colored_heightmap, 1)/2 size(colored_heightmap, 1)*2/3 size(colored_heightmap, 1)*5/6 size(colored_heightmap, 1)])
        yticklabels({'-90','-60','-30','0','30','60','90'})
        xticks([1 size(colored_heightmap, 2)/6 size(colored_heightmap, 2)/3 size(colored_heightmap, 2)/2 size(colored_heightmap, 2)*2/3 size(colored_heightmap, 2)*5/6 size(colored_heightmap, 2)])
        xticklabels({'-180','-120','-60','0','60','120','180'})
        pause(1); % pause for 1 second
    end
end

plot original map

function plot_original_map(colored_heightmap)
figure;
imagesc(colored_heightmap);
colormap([0 0 1; 0 1 0; 1 1 1]);

axis tight;
axis xy;
xlim([1 size(colored_heightmap, 2)]);
ylim([1 size(colored_heightmap, 1)]);
axis equal;
title("Original Map");

yticks([1 size(colored_heightmap, 1)/6 size(colored_heightmap, 1)/3 size(colored_heightmap, 1)/2 size(colored_heightmap, 1)*2/3 size(colored_heightmap, 1)*5/6 size(colored_heightmap, 1)])
yticklabels({'-90','-60','-30','0','30','60','90'})
xticks([1 size(colored_heightmap, 2)/6 size(colored_heightmap, 2)/3 size(colored_heightmap, 2)/2 size(colored_heightmap, 2)*2/3 size(colored_heightmap, 2)*5/6 size(colored_heightmap, 2)])
xticklabels({'-180','-120','-60','0','60','120','180'})
end

and I added the map in link

Thanks for any help .

https://files.fm/u/46euuevrp

https://files.fm/u/fcs5s8q6k



Issues I have, :

the heatmap not showing correct values of heat and difference among sea and land,

size is not matching the original one

the representation of the heat seems wrong as the winter in north emisphere looks like its equal and when its june it becomes hotter then reverts back to normal but doesn't get cold .

[waldronate]

It looks like you're calculating the temperature at the top of the atmosphere (or for an airless world), which can be a useful parameter to input into a climate model. I didn't try to verify your math because that sort of thing is well-defined in a lot of places already and I'm not up to that level of detail work at the moment, but make sure that you check your units.

Switching the base temperature based on the underlying land type is certainly one way to go about getting a temperature difference between sea, land, and mountainous provinces, but it won't get you any gradients between those values. If you have a heightmap, using the altitude as an input to a temperature function will get you the underlying gradients of your heightmap.

There is still a long way to go (ITCZ calculations, currents in the atmosphere, heat flow around the oceans and atmosphere along with exchanges between them, cloud albedo effects, greenhouse effects, and so on) to get to a minimal climate model, but solar influx at the top of the atmosphere is a good starting point.

[Naima]

yes the idea was to model independent systems and then couple together for interaction , I wanted to start first with the heat transfer and balance of energy to determine the general heating and surface temperature, I yet have to implement differential types of emissivity for land , sea and in future ice too , to keep things simple, I didn't want to model 3dimensional verticality in atmosphere or in ocean depth , eventually I might add a a second layer of shallow water and deep water , butyet to be seen , mathlab is giving me headakes to make it work as I want and I am not understaiding yet how to fix the issue that the north emisphere doesnt get colder in winter and so warmer the southern emisphere, somethig is wrong .

[waldronate]

mathlab is giving me headakes to make it work as I want

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.

[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;

[waldronate]

I understood generally what the code was doing from the variable names, but I don't know why you wanted to do those things. When writing software, I believe that comments should describe why you're doing something rather than describe what the code is doing because the code itself always says what it's doing.

I'm interested in the meanings of the physical quantities behind the terms "solar declination" and "angle of incidence" as used here and why you are using them in the manner shown. When I did these sorts of calculations last century, I used some simplifying assumptions like a circular orbit, constant orbital phase, and a day much shorter than the orbital period to get some very simple results. I came up with a spreadsheet that gave the set of graphs below for the traditional seasons:

insolation.png

They aren't "accurate" by any stretch of the imagination, but certainly were good enough for my purposes at the time, especially considering the coarse level of fidelity in what I was doing elsewhere. Using Earth as an example, the 4% variance in energy between perihelion and aphelion is dwarfed by the error of sampling the orbit at only four points.

As far as tooling, I used a spreadsheet for this because the calculations are mostly a function of two values (tilt and latitude), which is really easy to do on a spreadsheet. And it's what I had at hand way back when, of course.

[Naima]

The idea would like to be to create something simlar to this: http://www.buildyourownearth.com/byo...2=0&c2=39&v=pm
Without implementing the whole GCM of layers and strata , eventually do it a simplier but kind of acceptable for information.

The comments I put are there for me to remember what is what .
Slar Declination shoudl stand for the declination angle, the angle between the equator and a line drawn from the centre of the Earth to the centre of the sun.
then the incidence of radiation angle and the tilting angle, I dunno if I messed something may be in the math there ...

cos(ζ)=sin(ψ)sin(δ)+cos(ψ)cos(δ)cos(h)

file:///C:/Users/Stefano/Downloads/850-Article%20Text-895-1-10-20190628.pdf

So perhaps the problem is in the equations I implemented?

file:///C:/Users/Stefano/Downloads/Acomprehensivesolaranglessimulationandcalculationu singMatlab.pdf

[waldronate]

I can't read those PDFs because you linked to files on your local filesystem, not on the internet.

One way to figure out the problem is with the divide and conquer method. Make sure each part is doing what you expect it to do before you combine them. Maybe even simplify a little.

Some examples: make sure that the Positional_Matrix has the values you expect. What you have here is the same value in every column with no variance in rows, so you might be able to simplify it to be an array until you're sure that the math is right (whether that would help in this context is questionable, but a 2D graph is often easier to understand then a 3D one).
Make sure that the angle and distance orbital calculations are varying smoothly. It's a little harder with an elliptical orbit than a circular one, but not too much.
Make sure that the angle of the sun at each latitude for each season/month varies smoothly from pole to pole (roughly cos(axial_tilt-latitude) to -cos(axial_tilt+latitude) where axial_tilt is the tilt angle pointed toward the sun).

Only after you've gotten to the point of having your insolation at the top of the atmosphere calculation working should you move on to trying to get things like surface heating.

[Naima]

I can't read those PDFs because you linked to files on your local filesystem, not on the internet.

One way to figure out the problem is with the divide and conquer method. Make sure each part is doing what you expect it to do before you combine them. Maybe even simplify a little.

Some examples: make sure that the Positional_Matrix has the values you expect. What you have here is the same value in every column with no variance in rows, so you might be able to simplify it to be an array until you're sure that the math is right (whether that would help in this context is questionable, but a 2D graph is often easier to understand then a 3D one).
Make sure that the angle and distance orbital calculations are varying smoothly. It's a little harder with an elliptical orbit than a circular one, but not too much.
Make sure that the angle of the sun at each latitude for each season/month varies smoothly from pole to pole (roughly cos(axial_tilt-latitude) to -cos(axial_tilt+latitude) where axial_tilt is the tilt angle pointed toward the sun).

Only after you've gotten to the point of having your insolation at the top of the atmosphere calculation working should you move on to trying to get things like surface heating.

Yes sorry I thought to link the online resource and I linked my downoaded one ...

https://files.fm/u/ftmchffaf

Anyway I have tested the matrix and I got what I think correct values (-1.57 to 1.57 ) but not sure why this is not represented in the change , The other values for longitude should be zeros if not wrong. may be I should implement some other kind of checks, basically my checks were to numerically print the values in the command window.