Position of the sun given time of day, latitude and longitude in php

Question

This is a php implementation of Josh r code to calculate the position of the sun for a given date and time :

This is the corrected code after MvG help :

function getSunPosition($lat, $long, $year, $month, $day, $hour, $min) {
  // From https://stackoverflow.com/questions/8708048/position-of-the-sun-given-time-of-day-latitude-and-longitude?rq=1

  // Get Julian date for date at noon
  $jd = gregoriantojd($month,$day,$year);

  //correct for half-day offset
  $dayfrac = $hour / 24 - .5;

  //now set the fraction of a day      
  $frac = $dayfrac + $min / 60 / 24;
  $jd = $jd + $frac;

  // The input to the Atronomer's almanach is the difference between
  // the Julian date and JD 2451545.0 (noon, 1 January 2000)
  $time = ($jd - 2451545);
  // Ecliptic coordinates

  // Mean longitude
  $mnlong = (280.460 + 0.9856474 * $time);
  $mnlong = fmod($mnlong,360);      
  if ($mnlong < 0) $mnlong = ($mnlong + 360);

  // Mean anomaly
  $mnanom = (357.528 + 0.9856003 * $time);
  $mnanom = fmod($mnanom,360);
  if ($mnanom < 0) $mnanom = ($mnanom + 360);
  $mnanom = deg2rad($mnanom);

  // Ecliptic longitude and obliquity of ecliptic
  $eclong = ($mnlong + 1.915 * sin($mnanom) + 0.020 * sin(2 * $mnanom));
  $eclong = fmod($eclong,360);
  if ($eclong < 0) $eclong = ($eclong + 360);
  $oblqec = (23.439 - 0.0000004 * $time);
  $eclong = deg2rad($eclong);
  $oblqec = deg2rad($oblqec);

  // Celestial coordinates
  // Right ascension and declination
  $num = (cos($oblqec) * sin($eclong));
  $den = (cos($eclong));
  $ra = (atan($num / $den));
  if ($den < 0) $ra = ($ra + pi());
  if ($den >= 0 && $num <0) $ra = ($ra + 2*pi());
  $dec = (asin(sin($oblqec) * sin($eclong)));

  // Local coordinates
  // Greenwich mean sidereal time
  //$h = $hour + $min / 60 + $sec / 3600;
  $h = $hour + $min / 60;
  $gmst = (6.697375 + .0657098242 * $time + $h);
  $gmst = fmod($gmst,24);
  if ($gmst < 0) $gmst = ($gmst + 24);

  // Local mean sidereal time
  $lmst = ($gmst + $long / 15);
  $lmst = fmod($lmst,24);
  if ($lmst < 0) $lmst = ($lmst + 24);
  $lmst = deg2rad($lmst * 15);

  // Hour angle
  $ha = ($lmst - $ra);
  if ($ha < pi()) $ha = ($ha + 2*pi());
  if ($ha > pi()) $ha = ($ha - 2*pi());

  // Latitude to radians
  $lat = deg2rad($lat);

  // Azimuth and elevation
  $el = (asin(sin($dec) * sin($lat) + cos($dec) * cos($lat) * cos($ha)));
  $az = (asin(-cos($dec) * sin($ha) / cos($el)));

  // For logic and names, see Spencer, J.W. 1989. Solar Energy. 42(4):353      
  if ((sin($dec) - sin($el) * sin($lat)) >00) {
    if(sin($az) < 0) $az = ($az + 2*pi());
  } else {
    $az = (pi() - $az);
  }

  $el = rad2deg($el);
  $az = rad2deg($az);
  $lat = rad2deg($lat);

  return array(number_format($el,2),number_format($az,2));
}

This has been tested with Congo (near Equateur) lat/long : -4.77867 / 11.86364 for date Sept 1st 2013 at 10h00. In this case, the correct answer is : elevation = 67.77503 azimuth = 54.51532

Thanks for your help debuging this php code !

Greg Fabre.

It's not the only problem, but on line 52 (`$h = $hour + $min / 60 + sec / 3600;`) of the sample you reference and undefined constant `sec` - I'm guessing from the original you should either have a parameter $sec or it could be set to zero? — Clart Tent, Oct 03 '13 at 18:00
Thanks daniel, it was a typo when I copy/paste the code. I removed the seconds because it's too precise for my use. I corrected the code, and I'm still looking for the (mathematical ?) mistake ! — iero, Oct 03 '13 at 20:04
Compare all intermediate results between your code and the R code, and spot the first difference. — MvG, Oct 03 '13 at 22:49
I'll second @MvG's suggestion. In R, you can step through the code by: (1) reading in `sunPosition`'s definition; (2) doing `debugonce(sunPosition)`; and then (3) calling `sunPosition()` for a sample location. Do the same in parallel with the php equivalent of `debug` or `debugonce`, and it should take just a few minutes to find where the results diverge. Best of luck! — Josh O'Brien, Oct 04 '13 at 00:20
[here](http://ideone.com/xGlopE) is a run of [the R code by Josh O'Brien](http://stackoverflow.com/a/8764866/1468366) for your example location. The code was augmented to print all assignments along the way. I notice that the result is different from what you claim; I read elevation 67.77503 and azimuth 54.51532. Did you accidentially swap the result variables when quoting the result of the R code? — MvG, Oct 04 '13 at 07:05
[Here](http://ideone.com/eTBUh6) is a similar printing rung for your PHP code. One thing I notice is that the `$time` value already differs by almost a complete day. [Apparently](http://ideone.com/gqkFle) this is not a rounding issue. So I'd check that more carefully. By the way, [Stellarium](http://stellarium.org/) states Az=45°30'60"=45.5167, Alt=67°46'30"=67.775 for your location, 1m [AMSL](https://en.wikipedia.org/wiki/Height_above_mean_sea_level). So I'd guess the R implementation should be more correct, and the PHP date computation might be wrong in some way. — MvG, Oct 04 '13 at 07:54
Thank you for this step by step output, it helped me a lot. I didn't know this ideone.com And you're right, the dayfrac +=1 was wrong in my code ! — iero, Oct 04 '13 at 11:07

MvG · Accepted Answer · 2013-10-04T11:22:51.960

2

I believe the line

if ($dayfrac < 0) $dayfrac += 1;

is in error. If you are before noon, you don't want to refer to the same time one day later, but instead you want to specify a time before noon, i.e. subtract from the julian date which represents noon.

Removing that line, your example date corresponds to the one computed using http://www.imcce.fr/en/grandpublic/temps/jour_julien.php, namely 2456536.9166666665. The resulting

$el = 67.775028608168
$az = 54.515316112281

looks pretty good to me. In particular, it agrees with the R run

elevation = 67.77503
azimuth = 54.51532

and also with what Stellarium says (although I quoted this incorrectly in a comment above):

Alt = 67°46'30" = 67.775
Az  = 54°30'60" = 45.5167

It also (almost) agrees with sunearthtools.com, so I guess you made a mistake when first entering the data there:

Screenshot from sunearthtools

So I'd say that solves the problem.

edited Oct 04 '13 at 11:22

answered Oct 04 '13 at 08:14

MvG

57,380
22
148
276

Thanks, you are right. The answer I expected was from http://www.sunearthtools.com and it seems to be buggy (maybe near equateur). Thank you very much for your help, I will correct the code it can help other people ! – iero Oct 04 '13 at 11:05
@iero: Can't reproduce your problem with sunearthtools.com; works for me as the screenshot above demonstrates. – MvG Oct 04 '13 at 11:23
I cannot upload image as i need '10 reputations'. But you can see the results I get on [this link](http://oi39.tinypic.com/ff1ymw.jpg) – iero Oct 04 '13 at 15:27
@iero: Found out how to reproduce this: You had DST checked, so the time is 10:00 in UTC+1, therefore 9:00 UTC. To get 10:00 UTC you either have to uncheck that box or enter 11:00. – MvG Oct 04 '13 at 19:14
Nice catch ! I didn't know what was DST but as I checked "GMT 0", I didn't expect to not be in UTC time ! Thanks again for all your help, I hope that this php code will help somebody else ! Have a nice (and sunny) day ! – iero Oct 05 '13 at 08:42

Position of the sun given time of day, latitude and longitude in php

1 Answers1