在PHP中根据时间、纬度和经度确定太阳的位置

Question

在PHP中根据时间、纬度和经度确定太阳的位置

4

这是一个以PHP实现的Josh r代码，用于计算给定日期和时间下太阳的位置：

这是经过MvG帮助之后做出的修正后的代码：

function getSunPosition($lat, $long, $year, $month, $day, $hour, $min) {
  // From https://dev59.com/lmoy5IYBdhLWcg3wMLSj

  // 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));
}

这已经在2013年9月1日10:00针对刚果（赤道附近）的纬度/经度：-4.77867 / 11.86364进行了测试。在此情况下，正确答案为：海拔高度 = 67.77503 方位角 = 54.51532

感谢您帮助调试这段PHP代码！

Greg Fabre.

- iero

这不是唯一的问题，但在您引用的示例的第52行（$h = $hour + $min / 60 + sec / 3600;）中存在未定义常量sec - 我猜原始代码应该要么有一个参数 $sec，要么可以将其设置为零？ - Clart Tent

谢谢Daniel，我复制/粘贴代码时打错了一个字。我删除了秒，因为它对我的使用来说太精确了。我已经纠正了代码，但仍在寻找（数学上的？）错误！ - iero

1

比较你的代码和 R 代码之间的所有中间结果，并找出第一个不同之处。 - MvG

我赞同@MvG的建议。在R中，您可以通过以下步骤逐步执行代码：（1）阅读“sunPosition”的定义；（2）执行“debugonce(sunPosition)”；然后（3）为样本位置调用“sunPosition()”。在与“debug”或“debugonce”相等的php版本上进行相同操作，应该只需要几分钟即可找到结果分歧的地方。祝好运！ - Josh O'Brien

这是Josh O'Brien的R代码运行结果，用于你提供的位置示例。该代码被增强以打印所有赋值过程。我注意到结果与你所说的不同；我读取高度为67.77503，方位角为54.51532。你在引用R代码结果时是否意外交换了结果变量？ - MvG

显示剩余2条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- MvG · Accepted Answer

我相信这句话的含义是

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

出现错误。如果当前时间在中午之前，你不应该指的是一天后同样的时间，而是要指定一个在中午之前的时间，也就是从代表中午的儒略日中减去相应的时间。

删除那行代码后，你的示例日期与使用http://www.imcce.fr/en/grandpublic/temps/jour_julien.php计算得到的日期相符，即2456536.9166666665。结果如下：

$el = 67.775028608168
$az = 54.515316112281

我认为这看起来相当不错。特别是，它与R运行结果一致。

elevation = 67.77503
azimuth = 54.51532

并且与Stellarium的说法相符（尽管我在上面的评论中引用有误）：

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

它（几乎）与sunearthtools.com的结果相符，因此我猜您在首次输入数据时出错了：sunearthtools.com。

因此，我认为这解决了这个问题。

备注：已保留HTML标签。