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 ...

Code:
% 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 ...

Code:
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

Code:
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

Code:
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

Code:
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 .