#!/usr/bin/env -S awk -f
# s2k - display position of the sun using Meeus Astronomical Algorithms.
# s2k is Copyright Daniel K. Allen, 2000-2024.
# All rights reserved.
#
# 6 Dec 2000 - Created by Dan Allen.
# 24 Dec 2000 - Updated to 2nd edition of Meeus, p. 163.
# 3 Jan 2001 - Fixed bug in azimuth calculation; equation of time added.
# 3 Jul 2001 - Time zone set for west coast by a month heuristic.
# 5 Jul 2001 - Ported to Mac OS X, added tz and current date/time.
# 2 Sep 2001 - More accurate Equation of Time.
# 21 Mar 2002 - eqt & sunLong shown mod 360.
# 14 Oct 2002 - Uses apparent longitude & obliquity of the ecliptic for greater accuracy.
# 19 Apr 2004 - North Bend location updated.
# 5 May 2004 - Prints sunLong, not appLong.
# 12 Dec 2004 - Default now Spring Lake; tz estimate improved.
# 9 Jun 2015 - Default now Kirkland.
# 21 Jun 2016 - Handles US time zones approximately; default now West Yarmouth.
# 13 Oct 2016 - Uses env.
# 18 Oct 2016 - No longer uses env due to -S being non-portable.
# 18 Oct 2017 - Default now 400 N Flathead Ln Idaho Falls ID.
# 23 Apr 2020 - Uses env -S. Cygwin loses.
# 31 Oct 2020 - Added Shift(); handles Munich, Yulara, Tristan.
# 2 Aug 2021 - Detabbed. (DRA died 10 years ago today!)
# 25 Aug 2021 - Perl is better for Xcode.
# 14 Jul 2024 - Added Tofino, Anacortes; use strftime rather than calc.
#
# Xcode: Perl for color syntax, Monaco 8 for PDF.
#
BEGIN { # all arguments and results are in decimal degrees
CONVFMT = OFMT = "%10.6f"
PI = 4*atan2(1,1)
D2R = PI/180
R2D = 180/PI
if (lat == "") lat = 43.49924
if (lon == "") lon = 111.94135
if (ARGC > 1 && match(ARGV[1],/^\-/)) { # s2k - gives help
print "# s2k [opts...] [mm.ddyyyy | mm/dd/yyyy] [hh[:mm[:ss]]] [tz] - Sun 2000"
print "# [-v lat=y.y] [-v lon=x.x] - set decimal lat/lon, N/W positive"
print "# [Bog|Ein|Mu] [Ay|Yu] [TdC|Tris] - defaults to IF or pick 1 of these"
exit 1
}
if (match(ARGV[1],/Tof/)) { # Viewing Deck at Middle Beach Lodge, Tofino, BC
lat = 49.1336; lon = 125.90965; Shift()
}
if (match(ARGV[1],/Ana/)) { # Anacortes Marina, slip F-81, Anacortes, WA
lat = 48.5025; lon = 122.602816666667; Shift()
}
if (match(ARGV[1],/Bog|Ein|Mu/)) { # Einsteinstrasse 130, Bogenhausen, Munich, DE
lat = 48.13655; lon = -11.61268; Shift()
}
if (match(ARGV[1],/Ay|Yu/)) { # Yulara Camel Farm at Ayers Rock, Northern Territory, AU
lat = -25.25617; lon = -130.98636; Shift()
}
if (match(ARGV[1],/TdC|Tris/)) { # Albatross Bar, Tristan da Cunha, South Atlantic
lat = -37.06754; lon = 12.31071; Shift()
}
if (ARGC == 1) {
split(strftime("%m/%d/%Y %H:%M:%S"),date)
t = J2000(date[1],date[2])
}
else if (ARGC == 2 && (index(ARGV[1],"/") || index(ARGV[1],".")))
t = J2000(ARGV[1],12)
else if (ARGC == 3 && (index(ARGV[1],"/") || index(ARGV[1],".")))
t = J2000(ARGV[1],ARGV[2])
else if (ARGC == 2 && !index(ARGV[1],"-")) {
split(strftime("%m/%d/%Y %H:%M:%S"),date)
t = J2000(date[1],ARGV[1])
}
else if (ARGV == 4) {
tz = ARGV[3]
t = J2000(ARGV[1],ARGV[2])
}
l = Mod(280.46646 + 36000.76983*t + 0.0003032*t*t,360) # mean longitude
m = Mod(357.52911 + 35999.05029*t + -0.0001537*t*t,360) # mean anomaly
e = (0.016708634 + -0.000042037*t + -0.0000001267*t*t) # eccentricity
c = (1.914602 + -0.004817*t + -0.000014*t*t) * Sin(m) # equation of center
c += (0.019993 + -0.000101*t) * Sin(m*2)
c += (0.000289) * Sin(m*3)
sunLong = Mod(l + c,360)
trueAnom = m + c
radius = (1.000001018 * (1 - e*e)) / (1 + e * Cos(trueAnom))
ecliptic = (84381.448 + -46.815*t - 0.00059*t*t + 0.001813*t*t*t) / 3600
omega = 125.04 - 1934.136*t
apparentLong = sunLong - 0.00569 - 0.00478*Sin(omega)
apparentEcliptic = ecliptic + 0.00256 * Cos(omega)
ra = Mod(ATan2(Cos(apparentEcliptic)*Sin(apparentLong),Cos(apparentLong)),360)
dec = ASin(Sin(apparentEcliptic)*Sin(apparentLong))
gst = (280.46061837 + 0.000387933*t*t + -2.58331180573E-8*t*t*t)
gst += t*36525*360.985647366
gst = Mod(gst,360)
lha = gst - lon - ra
eqt = (sunLong - ra) - (sunLong - l) # positive values occur before UT
if (eqt > 300) eqt = eqt - 360
alt = ASin(Sin(lat) * Sin(dec) + Cos(lat) * Cos(dec) * Cos(lha))
az = 180+ATan2(Sin(lha),(Cos(lha) * Sin(lat) - Cos(lat)*Tan(dec)))
print " Lon:",lon," Lat:",lat
print "SnLo:",sunLong, " R:",radius
print " GHA:",Mod(lha+lon,360)," EQT:",eqt
print " RA:",ra, " Dec:",dec
print " AZ:",az, " Alt:",alt
}
function Floor(x) { return x < 0 ? int(x) - 1 : int(x) }
function Mod(x,y) { return x - y * Floor(x/y) }
function Sin(x) { return sin(x*D2R) }
function Cos(x) { return cos(x*D2R) }
function Tan(x) { return Sin(x)/Cos(x) }
function ASin(x) { return atan2(x,sqrt(1 - x * x))*R2D }
function ATan2(y,x){ return atan2(y,x)*R2D }
function J2000(date,time, m,d,y,a) {
if (index(date,"/") > 0) { # mm/dd/yy or mm/dd/yyyy - Y2K compatible
split(date,a,"/")
m = a[1]
d = a[2]
y = a[3] < 50 ? 2000 + a[3] : a[3] < 100 ? 1900 + a[3] : a[3]
delete a
}
else if (index(date,".") > 0) { # mm.ddyyyy - HP calculator style format
m = int(date)
d = int((date - m)*100)
y = int(1000000*(date - int(date) - d/100)+0.5)
}
if (tz == "") { # ~ US DST tz code
tz = m > 3 && m < 11 ? 4 : 5 # EDT or EST
if (lon < -7.5) tz -= 6 # CEST or CET
if (lon < -130 && lon > -145) tz = -9.5 # ACST (no DST)
if (lon > 87) tz++ # CDT or CST
if (lon > 102) tz++ # MDT or MST
if (lon > 114.5) tz++ # PDT or PST
}
split(time,a,":")
print "Time:",y,m,d,a[1],a[2],a[3]," TZ:",tz
return (Julian(y,m,d,a[1]+tz,a[2],a[3]) - Julian(2000,1,1,12,0,0))/36525
}
function Julian(year,month,day,hr,min,sec, m,y,a,b,jd)
{
if (month <= 2) { y = year - 1; m = month + 12 }
else { y = year; m = month }
a = int(y/100)
if (y < 1582 || y == 1582 && (m < 10 || m == 10 && day <=4)) b = 0
else b = 2 - a + int(a/4)
jd = int(365.25*(y+4716)) + int(30.6001*(m+1)) + day + b - 1524.5
jd += (hr-12)/24.0 + min/(24.0*60) + sec/(24.0*3600)
return jd
}
function Shift()
{
ARGV[1] = ARGV[2]; ARGV[2] = ARGV[3]; ARGV[3] = ARGV[4]
ARGV[4] = ARGV[5]; ARGC--
}
Back to Dan Allen's home page.
Created: 27 Feb 2001
Modified: 14 Jul 2024