Results 1 to 9 of 9

Thread: Planetary heatmap creation in matlab

  1. #1
    Professional Artist Naima's Avatar
    Join Date
    Mar 2010
    Location
    Italy
    Posts
    1,556

    Default Planetary heatmap creation in matlab

    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 .

  2. #2
    Administrator waldronate's Avatar
    Join Date
    Mar 2007
    Location
    The High Desert
    Posts
    3,549

    Default

    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.

  3. #3
    Professional Artist Naima's Avatar
    Join Date
    Mar 2010
    Location
    Italy
    Posts
    1,556

    Default

    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 .

  4. #4
    Administrator waldronate's Avatar
    Join Date
    Mar 2007
    Location
    The High Desert
    Posts
    3,549

    Default

    Quote Originally Posted by Naima View Post
    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.

  5. #5
    Professional Artist Naima's Avatar
    Join Date
    Mar 2010
    Location
    Italy
    Posts
    1,556

    Default

    Quote Originally Posted by waldronate View Post
    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 :
    Code:
    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
    Code:
    "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
    Code:
    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

    Code:
    Eccentric_Anomaly = 2 * atan(sqrt((1 - Eccentricity) / (1 + Eccentricity)) * tan(Mean_Anomaly / 2));
    % Calculate the solar declination

    Code:
    delta = asin(sin(axial_tilt) * sin(2 * pi * (Eccentric_Anomaly - day_of_year) / Orbital_Period));
    % Calculate the angle of incidence

    Code:
    Theta = acos(sin(Positional_Matrix) * sin(delta) + cos(Positional_Matrix) * cos(delta));
    % Calculate the heat balance of the planet

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

  6. #6
    Administrator waldronate's Avatar
    Join Date
    Mar 2007
    Location
    The High Desert
    Posts
    3,549

    Default

    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:

    Click image for larger version. 

Name:	insolation.png 
Views:	10 
Size:	53.5 KB 
ID:	135196

    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.

  7. #7
    Professional Artist Naima's Avatar
    Join Date
    Mar 2010
    Location
    Italy
    Posts
    1,556

    Default

    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
    Last edited by Naima; 02-09-2023 at 02:56 AM.

  8. #8
    Administrator waldronate's Avatar
    Join Date
    Mar 2007
    Location
    The High Desert
    Posts
    3,549

    Default

    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.

  9. #9
    Professional Artist Naima's Avatar
    Join Date
    Mar 2010
    Location
    Italy
    Posts
    1,556

    Default

    Quote Originally Posted by waldronate View Post
    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.

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •