2
$\begingroup$

As the title says I am trying to calculate solar coordinates of the sun for a given location and time. I have spent many 10's of hours on this without success so I will sincerely appreciate help on this. My RA value seems to be right on but the Dec, Alt and Az are off.

Alt value is within 6 degrees but I cannot figure out why it is that close because the latitude is in radians instead of degrees like the other values. If I change latitude to degrees the values don't even make sense.

I have included all the values of the intermediate variables to help see at a glance what might be off. I hard coded the JD and the Local Hour Angle just to take those out of the equation for debugging. Links are shown as to where I got those values.

The C# code is below although the programming language doesn't matter. The formulas are the same in any language.

  //for formulas
            double D = 2458599.125000 - 2451545.0; //the  first  term is the jd of interest; 2nd term is jd of 1-1-2000 at 12:00:00.0
            double Latitude = 41.369341;
            double Longitude = 81.885493;

            double gTmp = 357.529 + 0.98560028 * D;  //mean anomaly of the sun  = 7310.0765751550007
            double gTmp1 = (gTmp % 360);     //subtract multiples of 360 (returns remainder in degrees)  =  110.07657515500068
            double g = gTmp1 * Math.PI / 180;   //Convert to radians = 1.9211986657737494
            double qTmp = 280.459 + 0.98564736 * D; //mean longitude of the sun = 7233.33868336
            double qTmp1 = (qTmp % 360); //subtract multiples of 360 (returns remainder in degrees) =  33.338683359999777
            double q = qTmp1;   //don't convert to radians here; convert after terms in degrees are added in the computation of L; =33.338683359999777
            double L_tmp = q + 1.915 * Math.Sin(g) + 0.020 * Math.Sin(2 * g);//geocentric apparent ecliptic longitude of the sun(adjusted for aberration) = 35.124421105295383
            var L = (L_tmp % 360) * Math.PI / 180;   //Convert to radians; //subtract multiples of 360 = 0.61303679614439033
            double eTmp = 23.4393 - 0.00000036 * D;  // = 23.436760515
            double e = eTmp * Math.PI / 180;   //Convert to radians;   = 0.40904863698815186

            //http://www.neoprogrammics.com/de405_usno_ae98/DE405_Sun.php
            // formulas per https://aa.usno.navy.mil/faq/docs/SunApprox.php
            double RA = (Math.Atan2(Math.Cos(e) * Math.Sin(L), Math.Cos(L))) * 180 / Math.PI; //right ascension converted from radians to degrees = 32.838762087228496
            double dec = (Math.Asin(Math.Sin(e) * Math.Sin(L))) * 180 / Math.PI;  //declination converted from radians to degrees =  3.228747909189217         
            string[] sunCoordinates = new string[2];   
            sunCoordinates[0] = ToHMSStr(RA);    // = "2h 11m 21.3s"
            sunCoordinates[1] = ToDMSstr(dec);  //  = "13d 13m 43.49s"

            //compute alt and az per https://aa.usno.navy.mil/faq/docs/Alt_Az.php and https://aa.usno.navy.mil/faq/docs/GAST.php
            //LHA for testing is from https://ssd.jpl.nasa.gov/horizons.cgi#results and https://www.vercalendario.info/en/how/convert-ra-degrees-hours.html
            //from http://www.neoprogrammics.com/de405_usno_ae98/DE405_Sun.php  
            double LHA = 307.39;   //Local Hour Angle        
            double lat = Latitude * Math.PI / 180; //Convert to radians; //Latitude is defined above; We get sign (below) after converted to Rad????  = 0.72203120983028346
            double alt = (Math.Cos(LHA) * Math.Cos(dec) * Math.Cos(lat) + Math.Sin(dec) * Math.Sin(lat)) * 180 / Math.PI; //convert from radians  = 53.266598643452909
            double az =Math.Atan2(-Math.Sin(LHA) * Math.PI / 180, (Math.Tan(dec)*Math.Cos(lat) * Math.PI / 180 - Math.Sin(lat) * Math.Cos(LHA)) * Math.PI / 180); //convert each term to radians   = 2.4585704411338072
            az = az * 180 / Math.PI;   //convert from radians   = 140.86570991258415

```
$\endgroup$
7
  • $\begingroup$ What program language are you using for this calculation? $\endgroup$
    – Fil
    May 15, 2019 at 19:07
  • $\begingroup$ If some resources in this answer or especially in this answer are helpful, feel free to go ahead and write a short answer to your own question. $\endgroup$
    – uhoh
    May 16, 2019 at 6:54
  • $\begingroup$ Is valuable. Actually I already was aware of the C sharp implementation of the sofa Library. It is meticulously crafted but even the guy who wrote it said he has not had time to actually use it. So the implementation documentation is a bit sketchy. $\endgroup$ May 18, 2019 at 14:24
  • $\begingroup$ The spreadsheet for which included a link, is a work of art. All the formulas are easily visible and can the implemented in your programming language of choice. Do you know if there's such a spreadsheet for lunar ephemeris.? $\endgroup$ May 18, 2019 at 14:33
  • $\begingroup$ Are you sure your longitude is positive? That puts the location in the western part of China. Nothing wrong with that, but negative that longitude puts you in NE Ohio. $\endgroup$
    – user21
    May 18, 2019 at 18:07

3 Answers 3

1
$\begingroup$

For the Sun at JD 2458599.125, RA = 32.8° and dec = 13.2° are correct regardless of geographic location. In the alt/az calculation,

  • LHA should be 323.6° and converted to radians, and
  • the expression for alt needs a call to Math.Asin.

Trigonometric expressions are easier to verify if they are free of inline unit conversions. Store angles in radians where possible; otherwise make it obvious with a name suffix such as _deg.

When computing the local hour angle, you can approximate GAST with GMST. The difference eqeq is smaller than the approximation error of the Sun position algorithm you're using, and the alternative GMST formula is valid for D within a few centuries of J2000.0.

$\endgroup$
0
$\begingroup$

I recommend you to use astropy, it's a python package with plenty of astronomical libraries.
Check this tutorial on Determining and plotting the altitude/azimuth of a celestial object.
You have the specific function

from astropy.coordinates import get_sun

Hope it helps.

$\endgroup$
2
  • $\begingroup$ Thanks but this will end up being used in a mobile app so C# is the language of choice. Do you think I would have any luck using Python just to get the formulas and implement them myself in C#. All I need are the formulas and an explanation of when and where to convert them to and from radians etc. I have never used Python but this would be pretty basic. $\endgroup$ May 15, 2019 at 21:39
  • $\begingroup$ @noCodeMonkeys for Python use Skyfield or the (now deprecated) PyEphem. Skyfield uses an ephemeris and it needs to download information regularly (a bit table to interpolate) but PyEphem uses something a little different. That will be fun on a computer, but not for a mobile app. $\endgroup$
    – uhoh
    May 16, 2019 at 6:52
0
$\begingroup$

For anyone looking for the same information, I would point you to an excellent paper by Zhang, Stackhouse, Macpherson, and Mikovitz where they do precisely what you are looking for, providing equations along with a code implementation, albeit in Fortran. I implemented it myself in TypeScript. Do keep in mind that the term $$T_{GMT}$$ is fractional days since UTC 00:00 on the day in question.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .