Jeff, thanks for that. => Would you know what the accuracy of the SunPosition Calculations would or should be?

Here's the US Naval Laboratory for another comparison http://aa.usno.navy.mil/ they only have per minute on the site, but the numbers agree with the NASA 17:13 20.2 279.0 17:14 20.0 279.2 17:15 19.8 279.4 17:16 19.6 279.5

As Sjoerd kindly pointed out there's Mathematica did not take Atmospheric Refraction into consideration which will be something like 10 arc minutes at 5 degrees of altitude 1 at 45 and 35 at the horizon, or at 20 degrees of altitude something like 3 arc minutes (using Bennett's 1982 calculation) which is about 0.05 degrees. However this should increase the perceived altitude. So for the SunPosition Calculation this correction is big enough in magnitude, nor is it in the right direction to help with the reconciliation, as the SunPosition is supposed to exclude it.

=> Any thoughts how I reconcile these differences or am I missing something obvious?

=> Would you know what is supposed to be the most accurate estimate (I understand some approximation would be made for atmospheric refraction)?

Are you sure the output is ok? I've refreshed with RebuildPacletData[]

If http : // ssd.jpl.nasa.gov/horizons.cgi?s_loc = 1 #top is the gold standard Then for the sample point I looked at {43.953074,4.472485} and at 2014-07-22 17:15:00 UTC AstronomicalData["Sun",.... does not agree with SunPosition[GeoPosition and both are different from the nasa calculation e.g. for the altitude for this point respectively I get 19.85892, 20.06311334010335, and 19.7609

Unless I've done some thing stupid I can't see that 0.3 degrees out is much of a success

copied from mathematica 10 notebook

In[26]:= {AstronomicalData["Sun",{"Azimuth",{2014,7,22,18,15,0},{43.953074,4.472485}},TimeZone->1], AstronomicalData["Sun",{"Altitude",{2014,7,22,18,15,0},{43.953074,4.472485}},TimeZone->1]} Out[26]= {279.28408[Degree],19.85892[Degree]} In[27]:= SunPosition[GeoPosition[{43.953074,4.472485}], DateObject[{2014,07,22,19,15, 0}, TimeZone -> IntegerPart[LocalTimeZone[Entity["City",{"Paris","IleDeFrance","France"}], DateObject[{2014,07,22,19,15, 0},TimeZone-> 0]]]] ]//FullForm Out[27]//FullForm= List[Quantity[279.096879879350734008354638938468698465083.6876323489782243,"AngularDegrees"],Quantity[20.063113340103351466195182502734913517392.5442756941572635,"AngularDegrees"]]

Result of values in the jpl.nasa site 279.3737 19.7609

I suppose I'm still not clear why the two different methods in Mathematica (AstronomicalData and SunPostion) don't agree. If I look at the Azimuth for which refraction presumably doesn't make much difference then I get

AstronomicalData 279.28408
SunPosition 279.09687
and NASA 279.3737

which there's roughly 0.1 degrees between AstronomicalData and NASA and 0.2 degrees between AstronomicalData and SunPosition

=> What would you consider to be the gold standard?

I'm interested as I was wondering how accurate one could plausibly make a sundial for for example the next 100 years and thought it would be fairly easy to set up on Mathematica, but am now not sure. the timezone correction is 4 minutes a degree, i.e. 0.1 degrees is 24 seconds 0.1 degrees is something like 3.5 mm at 2 metres

I spoke too soon. It hasn't pushed to the production server just yet. Should hopefully be on production tomorrow. Sorry if I caused additional confusion.

I have found that, in addition to SunPosition returning the wrong result, DaylightQ and SiderealTime also produce the wrong result. SiderealTime code is almost certainly called by the SunPosition code, so this might be a good place to start looking.

Users here get the same results I do. I'm investigating for any issues I can think of that might cause such behavior.

I checked by setting my time zone to Champaign's

Block[{$TimeZone = -5}, 
  SunPosition[GeoPosition[{52.37`, 4.89`}], DateObject[{2013, 3, 1, #, 0, 0}, TimeZone -> 1]] & /@   Range[0, 23]
]

I now get your results.

So, you get the correct results because part of the calculation somehow assumes your time zone, which yields incorrect results for anyone not in your time zone. You are among the few who can't see the incorrect results ;-)

I saw your analysis. The API call includes a fair amount of information that may or may not be needed for a specific input. I see nothing obviously wrong with the API arguments being sent assuming for you:

$TimeZone -> 2
$GeoLocation -> GeoPosition[{52.09, 5.12}]

You only gave a single API call, but did not specify which value of date was used as input for that call.

Your input date specifies that you are trying to feed in a different timezone (1) from where you are (2). TimeZone conversions need this information to be correct or you could get confusing results.

Concerning SunPosition returning different results from AstronomicalData, I am unable to reproduce this in any significant way. Using the exact input you used:

SunPosition[GeoPosition[{52.37, 4.89}], DateObject[{2013, 3, 1, #, 0, 0}, TimeZone -> 1]] & /@ Range[0, 23]

I get the same results as with AstronomicalData (possible differences in precision and also results coming back as Quantities ignored.)

I definitely am unable to get any significant difference such as what you showed in your original post.

If you are using a pre-release version of Mathematica 10 or something related, that might be to blame.

In[1]:= SunPosition[GeoPosition[{52.37, 4.89}], DateObject[{2013, 3, 1, #, 0, 0}, TimeZone -> 1]] & /@ Range[0, 23]

Out[1]= {{Quantity[341.1, "AngularDegrees"], Quantity[-43.9, "AngularDegrees"]}, {Quantity[1.9, "AngularDegrees"], Quantity[-45.2, "AngularDegrees"]}, {Quantity[22.6, "AngularDegrees"], Quantity[-43.3, "AngularDegrees"]}, {Quantity[41.2, "AngularDegrees"], Quantity[-38.4, "AngularDegrees"]}, {Quantity[57.3, "AngularDegrees"], Quantity[-31.4, "AngularDegrees"]}, {Quantity[71.2, "AngularDegrees"], Quantity[-23.2, "AngularDegrees"]}, {Quantity[83.8, "AngularDegrees"], Quantity[-14.3, "AngularDegrees"]}, {Quantity[95.7, "AngularDegrees"], Quantity[-5.1, "AngularDegrees"]}, {Quantity[107.6, "AngularDegrees"], Quantity[3.9, "AngularDegrees"]}, {Quantity[120.0, "AngularDegrees"], Quantity[12.3, "AngularDegrees"]}, {Quantity[133.5, "AngularDegrees"], Quantity[19.6, "AngularDegrees"]}, {Quantity[148.3, "AngularDegrees"], Quantity[25.4, "AngularDegrees"]}, {Quantity[164.5, "AngularDegrees"], Quantity[29.1, "AngularDegrees"]}, {Quantity[181.6, "AngularDegrees"], Quantity[30.2, "AngularDegrees"]}, {Quantity[198.6, "AngularDegrees"], Quantity[28.6, "AngularDegrees"]}, {Quantity[214.6, "AngularDegrees"], Quantity[24.5, "AngularDegrees"]}, {Quantity[229.2, "AngularDegrees"], Quantity[18.4, "AngularDegrees"]}, {Quantity[242.5, "AngularDegrees"], Quantity[10.9, "AngularDegrees"]}, {Quantity[254.8, "AngularDegrees"], Quantity[2.4, "AngularDegrees"]}, {Quantity[266.6, "AngularDegrees"], Quantity[-6.7, "AngularDegrees"]}, {Quantity[278.6, "AngularDegrees"], Quantity[-15.8, "AngularDegrees"]}, {Quantity[291.4, "AngularDegrees"], Quantity[-24.6, "AngularDegrees"]}, {Quantity[305.7, "AngularDegrees"], Quantity[-32.6, "AngularDegrees"]}, {Quantity[322.2, "AngularDegrees"], Quantity[-39.2, "AngularDegrees"]}}

The speed overhead is a known issue that we are actively planning to address in a future release.

Very good posts, Sjoerd. For my projects, I've started doing the calculations in other languages, using well-known datasets and algorithms, and then finally feeding the various chunks of calculated data to Mathematica for higher-level assembly and presentation. It's not the most desirable solution, but as Mathematica becomes a higher and higher level language, it's a pragmatic approach, at the moment. Of course, the temptation is to try and do it in just one environment - but then, which one?! :)

For one-off questions or exploring a new problem space, firing off a few dozen WolframAlpha queries is fine, but, like elsewhere in computing, relying on the cloud isn't always going to be the ideal solution. There are some obvious disadvantages, too: you don't know where it comes from or how it's been calculated, how reliable it is, or when or whether it's been updated or modified, missing and unexplained data points can't be tracked down and repaired, and so on. Unless the cloud service is quick and very reliable, and well documented, it's easy to end up with problems that aren't easily fixed.

I can't answer the other issues since I do not use these functions. But on this point you raise

It needs an Internet connection

How would propose that Mathematica get the current ephemeris data for solar system planets and the sun without using an updated database somewhere on the internet?

I use JPL Horizon website myself for this information. You might want to see if there is a way to use this site via an API call to obtain this information. It is free.

I converted to Mathematica a VBS macro that you can find in this Excel sheet. It is written by Greg Pelletier and based on a NOAA Javascript program which was used to drive their sun calculation page. I performed a direct and ugly conversion without any optimization whatsoever. Testing it with the method above I found that it is about 30 times faster than V9 AstronomicalData and about 1100 times faster than V10's SunPosition.

Testing the output and taking the NASA website you refer to as the golden standard I get for the sun positions starting at 1 am (using AstronomicalData for V10):

Azimuth:

Elevation:

I declare winner the one with the smallest absolute deviation from the NASA data.

As to the need for ephemeris: The following plot compares the results for Pelletier's and V10's AstronomicalData in the year 3013, thousand years in the future (CircleDot is Mathematica):