Where is my calculation of sun angle incorrect?

[AWoods]

Hi there, I hope this is the right place for this question and furthermore doesn't appear rude that my first post is a question! If so, I do apologise, I haven't been able to find an appropriate place to ask this.

I'm working on calculating the position of the sun at a given day of the year for a given location. I have all the calculations laid out correctly I believe, but my elevation and azimuth are coming back incorrect.
I think this is due to me not using radians/degrees correctly, but honestly, I've dived right into the deep end given that I haven't studied maths and physics in maybe 5 years.

Just to clarify given the guidelines I just read, this isn't homework, it's a personal project I'm working on for simulating a dynamic time of day system for games, if that still classifies, sorry and feel free to remove my thread.

Any help you could give, or even just a nudge in the right direction would be fantastic, I'll post my current calculations below:

These are taken directly from my code, I'll try to simplify it where possible in the hopes that it's simple enough to read. I can provide further clarification if needed.

Given inputs are: current time decimalTimeHourMin(format: hour.minute), date, GMTOffset (in hours), latitude and longitude.
Acos and Asin signify Arccos and Arcsine inverse functions.

localStandardTimeMeridian = 15 * GMTOffset;

fractionalYear = (360 / 365) * ( (date.DayOfYear - 1) - 81);

*** This is where is starts to get harder to read admittedly, the Deg2Rad bits equate to (PI * 2) / 360 ***

equationOfTime = (9.87f * Sin( (2 * fractionalYear) * Deg2Rad) ) - (7.53f * Cos(fractionalYear * Deg2Rad) ) - (1.5f * Sin(fractionalYear * Deg2Rad) );

timeCorrection = 4*(longitude - localStandardTimeMeridian) + equationOfTime;

localSolarTime = decimalTimeHourMin + (timeCorrection / 60);

hourAngle = 15 * (localSolarTime - 12);

declination = 23.45f * Sin(fractionalYear * Deg2Rad);

*** These are the 2 calculations that do not work out, all other calculations above are correct to an acceptable degree of accuracy ***

elevation = Asin( (Sin(declination) * Sin(lattitude) ) + (Cos(declination) * Cos(lattitude) * Cos(hourAngle) ) );

azimuth = Acos( (Sin(declination) * Sin(lattitude) ) - (Cos(declination) * Cos(lattitude) * Cos(hourAngle) ) / Cos(elevation) );

I'm fairly positive that I need to handle some of these factors as radians, though I cannot work out where.
Have I done something terribly wrong?
I also have an Excel Spreadsheet with appropriate formulas in, which does the same thing, I can post a link to that if it would be helpful.

Thanks for reading

 

[Filip Larsen]

I'm puzzled.

Assuming everything above the two last equations is correct (I haven't checked) it should be a straight forward thing to check if each argument to a trigonometric function is indeed in radians, and if not convert it. You reference three variables in your last two equations, so perhaps you should start check if those three variables are in radians or not?

Later: You may also want to check if the last two equations are correct (e.g., compare with relevant equations on [1]), and, if your environment supports it as most do, you may also want to use atan2 [2] instead of atan.

[1] http://en.wikipedia.org/wiki/Celestial_coordinate_system
[2] http://en.wikipedia.org/wiki/Atan2

 

[AWoods]

Hi there, thanks for responding.

So I should have probably included this in my first post, but this is a calculator for the calculations I'm trying to replicate:
http://www.pveducation.org/pvcdrom/properties-of-sunlight/sun-position-calculator

I've double checked the equations and believe they are correct. (for reference, they are at the bottom of this page: http://www.pveducation.org/pvcdrom/properties-of-sunlight/suns-position )

I'm entering these values:
Date: 1st Jan
Time: 13:00
Lat: 52
Long: -1
GMT Offset: 1

As the results I'm looking for are the altitude and azimuth, the results should be:
Altitude: 14.97
Azimuth: 178.16

However in my calculations I get:
Altitude (elevation): ~0.9315
Azimuth: ~0.5125

I have tried combinations of radian conversions on the variables, but cannot get close to the correct result.
I will look into other sources for the equations, it's possible that this is not the correct way to get the results I'm looking for.

Thanks for your help

 

[Filip Larsen]

The last equation on that calculate page is not correct. I made a quick spreadsheet and got the same results as the calculator up until that equation. If I use a the correct equation (like that on [1]) I get same result as the calculator, that is, 178.16 deg.

Make sure that every argument you pass to sin or cos is converted to radians if it was in degrees. In my sheet I note the unit after the number and make a new equation converting deg to radians so there is no doubt. Also remember that the result from the arcus functions (asin, acos) are in radians and needs to be converted to degree if you want them in that unit.

[1] http://www.stargazing.net/kepler/altaz.html

 

[willem2]

One problem is that the results of acos and asin in the last two equations are in radians, and you should convert them to degrees.
In the equation for latitude, you should take the arcsine of the entire expression, and not just of Sin(declination) * Sin(latitude). Watch which closing parenthesis belongs to asin.
More or less the same thing happens in the second expression. However, as mentioned by Filip Larsen, the last equation on that page is wrong, it does not give the same answer as the solar position calculator, also on the pveducation site, which seems to be ok.

 

[D H]

The last equation on that calculate page is not correct. I made a quick spreadsheet and got the same results as the calculator up until that equation.

I made a quick spreadsheet and got the same results as the calculator, including that equation. Perhaps you made a mistake entering that last equation. The denominator in the last equation is $\cos(\alpha)$, where $\alpha$ is the elevation (or altitude, hence the use of $\alpha$).

That said, the equation as printed on that page is not quite correct. It's only valid for negative hour angles (i.e., morning). The previous page on the calculations page gives a fully correct formula. If the hour angle is positive, you need to subtract the value calculated using inverse cosine from 360° (i.e., use 360°-acos(...) instead of acos(...)).

@AWoods : All of the angles on that page and related pages are in degrees. Your calculator can easily switch between degrees and radians. Your spreadsheet provides an easy but more verbose mechanism. Use the functions radians to convert degrees to radians, degrees to convert from radians to degrees.

 

[AWoods]

Firstly thanks for all the responses and time you have put into helping me.

I haven't managed to get a correct result yet but just want to try and break this down a bit further so I can try and work out what should be happening.
I'd like to focus on the elevation calculation as I think I should be able to get azimuth fairly easily after that point.

So, the parameters for the elevation are the declination, latitude and hour angle, all of which are in degrees.
For this reason I've converted them all to radians in the equation for elevation.

I'm still getting incorrect results from this in my spreadsheet so far. But I feel like I'm close.
I've split the 2 'sides' of the equation and get odd results in both.

(including converting to radians)
Sin(declination) * Sin(lattitude) = -0.308
Cos(declination) * Cos(lattitude) * Cos(hourAngle) = 0.566

The whole equation (with the Asin included) came out to 0.261 (ish).

I'll include my whole excel equation for the elevation just for completeness:
=ASIN(((SIN(RADIANS(B13))*SIN(RADIANS(B4)))+(COS(RADIANS(B13))*COS(RADIANS(B4))*COS(RADIANS(B12)))))

I can't work out quite why the values I'm getting are so low.

 

[Filip Larsen]

I made a quick spreadsheet and got the same results as the calculator, including that equation. Perhaps you made a mistake entering that last equation.

Thanks for spotting that, DH. You are, as usual, spot on.

Now I try again I get same result using with both forms of equations for altitude. Even though I typed it in and verified it several times I must have made the same mistake each time. As you note, it is very likely that I got the use of $\alpha$ wrong as the use of this symbol was the main reason I ended up concluding the equation had typo, which was further confirmed when a different equation for elevation yielded the expected result. My bad.

 

[Filip Larsen]

I can't work out quite why the values I'm getting are so low.

Note, that the values from Asin is in radians. If you convert that back to degrees you get 14.97 as expected. I get same values as you do.

 

[AWoods]

Note, that the values from Asin is in radians. If you convert that back to degrees you get 14.97 as expected. I get same values as you do.

That was exactly it!
I feel daft now, thank you so much.
That does make a lot of sense though, it's just difficult to keep track of what values are in what format.

Thanks everyone for the help, I'll maybe post some pictures if I get the whole thing working as intended, no idea how long that will take though.

 

[AWoods]

I have one small further question to add.

I have an issue where the azimuth jumps around 50 degrees when the solar time reaches 12.
This is due to having implemented the fix to the equation suggested earlier, Azimuth = 360° - Azi, for LST > 12 or HRA >0.
The azimuth at the point of the switch is around 140 degrees (this does vary depending on day of the year).

My question is, is that azimuth calculation entirely correct?
It strikes me that in order for it to be correct, azimuth must be 180 or 0 at the point solar time reaches 12?

If that is correct, could it be that I'm getting incorrect values for local solar time?
Because of the time correction applied to it, there is of course an offset from the normal 0-24 range that takes it into the negatives at times.
I'm not sure if that has anything to do with this issue, as there are a lot of factors at work in these calculations.

If nothing else, does anyone know where I can find more information about this kind of application of azimuth calculations, the Wikipedia article and other bits I've found so far haven't helped much.

 

[Filip Larsen]

When using the equations in the form presented on the pveducation.org site it sure seems you should be able to use both criteria for determining if the Sun is east or west of the meridian. I use HRA > 0 and, for the few samples tested, I get same result as the calculator on the site.

You should get a "raw" azimuth at 180 from Degrees(Asin(-1)) when LST = 12 or HRA = 0. You may want to check that you get that.

 

[AWoods]

It turns out it was me being daft again, and was to do with an incorrect azimuth calculation, the solar morning/night switch does work correctly.
I was using sine for latitude where I should have been using cosine, it appears to all work as expected, and quite smoothly too.

Error is fairly low, which is great.
The cost of performing these calculations is incredibly low, which is surprising given that sine and cosine tend to be expensive for real-time application.

Now I can move on to more visual changes, like intensity and colour of the sun based on time of day.

Thanks for the help!

 

[tadchem]

I could not help but notice that you are using the term "latitude' as an input, but the equations are citing "lattitude" (with a double-t). Are the equations with the problematic results really using "lattitude" as an argument in lieu of "latitude"?

 

[AWoods]

I could not help but notice that you are using the term "latitude' as an input, but the equations are citing "lattitude" (with a double-t). Are the equations with the problematic results really using "lattitude" as an argument in lieu of "latitude"?

Ah yes, that's just bad spelling on my part, I was using the incorrect spelling but have now corrected all instances to "latitude".

 

[Shahid Naeem]

I just wanted to thank you all for this conversation. I've been struggling with a model I was developing for estimating leaf temperature using the models developed in the photovoltaics business for estimating the radiant energy incident on a panel of any angle, azimuth, at any location on Earth, at any time - but kept getting poor results. Out of curiosity, I wondered if anyone else was having similar problems and, much to my amazement, you folks have described and very nicely considered the likely problems leading to erroneous calculations. So, again, many thanks to you all for this public exchange.