#!/usr/bin/perl -w
#
# Solar elevation calculator v1.4
# by The Almighty Pegasus Epsilon <pegasus@pimpninjas.org>
# (C) 2010-2013 All rights reserved
#
# Outputs a message like /usr/games/pom, except about the sun's location.
#
# CHANGELOG
# v1.0 - ?? 2010 - First working version.
# v1.1 - Sep 13 2011 - Fixed (unused, commented) Solar Azimuth calculation.
# v1.2 - Sep 13 2011 - Updated to allow selection of home planet.
# v1.3 - Mar 14 2012 - Added command-line parsing (still no time setting yet).
# v1.4 - Dec 10 2013 - Fixed Sol's Angular Radius.
# v1.5 - Jan 22 2021 - POSIX module no longer exports isdigit()
### CONFIG ###
# local degrees longitude west
#my $longitude_west = 94.5471387;
my $longitude_west = 80.7141434;
# local degrees latitude north (phi)
#my $phi = 38.917505;
my $phi = 28.0270573;
# display precision
my $precision = 4;
# default planet
my $planet = "earth";
### END CONFIG ###
use strict;
use Math::Trig;
use Time::Local;
use POSIX qw(fmod);
# Days since Jan 01 2000 (the Epoch) -- high precision
my $day = ((time - timegm(0,0,0,1,0,100)) / 86400) - .5;
#print "current julian day (J): ".($day + 2451545)."\n";
# replace with [[[[[cc]yy]mm]dd]HH] ala pom
sub usage () {
print "Usage: $0 [planet] [-l lat lon] [[[[[cc]yy]mm]dd]HH]\n";
exit
}
# parse command line
while ($_ = shift @ARGV) {
s/^(-l)(\d+)/$1/ && unshift @ARGV, $2;
/^-l$/ && do {
usage unless defined ($phi = shift);
usage unless defined ($longitude_west = shift);
usage unless /\d/ =~ $phi and /\d/ =~ $longitude_west;
$phi %= 360; $phi -= 360 if $phi > 180; $phi %= 180;
$longitude_west %= 360;
$longitude_west -= 360 if $longitude_west > 180;
$longitude_west %= 180;
next;
};
if (/\d/ !~ $_) { $planet = lc $_ }
else { die "setting the time is not supported yet.\n" }
}
#print "local latitude north (phi): $phi degrees\n";
#print "local longitude west: $longitude_west degrees\n";
#print "on planet: $planet\n";
# degrees to radians
sub d2r ($) { $_[0] * pi / 180; }
# radians to degrees
sub r2d ($) { $_[0] * 180 / pi }
# pluto is a planet, goddamnit.
my %planet = (
# All planetary data is taken from the Epoch Jan 01 2000.
# "anomaly" below is data to solve for the current mean anomaly for
# each planet. More info: http://en.wikipedia.org/wiki/Mean_anomaly
# "angle" is the angle of the orbit of the given planet. Zero degrees
# is perihelion, closest possible position to the solar system
# barycenter. More info: http://en.wikipedia.org/wiki/Periapsis
# http://en.wikipedia.org/wiki/Center_of_mass
# "precession" is the amount of change in the eccentricity angle per
# earth day. Note that for earth it is roughly equivalent to <year
# in days> / 360. The difference between <year in days> / 360 and
# values below accounts for precession of the orbital path itself.
# More info: http://en.wikipedia.org/wiki/Perihelion_precession
# "center" is data to solve the Equation of the center to find the
# true anomaly from the mean anomaly in accordance with the Equation
# of Kepler. The latter is unsolvable, though (thanks, pi!), so
# we'll just get decently close. More info:
# http://en.wikipedia.org/wiki/Gravitational_two-body_problem
# "ecliptic" is data to solve for the ecliptical coordinates of the
# sun from the planet's perspective (planetocentric).
# http://en.wikipedia.org/wiki/Ecliptic
# "perihelion" is the ecliptic longitude of the sun.
# "obliquity" http://en.wikipedia.org/wiki/Obliquity
# to be continued...
mercury => {
anomaly => {
angle => 174.7948,
precession => 4.09233445
},
center => [ 23.4400, 2.9818, 0.5255, 0.1058, 0.0241, 0.0055 ],
ecliptic => {
perihelion => 111.5943,
obliquity => 0.02
},
ascension => [ 0 ],
declination => [ 0 ],
sidereal => [ 13.5964, 6.1385025 ]
},
venus => {
anomaly => {
angle => 50.4161,
precession => 1.60213034
},
center => [ 0.7758, 0.0033 ],
ecliptic => {
perihelion => 73.9519,
obliquity => 2.64
},
ascension => [ -0.0305 ],
declination => [ 2.6427, 0.0009 ],
sidereal => [ 215.2995, -1.4813688 ]
},
earth => {
anomaly => {
angle => 357.5291,
precession => 0.98560028
},
center => [ 1.9148, 0.0200, 0.0003 ],
ecliptic => {
perihelion => 102.9372,
obliquity => 23.45
},
ascension => [ -2.4680, 0.053, -0.0014 ],
declination => [ 22.8008, 0.5999, 0.0493 ],
sidereal => [ 280.1600, 360.9856235 ]
},
mars => {
anomaly => {
angle => 19.3730,
precession => 0.52402068
},
center => [ 10.6912, 0.6228, 0.0503, 0.0046, 0.0005 ],
ecliptic => {
perihelion => 70.9812,
obliquity => 25.19
},
ascension => [ -2.8605, 0.0712, -0.0022 ],
declination => [ 24.387, 0.7331, 0.0706 ],
sidereal => [ 313.4803, 350.89198226 ]
},
jupiter => {
anomaly => {
angle => 20.0202,
precession => 0.08308529
},
center => [ 5.5549, 0.1683, 0.0071, 0.0003 ],
ecliptic => {
perihelion => 237.2074,
obliquity => 3.12
},
ascension => [ -0.0424 ],
declination => [ 3.1151, 0.0015 ],
sidereal => [ 146.0727, 870.5366420 ]
},
saturn => {
anomaly => {
angle => 317.0207,
precession => 0.03344414
},
center => [ 6.3585, 0.2204, 0.0106, 0.0006 ],
ecliptic => {
perihelion => 99.4571,
obliquity => 26.74
},
ascension => [ -3.2364, 0.0911, -0.0031 ],
declination => [ 25.779, 0.8649, 0.0951 ],
sidereal => [ 174.3479, 810.7939024 ]
},
uranus => {
anomaly => {
angle => 141.0498,
precession => 0.01172834
},
center => [ 5.3042, 0.1534, 0.0062, 0.0003 ],
ecliptic => {
perihelion => 5.4639,
obliquity => 82.22
},
# This planet lays nearly on its side.
# The following numbers are nearly useless.
ascension => [ -42.5725, 12.8039, -2.6057 ],
declination => [ 56.9067, -0.8355, 26.1482 ],
sidereal => [ 17.9705, -501.1600928 ]
},
neptune => {
anomaly => {
angle => 256.2250,
precession => 0.00598103
},
center => [ 1.0302, 0.0058 ],
ecliptic => {
perihelion => 182.1957,
obliquity => 27.84
},
ascension => [ -3.5195, 0.1077, -0.0039 ],
declination => [ 26.7577, 0.9662, 0.1164 ],
sidereal => [ 52.3996, 536.3128492 ]
},
pluto => {
anomaly => {
angle => 14.882,
precession => 0.00396
},
center => [ 28.3150, 4.3408, 0.9214, 0.2235, 0.0627, 0.0174 ],
ecliptic => {
perihelion => 4.5433,
obliquity => 57.46
},
# This planet lays nearly on its side.
# The following numbers are nearly useless.
ascension => [ -17.1633, 2.4178, -0.3035 ],
declination => [ 48.3114, 4.788, 4.3582 ],
sidereal => [ 56.3183, -56.3623195 ]
}
);
my %data = %{$planet{$planet}};
# mean anomaly (mu) -- assumes constant orbital speed
my $mu = fmod(
$data{anomaly}{angle} +
$data{anomaly}{precession} * $day,
360
);
#print "mean anomaly (mu): $mu degrees\n";
# true anomaly (nu) -- approximate
my $nu = $mu;
for (my $i = 0; @{$data{center}} > $i; $i++) {
$nu += $data{center}[$i] * sin(($i+1)*d2r($mu));
}
#print "true anomaly (nu): $nu degrees\n";
# ecliptic longitude of sol (lambda)
my $lambda = fmod($nu + $data{ecliptic}{perihelion} + 180, 360);
#print "ecliptic longitude (lambda): $lambda degrees\n";
# right ascension of sol (alpha)
my $alpha = $lambda;
for (my $i = 0; @{$data{ascension}} > $i; $i++) {
$alpha += $data{ascension}[$i] * sin(($i + 1) * 2 * d2r($lambda));
}
#print "right ascension (alpha): $alpha degrees\n";
# declination of sol (delta)
my $delta = r2d(asin(sin(d2r($lambda)) * sin(d2r($data{ecliptic}{obliquity}))));
#print "declination (delta): $delta degrees\n";
# sidereal angle (theta) -- needed for hour angle (Eta)
my $theta = fmod(
$data{sidereal}[0] +
$data{sidereal}[1] * $day - $longitude_west,
360
);
#print "sidereal angle (theta): $theta degrees\n";
# local sidereal time -- needed because it's vexed me for years
#my $lsrtime = $theta * 240;
#print "current local sidereal time: ";
#printf "%02d:%02d:%02d\n", int($lsrtime/3600), ($lsrtime/60%60), ($lsrtime%60);
# hour angle (Eta) -- needed for elevation
my $Eta = fmod($theta - $alpha + 360, 360);
#print "hour angle (Eta): $Eta degrees\n";
# azimuth (Alpha) -- not actually needed
#my $Alpha = fmod(r2d(atan2(r2d(sin(d2r($Eta))), r2d(cos(d2r($Eta)) * sin(d2r($phi)) - tan(d2r($delta)) * cos(d2r($phi))))) + 180, 360);
#print "azimuth (A): $Alpha degrees\n";
# elevation (eta) -- this is the whole point of this file
my $eta = r2d(asin(sin(d2r($phi)) * sin(d2r($delta)) + cos(d2r($phi)) * cos(d2r($delta)) * cos(d2r($Eta))));
#print "true elevation (eta): $eta degrees\n";
#printf "true elevation (eta): %0.4f degrees\n", $eta;
# this bit previously assumed that the sun's angular radius (from earth) was 0.833 degrees.
# wolfram alpha and wikipedia both suggest a value closer to 0.266 degrees (ish).
# perhaps one day i'll bother calculating it based on the actual size of Sol and the
# distance from the planet in question at the current orbital eccentricity.
print "Sol is currently ".(
abs($eta) > 0.266 ? ($eta < 0 ? "not " : '')."visible" :
((localtime)[2] > 12 ? "setting" : "rising")." ".
($eta > 0 ? "above" : "below")." the horizon"
);
printf " (True Elevation %0.${precision}f degrees)\n", $eta;