/*
 * $Id: pom.c,v 1.8 1995/01/01 19:31:36 ceder Exp $
 * Copyright (C) 1991, 1993, 1994  Lysator Academic Computer Association.
 *
 * This file is part of the LysKOM server.
 * 
 * LysKOM is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by 
 * the Free Software Foundation; either version 1, or (at your option) 
 * any later version.
 * 
 * LysKOM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with LysKOM; see the file COPYING.  If not, write to
 * Lysator, c/o ISY, Linkoping University, S-581 83 Linkoping, SWEDEN,
 * or the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, 
 * MA 02139, USA.
 *
 * Please mail bug reports to bug-lyskom@lysator.liu.se. 
 */
/*
 *  pom.c  --  Calculate the Phase Of the Moon
 *
 *  Written by G|ran ]hlstr|m 1984-04-05, revision 5 1986-03-03
 *
 *  Modified by Thomas Bellman 1991-01-29
 */

static char *rcsid = "$Id: pom.c,v 1.8 1995/01/01 19:31:36 ceder Exp $";

#include <sys/types.h>
#include <time.h>
#include <math.h>
#if defined(__svr4__) && defined(__sparc__)
#include <values.h>
#endif
#include <stdio.h>

#include "pom.h"
#include "rcs.h"
USE(rcsid);


#define EXPORT

#ifndef M_PI
#define M_PI 3.141592653
#endif
#define RADS	(M_PI/180.0)

#if 0
static	char		* pmt[] = { "NM", "FQ", "FM", "LQ" };
#endif
static	short		  ix;
static	double		  m;		/* Sun's Mean Anomaly */
static	double		  mp;		/* Moon's Mean Anomaly */
static	double		  f;		/* Moon's Argument of Latitude	*/
static	struct tm	  dt;


/* New & Full Moon correction */
static  double
mcorr (double	t)

{
    double	  tmp;

    tmp = (0.1734-3.93E-4*t) * sin(m) + 0.0021 * sin(2.0*m) - 0.4068 * sin(mp);
    tmp += 0.0161 * sin(2.0*mp) - 0.0004 * sin(3.0*mp) + 0.0104 * sin(2.0*f);
    tmp += 0.0004 * sin(2.0*f + m) - 0.0051 * sin(m+mp) - 0.0074 * sin(m-mp);
    tmp -= 0.0004 * sin(2.0*f - m) + 0.0006 * sin(2.0*f + mp);
    return  tmp += 0.0010 * sin(2.0*f - mp) + 0.0005 * sin(m + 2.0*mp);
}



/*  First & Last Quarter correction  */
static  double
qcorr (double	t)

{
    double	  tmp, tmp2;

    tmp = (0.1721 - 0.0004*t) * sin(m) + 0.0021 * sin(2.0*m) - 0.628 * sin(mp);
    tmp +=  0.0089 * sin(2.0*mp) - 0.0004 * sin(3.0*mp) + 0.0079 * sin(2.0*f);
    tmp +=  0.0003 * sin(2.0*f + m) - 0.0119 * sin(m+mp) - 0.0047 * sin(m-mp);
    tmp +=  0.0021 * sin(2.0*f - mp) - 0.0004 * sin(2.0*f - m);
    tmp += -0.0006 * sin(2.0*f + mp) + 0.0003 * sin(m + 2.0*mp);
    tmp +=  0.0004 * sin(m - 2.0*mp) - 0.0003 * sin(2.0*m + mp);
    tmp2 = (0.0028 - 0.0004 * cos (m) + 0.0003 * cos (mp));
    if (ix == 1)
	tmp += tmp2;
    else
	tmp -= tmp2;

    return  tmp;
}



/* When is nearest phase? */
static  double
pcon (void)

{
    double	  frc;
    double	  year;
    double	  k;
    double	  t;
    double	  tmp;
    int		  leap;


    year = (dt.tm_hour + ((dt.tm_min + (dt.tm_sec / 60.0)) / 60.0)) / 24.0;
    year += dt.tm_yday;

    /* Leap year? */
    leap = dt.tm_year + 1900;
    leap = (leap % 4 == 0  &&  (leap % 100 != 0  ||  leap % 400 == 0));

    /* This year with decimals */
    year = dt.tm_year + year / (365.0 + leap);

    k = year * 12.3685;
    k -= (frc = k - floor (k));

    /* Determine nearest phase */
    for ( ix = 0 ;  frc > ix * 0.25 + 0.125 ;  ix++ )
	;
    t = (k += ix * 0.25) / 1236.85;
    tmp = 2415020.75933 + 29.53058868 * k
	+ (1.178E-4 - 1.55E-7 * t) * pow (t, 2.0);
    tmp += 3.3E-4 * sin (166.56*RADS + 132.87*RADS * t
			 - 9.173E-3*RADS * pow(t, 2.0));
    m = RADS * fmod (359.2242 + 29.10535608 * k
		     - (3.47e-6*t + 3.33e-5) * pow (t, 2.0), 360.0);
    mp = RADS * fmod (306.0253 + 385.81691806 * k
		      + (1.236e-5*t + 0.0107306) * pow (t, 2.0), 360.0);
    f = RADS * fmod (21.2964 + 390.67050646 * k
		     - (2.39e-6*t + 1.6528e-3) * pow (t, 2.0), 360.0);
    return  tmp += ((ix % 2 == 0) ? mcorr (t) : qcorr (t));
}



/* JD at this instant */
static  double
dcon (void)

{
    double	  a, b, d, m, y;


    m = dt.tm_mon;
    y = dt.tm_year + 1900.0;
    d = dt.tm_mday + dt.tm_hour / 24.0 + dt.tm_min / 1440.0
	+ dt.tm_sec / 86400.0;
    if (m <= 2.0)
    {
	y -= 1.0;
	m += 12.0;
    }
    a = floor (y/100.0);	/* Gregorian Calendar correction */
    b = 2.0 - a + floor (a / 4.0);
    return  (floor (365.25 * y) + floor (30.6001 * (m+1.0)) + d + 1720994.5 + b);
}



EXPORT  struct pom
phase_of_the_moon (struct tm	* when)

{
    int		  p1, p2, p3, p4;
    char	  sgn;
    double	  tmp;
    struct pom	  result;


    dt = *when;
    dt.tm_mon++;
    dt.tm_yday++;
    tmp = pcon() - dcon();
    sgn = (tmp > 0)  ?  '-'  :  '+';
    tmp = fabs (tmp);
    tmp -= (p1 = floor (tmp));
    tmp = tmp * 24.0 - (p2 = floor (tmp * 24.0));
    tmp = tmp * 60.0 - (p3 = floor (tmp * 60.0));
    p4 = floor (tmp * 60.0);

    result.quarter = ix % 4;
    result.days    = (sgn == '+')  ?  p1  :  -p1;
    result.hours   = (sgn == '+')  ?  p2  :  -p2;
    result.minutes = (sgn == '+')  ?  p3  :  -p3;
    result.seconds = (sgn == '+')  ?  p4  :  -p4;

    return  result;
}