phase.java

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

	t2 = t * t;		   /* square for frequent use */
	t3 = t2 * t;		   /* cube for frequent use */
	pt = 2415020.75933	   /* mean time of phase */
	     + synmonth * k
	     + 0.0001178 * t2
	     - 0.000000155 * t3
	     + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2);

        m = 359.2242               /* Sun's mean anomaly */
	    + 29.10535608 * k
	    - 0.0000333 * t2
	    - 0.00000347 * t3;
        mprime = 306.0253          /* Moon's mean anomaly */
	    + 385.81691806 * k
	    + 0.0107306 * t2
	    + 0.00001236 * t3;
        f = 21.2964                /* Moon's argument of latitude */
	    + 390.67050646 * k
	    - 0.0016528 * t2
	    - 0.00000239 * t3;
	if ((phase < 0.01) || (Math.abs(phase - 0.5) < 0.01))
	    {
	    /* Corrections for New and Full Moon. */
	    pt +=     (0.1734 - 0.000393 * t) * dsin(m)
		     + 0.0021 * dsin(2 * m)
		     - 0.4068 * dsin(mprime)
		     + 0.0161 * dsin(2 * mprime)
		     - 0.0004 * dsin(3 * mprime)
		     + 0.0104 * dsin(2 * f)
		     - 0.0051 * dsin(m + mprime)
		     - 0.0074 * dsin(m - mprime)
		     + 0.0004 * dsin(2 * f + m)
		     - 0.0004 * dsin(2 * f - m)
		     - 0.0006 * dsin(2 * f + mprime)
		     + 0.0010 * dsin(2 * f - mprime)
		     + 0.0005 * dsin(m + 2 * mprime);
	    apcor = true;
	    }
	else if ((Math.abs(phase - 0.25) < 0.01 || (Math.abs(phase - 0.75) < 0.01)))
	    {
	    pt +=     (0.1721 - 0.0004 * t) * dsin(m)
		     + 0.0021 * dsin(2 * m)
		     - 0.6280 * dsin(mprime)
		     + 0.0089 * dsin(2 * mprime)
		     - 0.0004 * dsin(3 * mprime)
		     + 0.0079 * dsin(2 * f)
		     - 0.0119 * dsin(m + mprime)
		     - 0.0047 * dsin(m - mprime)
		     + 0.0003 * dsin(2 * f + m)
		     - 0.0004 * dsin(2 * f - m)
		     - 0.0006 * dsin(2 * f + mprime)
		     + 0.0021 * dsin(2 * f - mprime)
		     + 0.0003 * dsin(m + 2 * mprime)
		     + 0.0004 * dsin(m - 2 * mprime)
		     - 0.0003 * dsin(2 * m + mprime);
	    if (phase < 0.5)
	        /* First quarter correction. */
	        pt += 0.0028 - 0.0004 * dcos(m) + 0.0003 * dcos(mprime);
	    else
	        /* Last quarter correction. */
	        pt += -0.0028 + 0.0004 * dcos(m) - 0.0003 * dcos(mprime);
	    apcor = true;
	    }
	if (!apcor)
            throw new InternalError( "Acme.Phase.truephase() called with invalid phase selector" );
	return pt;
	}

    /// Find time of phases of the moon which surround the current
    // date.  Five phases are found, starting and ending with the
    // new moons which bound the current lunation.
    public static void phasehunt5( double sdate, double[] phases )
	{
	double adate, nt1, nt2;
	RefDouble k1 = new RefDouble();
	RefDouble k2 = new RefDouble();

	adate = sdate - 45;
	nt1 = meanphase(adate, 0.0, k1);
	for (;;)
	    {
	    adate += synmonth;
	    nt2 = meanphase(adate, 0.0, k2);
	    if (nt1 <= sdate && nt2 > sdate)
	        break;
	    nt1 = nt2;
	    k1.val = k2.val;
	    }
	phases[0] = truephase(k1.val, 0.0);
	phases[1] = truephase(k1.val, 0.25);
	phases[2] = truephase(k1.val, 0.5);
	phases[3] = truephase(k1.val, 0.75);
	phases[4] = truephase(k2.val, 0.0);
	}

    // phasehunt2 - find time of phases of the moon which surround the current
    //     date.  Two phases are found.
    void phasehunt2( double sdate, double[] phases, double[] which )
	{
	double adate, nt1, nt2;
	RefDouble k1 = new RefDouble();
	RefDouble k2 = new RefDouble();

	adate = sdate - 45;
	nt1 = meanphase(adate, 0.0, k1);
	for (;;)
	    {
	    adate += synmonth;
	    nt2 = meanphase(adate, 0.0, k2);
	    if (nt1 <= sdate && nt2 > sdate)
	        break;
	    nt1 = nt2;
	    k1.val = k2.val;
	    }
	phases[0] = truephase(k1.val, 0.0);
	which[0] = 0.0;
	phases[1] = truephase(k1.val, 0.25);
	which[1] = 0.25;
	if ( phases[1] <= sdate )
	    {
	    phases[0] = phases[1];
	    which[0] = which[1];
	    phases[1] = truephase(k1.val, 0.5);
	    which[1] = 0.5;
	    if ( phases[1] <= sdate )
		{
	        phases[0] = phases[1];
	        which[0] = which[1];
	        phases[1] = truephase(k1.val, 0.75);
	        which[1] = 0.75;
	        if ( phases[1] <= sdate )
		    {
		    phases[0] = phases[1];
		    which[0] = which[1];
		    phases[1] = truephase(k2.val, 0.0);
		    which[1] = 0.0;
		    }
		}
	    }
	}

    // kepler - solve the equation of Kepler
    private static double kepler( double m,  double ecc )
	{
	double e, delta;

	e = m = torad(m);
	do
	    {
	    delta = e - ecc * Math.sin(e) - m;
	    e -= delta / (1 - ecc * Math.cos(e));
	    }
	while (Math.abs(delta) > EPSILON);
	return e;
	}

    /// Calculate phase of moon as a fraction.
    // <P>
    // @param pdate time for which the phase is requested, as from jtime()
    // @param pphaseR Ref for illuminated fraction of Moon's disk
    // @param mageR Ref for age of moon in days
    // @param distR Ref for distance in km from center of Earth
    // @param angdiaR Ref for angular diameter in degrees as seen from Earth
    // @param sudistR Ref for distance in km to Sun
    // @param suangdiaR Ref for Sun's angular diameter
    // @return terminator phase angle as a fraction of a full circle (i.e., 0 to 1)
    // 
    public static double phase(
	double pdate, RefDouble pphaseR, RefDouble mageR, RefDouble distR,
	RefDouble angdiaR, RefDouble sudistR, RefDouble suangdiaR )
	{
	double Day, N, M, Ec, Lambdasun, ml, MM, Ev, Ae, A3, MmP,
	       mEc, A4, lP, V, lPP,
	       MoonAge, MoonPhase,
	       MoonDist, MoonDFrac, MoonAng,
	       F, SunDist, 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);

	// 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;

	// 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;

	pphaseR.val = MoonPhase;
	mageR.val = synmonth * (fixangle(MoonAge) / 360.0);
	distR.val = MoonDist;
	angdiaR.val = MoonAng;
	sudistR.val = SunDist;
	suangdiaR.val = SunAng;
	return torad(fixangle(MoonAge));
	}

    }