Celestia Forums

I'm having an inordinate amount of trouble trying to replicate a graph I found in a textbook on climate, and I'm wondering if anyone here can help because this is just about the only place I can think of that I could possibly find someone who can figure this out! The problem is mostly trigonometric in nature, but is rather complex, so bear with me.

It's a graph that shows the amount of solar insolation at the top of Earth's atmosphere, in terms of time of year and latitude:

Here's what the graph looks like, just so you know what I'm talking about (he says, hoping that someone here recognises it).

Now, the formula to figure this out is rather complex:

Qo = (P/pi).S.((dm/d)^2).[(H.sin(phi).sin(delta))+(cos(phi).cos(delta).sin(H))]

where:

P = rotation period of planet in seconds
pi = 3.14159etc
S = solar constant (W/m2)
dm = mean orbital distance from sun in metres
d = current orbital distance from sun in metres (varies over a year)
H = Solar hour angle (radians)
phi = latitude of observer on planet (radians)
delta = solar declination (radians)

where:
H = arccos (-tan(phi).tan(delta))

As best as I can tell, the Solar hour angle is generally defined as "the angle measured at the celestial pole between the observer's meridian and the solar meridian". Counting from midday, it changes by 15 degrees per hour. However, in this equation H is actually the hour angle at local sunrise or sunset, where the sun intersects the observer's horizon.

(following me so far? )

Now, if I plug this into excel, with the following numbers:

P = 86400 s
S = 1360 W/m2
dm= 1.496e+11 m
d = 1.47E+11 in January to 1.52E+11 in June
delta = varies between -23 degrees in January and +23 degrees in June. (express as radians in formula)
phi = varies between 90 and -90 degrees, depending on location. (express as radians in formula)

, then I get the numbers shown on that graph. So far so good.

BUT...

there seems to be a problem with the Hour angle - when I figure it out it's symmetric about a latitude of 0 degrees. Which means that in January, as I calculate it both the north AND south poles don't see any sunlight, which of course is impossible (when one points away from the sun and is in darkness, the other pole has to be pointing toward the sun). Furthermore, the symmetry is there for all the other months too, so what I'm actually getting is that graph but with all the shaded areas reflected along the 0 degree latitude line.

I'm sure the equation for Hour Angle is correct, but something is clearly going screwy somewhere, and I'm buggered if I can spot where the problem is... but it seems to have something to do with the hour angle. I seem to be getting Hour Angles that are above 90 degrees, and I'm not convinced that this is even possible. So I think there's a trigonometry problem here.

So does anyone have any idea where I could be going wrong here? (I hope I've given enough information - if you need more I can provide it). I basically need to remove that symmetry between the poles somehow, and I'm pretty sure the formulae are all correct, so maybe my numbers are screwy? But can anyone just generally see where I'm going wrong here based on what I've described?

Are you trying to use arccos(a) with abs(a) > 1?

- Hank

Oh Evil One,

You wrote

when I figure it out it's symmetric about a latitude of 0 degrees. Which means that in January, as I calculate it both the north AND south poles don't see any sunlight,

Well, a latitude of 0 is on the equator. so it should be symmetric about that. However, I don't understand why you say "that means .. January". Delta is -23 in January, so the south pole gets light and the north pole doesn't. The poles should be (approximately) equally illuminated (or dark) at delta=0, in spring and fall.

(A quibble: why +/- 23 degrees? shouldn't it be +/- 24.5 degrees? Or is it taking into account the angular diameter of the sun, maybe?)

hank wrote:Are you trying to use arccos(a) with abs(a) > 1?

- Hank

Yes, that's what seems to be happening. What I can't figure out is why (a) is > 1 at both poles when it doesn't appear to be the case at the south pole in January from the graph.


[quote=selden]Well, a latitude of 0 is on the equator. so it should be symmetric about that. However, I don't understand why you say "that means .. January". Delta is -23 in January, so the south pole gets light and the north pole doesn't. The poles should be (approximately) equally illuminated (or dark) at delta=0, in spring and fall.[/quote]

Well yes, that's exactly the problem . That symmetry in the hour angle calculation seems to imply that in January, it's not just the North pole that is in darkness, but also the south pole (which should be illuminated, as you said and as you can see in the graph). Which makes me suspect that there's a term that I'm missing here that adds illumination to the south pole, or maybe I'm not understanding how the terms work. I mean: what, mathematically, would I have to do to get illumination going all the way to the south pole in January?

(A quibble: why +/- 23 degrees? shouldn't it be +/- 24.5 degrees? Or is it taking into account the angular diameter of the sun, maybe?)

It's because I'm just using a generally accurate number for the tilt . It doesn't make any different to the problem though.

I'll try to clean up the excel file that I've done and post that later on, maybe someone can spot the problem then.

Evil Dr Ganymede wrote:

hank wrote:Are you trying to use arccos(a) with abs(a) > 1?

- Hank

Yes, that's what seems to be happening. What I can't figure out is why (a) is > 1 at both poles when it doesn't appear to be the case at the south pole in January from the graph.

What value are you using for tan(phi) when phi=90 deg?

- Hank

Darn. I can't FTP from work. I'll have to paste something here...

So, the table below shows the calculations for January for the Hour Angle - remember that H = arccos [-tan(phi)*tan(delta)] :

The first column is Phi (Latitude) in degrees.
The second column in Phi in radians. (all angles below are in radians)
The third column is -tan(phi)
The fourth column is tan(delta)
The fifth column is -tan(phi)*tan(delta)
The sixth column is arccos[-tan(phi)*tan(delta)] (this is H)





As you can see, the value for [-tan(phi)*tan(delta)] becomes higher than 1 at both poles (above |67| degrees, to be exact), which means that I can't get an Hour angle at either pole.


Alternatively, here's the sixth column, but expressed in degrees to make it easier to visualise:




Are you even supposed to get Solar hour angles that are greater than 90 degrees? Don't forget, H is supposed to the the solar hour angle at local sunrise or sunset.... Do these numbers even make sense (what does an H angle of 90 degrees mean at the equator, given this is in January?!)

Does this help make anything clearer?

Try

H = 0 deg if abs(phi-delta) > 90 deg
H = 180 deg if abs(phi+delta) > 90 deg

(i.e. clamp the arccos argument at 1 and -1.)

- Hank

Evil Dr Ganymede wrote:Are you even supposed to get Solar hour angles that are greater than 90 degrees? Don't forget, H is supposed to the the solar hour angle at local sunrise or sunset.... Do these numbers even make sense (what does an H angle of 90 degrees mean at the equator, given this is in January?!)

I think H of 90 degrees means the sun rises and sets due east and west. H of 0 degrees means the sun never rises. H of 180 degrees means the sun never sets.

- Hank

Are you even supposed to get Solar hour angles that are greater than 90 degrees? Don't forget, H is supposed to the the solar hour angle at local sunrise or sunset.... Do these numbers even make sense (what does an H angle of 90 degrees mean at the equator, given this is in January?!)

Thinking about it, I guess it is possible - an hour angle of > 90 degrees means that there's more than half of the ellipse that the sun traces in the sky above the horizon. So presumably if you're further south than -67 degrees in January, the sun's ellipse just doesn't intersect the horizon at all, which means it's always visible. But then wouldn't that mean that the total solar insolation would be the same for those extreme southern latitudes? How could one account for that in the equation?

hank wrote:

Evil Dr Ganymede wrote:Are you even supposed to get Solar hour angles that are greater than 90 degrees? Don't forget, H is supposed to the the solar hour angle at local sunrise or sunset.... Do these numbers even make sense (what does an H angle of 90 degrees mean at the equator, given this is in January?!)

I think H of 90 degrees means the sun rises and sets due east and west. H of 0 degrees means the sun never rises. H of 180 degrees means the sun never sets.

- Hank

I think that H is the angle between the highest point that the sun reaches in the sky from a given latitude and where it intersects the horizon on the ellipse that its path traces through the sky.

If that's correct, then H=0 means that the highest that the sun gets in the sky at the horizon. So I guess it's permanently sunrise, or something. H=180 would mean that the bottom of the ellipse that the sun traces in the sky is at the horizon (ie the horizon is a tangent to one end of the ellipse), so the sun wouldn't set.

So you're right - but I think what's happening here is that there's no way of addressing what happens when the ellipse that the sun traces in the sky doesn't intersect the horizon at all. It's like asking at what angle does an object on an elliptical path intersect a line that is totally outside the ellipse - the answer of course is that it doesn't intersect it at all, so the angle is meaningless in that instance.

I've got a nice picture illustrating this that I can't upload from here

To put it another way, I think H measures (half of) that part of the sun's apparent path that is above the horizon. So if the sun never rises H = 0, and if the sun makes a complete circle above the horizon, H = 180.

- Hank

hank wrote:To put it another way, I think H measures (half of) that part of the sun's apparent path that is above the horizon. So if the sun never rises H = 0, and if the sun makes a complete circle above the horizon, H = 180.

- Hank

Ah, that would make more sense...

Here's that graphic I made, in case anyone else is having trouble figuring out what Hour angle actuallly means. The black line is the horizon as seen by the observer on the planet, the ellipses are the path that the sun seems to take in the sky.

EDG,

My understanding is as follows:

The sun's apparent path on a given day is essentially a circle centered at the celestial pole. The part of the circle which is above the horizon is determined by the sun's declination (which determines the radius of the sun's circular path) and observer's latitude (which determines the pole's elevation above the horizon). This is no different than the apparent motion of any star, except that over time the sun's declination changes. But the change is slow and practically negligible on a given day.

The observer's meridian is the great circle passing thru the celestial poles and the zenith. It intercepts the horizon due north and south. The sun crosses the meridian at local noon. The hour angle H in the formula is the angle between the meridian (the sun's position at noon) and the sun's position at sunrise/sunset, measured at the celestia pole. It is one half of the length of the part of the sun's path which is above the horizon. Thus, it indicates the proportional duration of daily insolation.

In some seasons, at sufficiently high latitudes, the sun may not rise (winter) or may not set (summer). In that case the hour angle at sunrise/sunset is undefined. However, the (half) length of the part of the sun's path which is above the horizon is still meaningful (0 or 180 degrees), and still indicates the proportional duration of daily insolation.

Your graphic incorrectly depicts the sun's path as elliptical rather than circular. Also, bear in mind that it only illustrates the changes in the sunrise/sunset hour angle on a given day due to changes in the observer's latitude (which changes the elevation of the pole). Seasonal changes in the hour angle at a given latitude due to changes in the sun's declination are caused by changes in the radius of the sun's apparent path, not the position of its center.

Does that help?

- Hank

Ok, but if I replace the ellipses with circles, then my diagram is correct, right?

And yes, the diagram shows the path the sun takes at different latitudes - but the point was just to show what the hour angle of sunset/sunrise actually was.

Anyway, I'm agreeing with you here and I'm following what you're saying, Hank . In fact, you've already solved the problem for me (thanks!). You're right, it was just a case of saying that if arccos(A) was > 1 , then the hour angle is 0 degrees ; and if arccos(A) was < -1, then the hour angle is 180 degrees. What I find annoying is that none of the texts that I referred to for this actually explicitly stated that...

EDG,

As you correctly noted in your original post, the solar hour angle is "the angle measured at the celestial pole between the observer's meridian and the solar meridian". The "observer's meridian" is the great circle passing through the celestial poles and the zenith, and the "solar meridian" is the great circle passing through the celestial poles and the sun. Thus defined the solar hour angle is always well determined. Actually, hour angle is a general term that applies to any celestial object, and is determined by the object's right ascension and the observer's local sidereal time.

Again as you correctly noted, in your formula H is supposed to be the solar hour angle at the time of local sunrise or sunset. That time depends on the sun's declination and the observer's latitude. And as you pointed out, it is undefined if the sun doesn't rise or set. So the formula fails in those circumstances.

Evidently the authors of the texts you consulted decided to leave the resolution of this detail as an exercise for the reader. (Or more likely, they never imagined that anyone would actually try to use their formula!)

My thought was that probably the formula was using this parameter to factor in the proportional duration of daily insolation, which implies what values are appropriate to substitute when the sunrise/sunset hour angles are undefined. (But that's just a guess on my part.)

There is a further problem with the formula in that tan goes to infinity at 90 deg, so the arccos argument can't be calculated near the poles. That's why I suggested using:

H = 0 deg if abs(phi-delta) > 90 deg
H = 180 deg if abs(phi+delta) > 90 deg

- Hank