## -: pom.c; $NetBSD; pom.c,v 1.14 2004/01/27 20:30:30 jsm Exp $

# Copyright (c) 1989, 1993
#    The Regents of the University of California.  All rights reserved.
#
# This code is derived from software posted to USENET.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
# 1. Redistributions of source code must retain the above copyright
#    notice, this list of conditions and the following disclaimer.
# 2. Redistributions in binary form must reproduce the above copyright
#    notice, this list of conditions and the following disclaimer in the
#    documentation and/or other materials provided with the distribution.
# 3. Neither the name of the University nor the names of its contributors
#    may be used to endorse or promote products derived from this software
#    without specific prior written permission.

# THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
# OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
# SUCH DAMAGE.

# usage: "$daytime" http://antwrp.gsfc.nasa.gov/apod/
BEGIN {
    PI = 3.14159265358979323846
    EPOCH_MINUS_1970 = (20 * 365 + 5 - 1) # 20 years, 5 leaps, back 1 day to Jan 0
    EPSILONg= 279.403303    # solar ecliptic long at EPOCH
    RHOg    = 282.768422    # solar ecliptic long of perigee at EPOCH
    ECCEN   = 0.016713      # solar orbit eccentricity
    lzero   = 318.351648    # lunar mean long at EPOCH
    Pzero   = 36.340410     # lunar mean long of perigee at EPOCH
### Nzero   = 318.510107    # lunar mean long of node at EPOCH
    baseurl = "https://apod.nasa.gov/apod/"
    imageurl = "image url"
    now = tmpt = systime()
    if (1 < ARGC && ARGV[1] !~ /:\/\//) {
        for (i = 1; i < ARGC && ARGV[i] !~ /:\/\//; i++) {
            aday = aday ARGV[i] " "; ARGV[i] = ""	# (引用符で囲まれていない)引数を連結する
        }
        tmpt = parsetime(aday)
    }
    days = (tmpt - EPOCH_MINUS_1970 * 86400) / 86400.0
    today = potm(days) + .5
    if (tmpt < now)
        result = "The Moon was "
    else if (tmpt == now)
        result = "The Moon is "
    else
        result = "The Moon will be "
    if (int(today) == 100)
        result = result "Full"
    else if (!int(today))
        result = result "New"
    else {
        tomorrow = potm(days + 1)
        if (int(today) == 50)
            result = result ((tomorrow > today) ? "at the First Quarter" : "at the Last Quarter")
            # today is 0.5 too big, but it doesn't matter here
            # since the phase is changing fast enough
        else {
            today -= 0.5    # Now it might matter
            result = result ((tomorrow > today) ? "Waxing " : "Waning ")
            if (today > 50)
                result = result sprintf("Gibbous (%1.0f%% of Full)", today)
            else if (today < 50)
                result = result sprintf("Crescent (%1.0f%% of Full)", today)
        }
    }
### print result
}
# return phase of the moon
function potm(days,  N, Msol, Ec, LambdaSol, l, Mm, Ev,
        Ac, A3, Mmprime, A4, lprime, V, ldprime, D) {
    N = (360 * days / 365.242191) % 360         # sec 46 #3
    Msol = (N + EPSILONg - RHOg) % 360          # sec 46 #4
    Ec = 360 / PI * ECCEN * sin(dtor(Msol))     # sec 46 #5
    LambdaSol = (N + Ec + EPSILONg) % 360       # sec 46 #6
    l = (13.1763966 * days + lzero) % 360       # sec 65 #4
    Mm = (l - (0.1114041 * days) - Pzero) % 360 # sec 65 #5
### Nm = (Nzero - (0.0529539 * days)) % 360     # sec 65 #6
    Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)) # sec 65 #7
    Ac = 0.1858 * sin(dtor(Msol))               # sec 65 #8
    A3 = 0.37 * sin(dtor(Msol))
    Mmprime = Mm + Ev - Ac - A3                 # sec 65 #9
    Ec = 6.2886 * sin(dtor(Mmprime))            # sec 65 #10
    A4 = 0.214 * sin(dtor(2 * Mmprime))         # sec 65 #11
    lprime = l + Ev + Ec - Ac + A4              # sec 65 #12
    V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))) # sec 65 #13
    ldprime = lprime + V                        # sec 65 #14
    D = ldprime - LambdaSol                     # sec 67 #2
    return 50.0 * (1 - cos(dtor(D)))            # sec 67 #3
}
# convert degrees to radians
function dtor(deg) {
    return deg * PI / 180
}
#
function parsetime(aday,  atime) {
    gsub(/["'_]/, " ", aday)
    atime = mktime(aday)
    if (-1 == atime) atime = systime()		# parse error
### printf("%s:  ", strftime("yyyy/MM/dd HH:mm:ss", atime))
    return atime
}
# NASA (Astronomy Picture of the Day)
	{
	if (sub(/.*<IMG SRC=\"/, "")) {
		baseurl = FILENAME; sub(/[^\/]+\.html$/, "", baseurl)
		sub(/\".*$/, ""); imageurl = $0
		exit(0)
	}
}
END	{
	print "<html><body link='cyan' vlink='cyan' style='background-color: black; " \
		"width: 1024px; height: 768px; background-image: url(" baseurl imageurl \
		");'><h4><a href='" FILENAME "'>" result "</a></h4></body></html>"
}
