phase.java

java高级使用教程 全书一共分六章

// Phase - phase of the moon calculations
//
// Adapted from "moontool.c" by John Walker, Release 2.0.
//
// Copyright (C)1996,1998 by Jef Poskanzer <[e-mail]>.  All rights reserved.
//
// 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.
//
// THIS SOFTWARE IS PROVIDED BY THE AUTHOR 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 AUTHOR 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.
//
// Visit the ACME Labs Java page for up-to-date versions of this and other
// fine Java utilities: http://www.acme.com/java/

package Acme;

import java.util.*;

/// Phase - phase of the moon calculations
// <P>
// Adapted from "moontool.c" by John Walker, Release 2.0.
// <P>
// <A HREF="/resources/classes/Acme/Phase.java">Fetch the software.</A><BR>
// <A HREF="/resources/classes/Acme.tar.gz">Fetch the entire Acme package.</A>

public class Phase
    {

    // Astronomical constants.

    // 1980 January 0.0
    private static final double epoch = 2444238.5D;

    // Constants defining the Sun's apparent orbit.

    // ecliptic longitude of the Sun at epoch 1980.0
    private static final double elonge = 278.833540D;

    // ecliptic longitude of the Sun at perigee
    private static final double elongp = 282.596403D;

    // eccentricity of Earth's orbit
    private static final double eccent = 0.016718D;

    // semi-major axis of Earth's orbit, km
    private static final double sunsmax = 1.495985e8D;

    // sun's angular size, degrees, at semi-major axis distance
    private static final double sunangsiz = 0.533128D;

    // Elements of the Moon's orbit, epoch 1980.0.

    // moon's mean lonigitude at the epoch
    private static final double mmlong = 64.975464D;

    // mean longitude of the perigee at the epoch
    private static final double mmlongp = 349.383063D;

    // mean longitude of the node at the epoch
    private static final double mlnode = 151.950429D;

    // inclination of the Moon's orbit
    private static final double minc = 5.145396D;

    // eccentricity of the Moon's orbit
    private static final double mecc = 0.054900D;

    // moon's angular size at distance a from Earth
    private static final double mangsiz = 0.5181D;

    // semi-major axis of Moon's orbit in km
    private static final double msmax = 384401.0D;

    // parallax at distance a from Earth
    private static final double mparallax = 0.9507D;

    // synodic month (new Moon to new Moon)
    private static final double synmonth = 29.53058868D;

    // base date for E. W. Brown's numbered series of lunations (1923 Jan 16)
    private static final double lunatbase = 2423436.0D;

    // Properties of the Earth.

    // radius of Earth in kilometres
    private static final double earthrad = 6378.16D;

    // Mathematical constants.
    private static final double EPSILON = 1E-6D;

    // Handy mathematical functions.

    // Fix angle.
    private static double fixangle( double a )
	{
	double b = a - 360.0 * Math.floor( a / 360.0D );
	// Can't use Math.IEEEremainder here because remainder differs
	// from modulus for negative numbers.
	return b;
	}

    // Degrees to radians.
    private static double torad( double d )
	{
	return d * Math.PI / 180.0D;
	}

    // Radians to degrees.
    private static double todeg( double r )
	{
	return r * 180.0D / Math.PI;
	}

    // Sin from degrees.
    private static double dsin( double d )
	{
	return Math.sin( torad( d ) );
	}

    // Cos from degrees.
    private static double dcos( double d )
	{
	return Math.cos( torad( d ) );
	}

    // jdate - convert internal GMT date to Julian day.
    private static long jdate( Date t )
	{
	long c, m, y;

	Calendar cal = new GregorianCalendar();
	cal.setTime( t );
	y = cal.get( Calendar.YEAR ) + 1900;
	m = cal.get( Calendar.MONTH ) + 1;
	if ( m > 2 )
	    m = m - 3;
	else
	    {
	    m = m + 9;
	    --y;
	    }
	c = y / 100L;		// compute century
	y -= 100L * c;
	return cal.get( Calendar.DATE ) +
	  ( c * 146097L ) / 4 + ( y * 1461L ) / 4 +
	  ( m * 153L + 2 ) / 5 + 1721119L;
	}

    /// Convert internal date and time to astronomical Julian
    // time (i.e. Julian date plus day fraction, expressed as
    // a double).
    public static double jtime( Date t )
	{
	int c;

	Calendar cal = new GregorianCalendar();
	cal.setTime( t );
	c = - cal.get( Calendar.ZONE_OFFSET );	// !!! should this be negative?
	return ( jdate( t ) - 0.5 ) + 
	    ( cal.get( Calendar.SECOND ) + 60 * ( cal.get( Calendar.MINUTE ) + c + 60 * cal.get( Calendar.HOUR_OF_DAY ) ) ) / 86400.0;
	}

    // jyear - convert Julian date to year, month, day, which are
    //     returned via integer pointers to integers
    private static void jyear( double td, RefInt yy, RefInt mm, RefInt dd )
	{
	double j, d, y, m;

	td += 0.5;	// astronomical to civil
	j = Math.floor(td);
	j = j - 1721119.0;
	y = Math.floor(((4 * j) - 1) / 146097.0);
	j = (j * 4.0) - (1.0 + (146097.0 * y));
	d = Math.floor(j / 4.0);
	j = Math.floor(((4.0 * d) + 3.0) / 1461.0);
	d = ((4.0 * d) + 3.0) - (1461.0 * j);
	d = Math.floor((d + 4.0) / 4.0);
	m = Math.floor(((5.0 * d) - 3) / 153.0);
	d = (5.0 * d) - (3.0 + (153.0 * m));
	d = Math.floor((d + 5.0) / 5.0);
	y = (100.0 * y) + j;
	if (m < 10.0)
	    m = m + 3;
	else
	    {
	    m = m - 9;
	    y = y + 1;
	    }
	yy.val = (int) y;
	mm.val = (int) m;
	dd.val = (int) d;
	}

    // meanphase - calculates mean phase of the Moon for a given base date
    //     and desired phase:
    //         0.0   New Moon
    //         0.25  First quarter
    //         0.5   Full moon
    //         0.75  Last quarter
    // Beware!!!  This routine returns meaningless results for any other
    // phase arguments.  Don't attempt to generalise it without understanding
    // that the motion of the moon is far more complicated that this
    // calculation reveals.
    private static double meanphase( double sdate, double phase, RefDouble usek )
	{
	RefInt yy = new RefInt();
	RefInt mm = new RefInt();
	RefInt dd = new RefInt();
	double k, t, t2, t3, nt1;

	jyear( sdate, yy, mm, dd );

	k = (yy.val + ((mm.val - 1) * (1.0 / 12.0)) - 1900) * 12.3685;

	// Time in Julian centuries from 1900 January 0.5.
	t = (sdate - 2415020.0) / 36525;
	t2 = t * t;		   // square for frequent use
	t3 = t2 * t;		   // cube for frequent use

	usek.val = k = Math.floor(k) + phase;
	nt1 = 2415020.75933 + synmonth * k
	      + 0.0001178 * t2
	      - 0.000000155 * t3
	      + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2);

	return nt1;
	}

    // truephase - given a K value used to determine the mean phase of the
    //     new moon, and a phase selector (0.0, 0.25, 0.5, 0.75),
    //     obtain the true, corrected phase time
    private static double truephase( double k, double phase )
	{
	double t, t2, t3, pt, m, mprime, f;
	boolean apcor = false;

	k += phase;		   /* add phase to new moon time */
	t = k / 1236.85;	   /* time in Julian centuries from
				        1900 January 0.5 */