/*
        The following routines come from Moontool for Windows.  Below is
        a statement from the author:
 
                          Moontool for Windows
                               Release 2.0
 
    Designed and implemented by John Walker in December 1987.
    Revised and updated in July of 1989 by Ron Hitchens.
    Converted to Microsoft Windows in January of 1992 by John Walker.
        Convert to Win32 in March of 1999 by John Walker.
 
    The  algorithms used in this program to calculate the positions of
    the Sun and Moon as seen from the Earth  are  given  in  the  book
    "Practical Astronomy With Your Calculator" by Peter Duffett-Smith,
    Second Edition, Cambridge University Press, 1981.  Ignore the word
    "Calculator"  in  the  title;  this  is  an essential reference if
    you're  interested  in  developing   software   which   calculates
    planetary  positions,  orbits,  eclipses, and the like.  If you're
    interested in pursuing such programming, you should  also  obtain:
 
    "Astronomical Formulae  for  Calculators"  by  Jean  Meeus,  Third
    Edition, Willmann-Bell, 1985.  A must-have.
 
    "Planetary  Programs  and  Tables  from  -4000 to +2800" by Pierre
    Bretagnon and Jean-Louis Simon, Willmann-Bell, 1986.  If you  want
    the  utmost  (outside of JPL) accuracy for the planets, it's here.
 
    "Celestial BASIC" by Eric Burgess, Revised Edition,  Sybex,  1985.
    Very cookbook oriented, and many of the algorithms are hard to dig
    out of the turgid BASIC code, but you'll probably want it anyway.
 
    Many  of these references can be obtained from Willmann-Bell, P.O.
    Box 35025, Richmond, VA 23235, USA.  Phone:  (804)  320-7016.   In
    addition  to  their  own  publications,  they  stock  most  of the
    standard references for mathematical and positional astronomy.
 
    This program was written by:
 
        John Walker
                http://www.fourmilab.ch/
 
    This  program is in the public domain: "Do what thou wilt shall be
    the whole of the law".  I'd appreciate  receiving  any  bug  fixes
    and/or  enhancements, which I'll incorporate in future versions of
    the program.  Please leave the  original  attribution  information
    intact so that credit and blame may be properly apportioned.
 
        History:
        --------
        June 1988       Version 2.0 for the Sun workstation posted
                        to usenet by John Walker
 
        June 1988       Modified by Ron Hitchens
 
        July 1989       Modified a little bit more to use an accurate
                        grey-scale moon face created by Joe Hitchens
                        on an Amiga.
                        Added The Apollo 11 Commemorative Red Dot, to show
                        where Neil and Buzz went on vacation 20 years ago.
 
        March 1992      Moontool for Windows 1.0 implemented by John Walker.
 
        April 1992      Bug fix update 1.01 correcting problems reported by
                        Michael Geary.
 
                March 1999		Win32 version, reverting to standard Julian Day
                                                definition after ill-conceived flirtation with "civil
                                                Julian Day" beginning at midnight.  Release 2.0.
 
 */
 
	/*  Astronomical constants  */
    
    var epoch  =     2444238.5;      /* 1980 January 0.0 */
    
    /*  Constants defining the Sun's apparent orbit  */
    
    var  elonge  =    278.833540;     /* Ecliptic longitude of the Sun
                                      at epoch 1980.0 */
    var  elongp  =    282.596403;     /* Ecliptic longitude of the Sun at
                                      perigee */
    var  eccent  =    0.016718;       /* Eccentricity of Earth's orbit */
    var  sunsmax  =   1.495985e8;     /* Semi-major axis of Earth's orbit, km */
    var  sunangsiz =  0.533128;       /* Sun's angular size, degrees, at
                                      semi-major axis distance */
    
    /*  Elements of the Moon's orbit, epoch 1980.0  */
    
    var  mmlong  =    64.975464;      /* Moon's mean longitude at the epoch */
    var  mmlongp  =   349.383063;     /* Mean longitude of the perigee at the
                                      epoch */
    var  mlnode =     151.950429;     /* Mean longitude of the node at the
                                      epoch */
    var  minc =       5.145396;       /* Inclination of the Moon's orbit */
    var  mecc =       0.054900;       /* Eccentricity of the Moon's orbit */
    var mangsiz =    0.5181;         /* Moon's angular size at distance a
                                      from Earth */
    var msmax  =     384401.0;       /* Semi-major axis of Moon's orbit in km */
    var mparallax =   0.9507;         /* Parallax at distance a from Earth */
    var synmonth =   29.53058868;    /* Synodic month (new Moon to new Moon) */
    var lunatbase  = 2423436.0;      /* Base date for E. W. Brown's numbered
                                      series of lunations (1923 January 16) */
    
    /*  Properties of the Earth  */
    
    var PI =3.14159265358979323846;  /* Assume not near black hole nor in
                                      Tennessee */
    
    var EPSILON = 1E-6;
    
	function torad(d) {
		return ((d) * (PI / 180.0))
	}

	function todeg(d) {
		return ((d) * (180.0 / PI)) 
	}
    
	/*  UCTTOJ  --	Convert GMT date and time to astronomical
					Julian time (i.e. Julian date plus day fraction,
					expressed as a double).  */

	function ucttoj(year, mon, mday, hour, min, sec)
	{

		// Algorithm as given in Meeus, Astronomical Algorithms, Chapter 7, page 61

		var a;
		var b;
		var m;
		var y;


		m = mon + 1;
		y = year;

		if (m <= 2) {
			y--;
			m += 12;
		}

		/* Determine whether date is in Julian or Gregorian calendar based on
		   canonical date of calendar reform. */

		if ((year < 1582) || ((year == 1582) && ((mon < 9) || (mon == 9 && mday < 5)))) {
			b = 0;
		} else {
			a = Math.floor(y / 100);
			b = 2 - a + (a / 4);
		}

		return ((Math.floor(365.25 * (y + 4716))) + (Math.floor(30.6001 * (m + 1))) +
					mday + b - 1524.5) +
				((sec + 60 * (min + 60 * hour)) / 86400.0);
	}


    /*  KEPLER  --   Solve the equation of Kepler.  */
    
    function kepler( m, ecc) {
        var e;
		var delta;
        
        m = torad(m);
		e = m
        do {
            delta = e - ecc * Math.sin(e) - m;
            e -= delta / (1 - ecc * Math.cos(e));
        } while (Math.abs(delta) > EPSILON);
        return e;
    }
    
/*  PHASE  --  Calculate phase of moon as a fraction:
 
    The  argument  is  the  time  for  which  the  phase is requested,
    expressed as a Julian date and fraction.  Returns  the  terminator
    phase  angle  as a percentage of a full circle (i.e., 0 to 1), and
    stores into pointer arguments  the  illuminated  fraction  of  the
    Moon's  disc, the Moon's age in days and fraction, the distance of
    the Moon from the centre of the Earth, and  the  angular  diameter
    subtended  by the Moon as seen by an observer at the centre of the
    Earth.
 */
    function fixangle(a) {
        return ((a) - 360.0 * (Math.floor((a) / 360.0)))  ;
    }
    

	
	function CalcPhase( pdate )   {                   /* Date for which to calculate phase */
//        PhaseData data = new PhaseData();
        var Day;
		var N;
		var M;
		var Ec;
		var Lambdasun;
		var ml;
		var MM;
		var MN;
		var Ev;
		var Ae;
		var A3;
		var MmP;
		var mEc;
		var A4;
		var lP;
		var V;
		var lPP;
		var NP;
		var y;
		var x;
		var Lambdamoon;
		var BetaM;
		var MoonAge;
		var MoonPhase;
		var MoonDist;
		var MoonDFrac;
		var MoonAng;
		var MoonPar;
		var F;
		var SunDist;
		var SunAng;
        
        /* Calculation of the Sun's position */
        
        Day = pdate - epoch;                    /* Date within epoch */
        N = fixangle((360 / 365.2422) * Day);   /* Mean anomaly of the Sun */
        M = fixangle(N + elonge - elongp);      /* Convert from perigee
                                               co-ordinates to epoch 1980.0 */
        Ec = kepler(M, eccent);                 /* Solve equation of Kepler */
        Ec = Math.sqrt((1 + eccent) / (1 - eccent)) * Math.tan(Ec / 2);
        Ec = 2 * todeg(Math.atan(Ec));               /* True anomaly */
        Lambdasun = fixangle(Ec + elongp);      /* Sun's geocentric ecliptic
                                               longitude */
        /* Orbital distance factor */
        F = ((1 + eccent * Math.cos(torad(Ec))) / (1 - eccent * eccent));
        SunDist = sunsmax / F;                  /* Distance to Sun in km */
        SunAng = F * sunangsiz;                 /* Sun's angular size in degrees */
        
        /* Calculation of the Moon's position */
        
        /* Moon's mean longitude */
        ml = fixangle(13.1763966 * Day + mmlong);
        
        /* Moon's mean anomaly */
        MM = fixangle(ml - 0.1114041 * Day - mmlongp);
        
        /* Moon's ascending node mean longitude */
        MN = fixangle(mlnode - 0.0529539 * Day);
        
        /* Evection */
        Ev = 1.2739 * Math.sin(torad(2 * (ml - Lambdasun) - MM));
        
        /* Annual equation */
        Ae = 0.1858 * Math.sin(torad(M));
        
        /* Correction term */
        A3 = 0.37 * Math.sin(torad(M));
        
        /* Corrected anomaly */
        MmP = MM + Ev - Ae - A3;
        
        /* Correction for the equation of the centre */
        mEc = 6.2886 * Math.sin(torad(MmP));
        
        /* Another correction term */
        A4 = 0.214 * Math.sin(torad(2 * MmP));
        
        /* Corrected longitude */
        lP = ml + Ev + mEc - Ae + A4;
        
        /* Variation */
        V = 0.6583 * Math.sin(torad(2 * (lP - Lambdasun)));
        
        /* True longitude */
        lPP = lP + V;
        
        /* Corrected longitude of the node */
        NP = MN - 0.16 * Math.sin(torad(M));
        
        /* Y inclination coordinate */
        y = Math.sin(torad(lPP - NP)) * Math.cos(torad(minc));
        
        /* X inclination coordinate */
        x = Math.cos(torad(lPP - NP));
        
        /* Ecliptic longitude */
        Lambdamoon = todeg(Math.atan2(y, x));
        Lambdamoon += NP;
        
        /* Ecliptic latitude */
        BetaM = todeg(Math.asin(Math.sin(torad(lPP - NP)) * Math.sin(torad(minc))));
        
        /* Calculation of the phase of the Moon */
        
        /* Age of the Moon in degrees */
        MoonAge = lPP - Lambdasun;
        
        /* Phase of the Moon */
        MoonPhase = (1 - Math.cos(torad(MoonAge))) / 2;
        
        /* Calculate distance of moon from the centre of the Earth */
        
        MoonDist = (msmax * (1 - mecc * mecc)) /
        (1 + mecc * Math.cos(torad(MmP + mEc)));
        
        /* Calculate Moon's angular diameter */
        
        MoonDFrac = MoonDist / msmax;
        MoonAng = mangsiz / MoonDFrac;
        
        /* Calculate Moon's parallax */
        
        MoonPar = mparallax / MoonDFrac;
        
//        data.pphase = MoonPhase;
//        data.mage = synmonth * (fixangle(MoonAge) / 360.0);
//		alert(synmonth * (fixangle(MoonAge) / 360.0));
//        data.dist = MoonDist;
//        data.angdia = MoonAng;
//        data.sudist = SunDist;
//        data.suangdia = SunAng;
//        data.perc = fixangle(MoonAge) / 360.0;
        return synmonth * (fixangle(MoonAge) / 360.0);
    }