Click images for more details



Recent comments
Recent posts
Currently discussing

A few sites I've stumbled across recently....

Powered by Squarespace

Discussion > Numerical calculation of Planetary Black Body Temperature

Having 5 minutes spare, I wrote the following noddy program to test what temperature you might get if you used a numerical solution to working out average temperature (instead of averages):

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace NandZ
   class Program
      static double SBconstant = 5.670373e-8;
      static double backgroundTemp = 2.72;

      static void Main(string[] args)
         CalcTemp(1366.0, 0.3, 0.95, "Earth");
         CalcTemp(1366.0, 0.12, 0.95, "Moon");

      static void CalcTemp(double TotalSolar, double Albedo, double Emissivity, string name)
         int steps = 10000;
         double allT = 0;
         int count = 0;

         for (double lat = 0f; lat <= Math.PI; lat = lat + (Math.PI / steps))
            for (double lon = 0f; lon <= Math.PI; lon = lon + (Math.PI / steps))
               //lit side
               double T = backgroundTemp + Math.Pow(((1 - Albedo) * TotalSolar * Math.Sin(lon) * Math.Sin(lat)) / (Emissivity * SBconstant), 0.25);
               allT += T;
               // dark side
               allT += backgroundTemp;

         double avgT = allT / count ;
         Console.WriteLine("{0} Average Temp Without Atmosphere: {1:f1} K", name, avgT);



When run this gives the following output

Earth Average Temp Without Atmosphere: 137.5 K
Moon Average Temp Without Atmosphere: 145.4 K

Which appears to approximately endorse N&Z's double integration of average temperature, and means that the earth is being warmed by around 150K and not the popular 33K. This is interetsing, but I'd like someone else to look over my code for howlers in code and/or technique.

Anyone spot any problems?

Jan 7, 2013 at 3:51 PM | Unregistered CommenterTheBigYinJames

On the way home I realised that this technique gives too much weight to the polar regions, just like a Mercator projection map

I may have to introduce a weighting factor of sin(longitude) to my count variable (instead of using 1 throughout) which will reduce the impact of the closer together polar areas.

I'll try it when I get a chance.

Jan 7, 2013 at 4:53 PM | Unregistered CommenterTheBigYinJames

BigYin, if you are not rotating, you are actually radiating from half the sphere as soon as you reach equilibrium. That calc only works if no part of the dark side ever reaches its asymptote (OK, near to it). And probably if no part on the warm side reaches SB temp. You'd need to take a slice and work out the incoming over the daily cycle, then integrate for latitude AND time of day, with provision for up and downside limits.. Then tell me if the daylength cancels out or stays in

Jan 7, 2013 at 5:04 PM | Registered Commenterrhoda

It occurs, without running a check, that if you had half as much radiating area you'd increase temps by root two. And get around 200K. Or is it fourth root of two? Or is it two?

Jan 7, 2013 at 5:07 PM | Registered Commenterrhoda

Thanks Rhoda, I was considering this for a next step once I was happy I was modelling the static un-rotating condition properly.

Currently in the middle of an email dialogue with Martin A trying to get a handle of exactly what rotation will do the temps at each point.

Jan 7, 2013 at 5:08 PM | Unregistered CommenterTheBigYinJames

Ouch, by adding the polar weighting, the temps become:

Earth Average Temp Without Atmosphere: 216.0 K
Moon Average Temp Without Atmosphere: 228.5 K

Which puts earth GHE effect at a more modest 72K.

Jan 7, 2013 at 5:22 PM | Unregistered CommenterTheBigYinJames


shouldn't you use Cos(latitude)?

Jan 7, 2013 at 5:30 PM | Unregistered CommenterPaul Dennis


your results now look closer to mine. I used the insolation averaged over a daily cycle which, I think, may lead to higher estimated S-B average temperatures. My reason for using the average was that there is heat flow into and out of the regolith, thus a component of the energy balance occurs by radiation loss during the lunar night. In the static model everything on the dark side is at the background temperature. This creates a very large surface area at 2K that probably has a significant effect on the average temperature

Jan 7, 2013 at 5:37 PM | Unregistered CommenterPaul Dennis


I would have used cosine if I was measuring my angle as a latitude. But since I'm iterating from 0 to pi in both directions, sin is fine for weights tending to zero near the top and bottom.

I will have to consider a time vector when doing rotation,and should really model the ecliptic tilt.

My worry is what is the starting condition for a daily insolation cycle, surely if we are taking regolith emission into account, I have to start off with some warm rock in the darkside? Or do I just run it for a few hundred cycles and hope it settles down?

Jan 7, 2013 at 5:52 PM | Unregistered CommenterTheBigYinJames

OK, I see I am mistaken about the radiating area, ignore my last.

Paul, not quiite getting what you are saying, but surely that method if I understand it correctly doesn't allow for a difference between no/slow rotation when the SB temp is maxed out compared to fast, where it is only being approached. But as I said on the other thread, if your max temp is circa 394K at one point, and your minimum is 2K for a hemisphere, average cannot be more than 198K and is likely to be much less. Have I got the SB temp for full sun wrong? No albedo, black body, no conduction.

Jan 7, 2013 at 5:57 PM | Registered Commenterrhoda

Oh, and BigYin, kudos for posting hastily-written code for all to see. Brave.

Jan 7, 2013 at 5:58 PM | Registered Commenterrhoda


I see what you've done now. Yes sin is fine as your integrating from pole to pole and between the illuminated edges.

I'm not sure how to make the model more sophisticated as the next steps are to model the thermal properties of regolith and underlying rocks.

My first attempt, rather like yours in approach but used average daily insolation as a function of latitude (Integral [0-Pi] (cos (lat) *sin (lon)) d(lon) / 2*Pi). If my written equation makes since. This simply reduces to cos(lat)/Pi

Using this I was gratified to find that my estimate of the S-B temperature in the equatorial region is not dissimilar to the temperature reported for the isothermal point just beneath the surface.

Jan 7, 2013 at 6:06 PM | Unregistered CommenterPaul Dennis


I can't fault your logic in what you have written. I think the ambiguity arises from the T^4 dependency of the radiated energy. If you determine an average insolation as applied to the whole sphere because of rotation you obviously reduce the max temperature but the temperature is spread out over a wider area. Consider the equatorial belt with a max temperature on the spot facing the sun for a non rotating sphere of 394K (your calculation was right). The temperatures drop off as to both illuminated edges and are effectively zero (well 2K) on the dark side. The average temperature of this equatorial belt has to be somewhat less than 198K.

However, if you allow the sphere to rotate and compute the average energy as 1366/Pi watts/m^2 then the temperature at the equatorial belt is 270K everywhere!

Does this help, or am I talking gobbledy gook?

Jan 7, 2013 at 6:17 PM | Unregistered CommenterPaul Dennis

I'm still not completely convinced that a rotation will make a major difference to an average temperature.

Martin A and rhoda are sure that a slowly spinning planet will have a different average temperature to one spinning more quickly, while I believe all that would happen is the differential between average lit and unlit temperatures would be less for a more rapidly spinning planet.

Martin's reasoning is that a planet with a lit side at a higher temperature (clower rortaion) loses energy at a rate proportional to temperature to the 4th power of a higher temperature. I have pointed out this effect will be limited by thermal equilibrium, and he agrees that under a certain rotational speed this effect tend to zero.

I await developments.

Jan 7, 2013 at 6:19 PM | Unregistered CommenterTheBigYinJames

I've just noticed that since I included emissivity estimates, this thread should actually be Numerical calculation of Planetary Grey Body Temp, but you get the idea.

Paul, would you mind letting me have your temp figures for earth and or moon, so I can work out where we differ.

I'd also like to try the 'washer' technique, and also polar weighting by reducing the number of points, rather than weighting them. See if any of them are different.

Jan 7, 2013 at 6:37 PM | Unregistered CommenterTheBigYinJames


for the moon I get 267K and for Earth I got an average temperature of 252K. I used an albedo of 0.12 (moon), 0.3 (Earth) and an emissivity of 1 for both the moon and Earth. I know these are just approximations but for the moon the Areveda paper estimates average albedos of about 0.12 (for normal incidence and decreasing considerably at lower angles of incidence which is not included in my model) and an emissivity of 0.98.

I don't think I'll get time tonight to look at this in any more detail but good luck with your estimates. I'm not sure that N&Z have a right estimate.

Jan 7, 2013 at 6:47 PM | Unregistered CommenterPaul Dennis

You only need to check the rotation effect on an equatorial belt of very low thickness. A whole lot of those discs make a sphere, but rotation effect will only be the same. Indeed any average temp will be the same too, in the ideal case, no albedo, no tilt, etc. and if you calculate the temp you will get a temp, in rotation, rising to a max and falling to a min. So long as either curve is at an extreme, it will not average out. At a rotational speed sufficient to prevent any equilibrium maybe it does average out. But I suspect a general formula will not have a rate component that cancels out. And it would have to cancel out, to not appear in a formula, for the rotational rate to be irrelevant. I'll see whether this is within my remaining maths, but I suspect it is not.

Jan 7, 2013 at 10:03 PM | Registered Commenterrhoda

Ok, tried that, my missing calculus would not let me produce a working formula, but let's see what happens. I postulated a rotating disk axis perpendicular to the incoming radiation. The formula for the dark side is that on entrance to the shade there is a starting temperature and each part of the edge radiates its way down to zero or to the exit to the shade. It has a minimum temp at that point. Which is the starting point for the hot side formula, which is not the same, it starts at dark min, and receives incoming radiation at an increasing rate until 'noon' then decreasing radiation until evening. Thus the formula for the light side is different. Plotted as a curve, the two profiles are different. Depending on the temps at east and west the average is different. And the temps at east and west are determined by speed of rotation. QED? If so, proof for a disc counts for a sphere.

Jan 7, 2013 at 10:28 PM | Registered Commenterrhoda

I suspect you cannot solve the equations mathematically, only numerically, thus this thread.

Jan 7, 2013 at 10:36 PM | Unregistered CommenterTheBigYinJames


I think you are describing is the temporal profile of temperature as the disk rotates and you've correctly pointed out that this is not quite symmetric. To calculate the average temperature over a diurnal cycle from this profile will be difficult because it depends on the thermal properties of the near surface layers.

I have taken a different approach which is to determine the diurnal average insolation as a function of latitude and from this and the S-B equation the average temperature as a function of latitude. It ranges from 287K at the equator to 0K at the poles (I didn't account for the background radiation). These 'latitudinal' temperatures should be close to the constant sub surface temperatures at the same latitudes. There's not much data out there. I've seen Moon average temperatures of 250K quoted somewhere, the recent Diviner data shows an equatorial equilibrium temperature of 240K (somewhat less than my estimate BUT they report that albedo changes significantly with incidence angle which accounts for the major difference).

I integrated this function over the surface area to give an overall average temperature, which for the moon is 266K, and for Earth is 252K. I've done this two different ways now. The first, rough and ready with a spreadsheet. The second by integration (using Mathematica). I get essentially the same result using both methods.

This leads me to question what N&Z have done. I think they may have used a 'static' sphere. i.e. incident radiation on one hemisphere only. It's not immediately clear from their paper at Tallbloke's site because they give little information as to how they set up their model.

Anyhow, I'm now satisfied that the estimates of average temperature for a grey body sphere at 1AU from the sun with an albedo of 0.3 (Earth) is about 250K, and for the moon is about 266K, albeit there is a temperature gradient across the surface.

Jan 7, 2013 at 11:05 PM | Unregistered CommenterPaul Dennis


I used Mathematica to integrate tonight and the integral of cos(theta)^0.25 is

(0.8 Cos[theta]^1.25 Hypergeometric2F1[0.5, 0.625, 1.625, Cos[theta]^2] Sin[

Excuse the formatting, but I think you get the idea! Mathematica is a wonderful tool...just put the function in, give it the limits and get the numerical result without having to worry about the analytical solution!

Jan 7, 2013 at 11:10 PM | Unregistered CommenterPaul Dennis

I've had a chance to sleep on my musings of the average S-B temperature on a sphere. I'm happy with my calculation that for Earth it's about 252K and for the moon 266K. I don't know where N&Z get their result from. It's actually difficult to follow their integration since they have no model and some undefined terms and are not sure what limits they are integrating between.

These averages are very close to that estimated from the 'standard' model used to calculate the S-B temperature.

For those that still doubt these averages you need to remember that the latitudinal temperatures need to be weighted by the latitudinal area which is simply the cosine of the latitude. i.e. there is a remarkably small proportion of the sphere that is at low temperature.

Given that I believe the N&Z treatment of the average S-B temperature is wrong then all their estimates of the supposed ATE effect are wrong.

Steps in calculating the average S-B temperature (just consider a single hemisphere):

1) Determine the functional dependence of the S-B temperature on latitude using the average solar insolation at each latitude
2) Multiply this function by cos (latitude)
3) Integrate between 0 (equator) and Pi/2 (pole)

My last word on the subject.

Have a good day everyone.

Jan 8, 2013 at 8:25 AM | Unregistered CommenterPaul Dennis

Last word on the N&Z subject I can understand (although I'd like to definitely pin where they went wrong) but I'd like to develop the idea of numeric modelling a bit more, add in some more realism, just to see where the numbers go, and your help or even just observation there would be invaluable.

Jan 8, 2013 at 9:16 AM | Unregistered CommenterTheBigYinJames

Thank you Paul. I'm sure there are many, like me, that understand more as a result.

Jan 8, 2013 at 10:01 AM | Registered CommenterRichard Drake

Just got a couple of minutes, but fascinating work BigYin & Paul! How do your results stack up with the Diviner data (I expect you would need to add rotation and model the thermal conductivity of the regolith).

If RKS is looking at this he may wish to relay it to N&Z for comment?

I wonder if the definition of "average" comes into play?

Intuition tells me that rotation is an important factor when calculating an "average", just thinking about Venus (which looks the same temperature from any direction) and Mercury, which exhibits temperature differences of up to 700K (from memory).

I completely agree with BigYin about the comment on numerical modelling.

Jan 8, 2013 at 10:36 AM | Unregistered CommenterRoger Longstaff