/*  GeochronApplet.java
 *
 *  Applet to display day/night regions of Earth
 *
 *  (c) 2003,2009 Rich Townsend (townsend@astro.wisc.edu)
 *
 *  Distributed under the GNU GPL; see:
 *   http://www.gnu.org/licenses/gpl.html 
 */

import java.util.*;
import java.applet.*;
import java.awt.*;
import java.awt.image.*;
import java.net.*;

public class GeochronApplet extends Applet {

    private static final double DtoR = Math.PI/180.;
    private static final double RtoD = 180./Math.PI;

    private static final double dayAltMin = 0.;    // Minimum solar altitude for daytime
    private static final double nightAltMax = -9.; // Maximum solar altitude for nighttime

    private BufferedImage dayImage;
    private BufferedImage nightImage;

    private BufferedImage blendedImage;

    private int n_x;
    private int n_y;

    public void init() {

	// Set up the image URLs

	URL dayImageURL;
	URL nightImageURL;

	try {
	    dayImageURL = new URL(this.getParameter("dayImageURL"));
	}
	catch (MalformedURLException e1) {
	    this.showStatus("dayImageURL is malformed");
	    return;
	}

	try {
	    nightImageURL = new URL(this.getParameter("nightImageURL"));
	}
	catch (MalformedURLException e2) {
	    this.showStatus("nightImageURL is malformed");
	    return;
	}

	// Load the images

	Image dayImage = this.getImage(dayImageURL);
	Image nightImage = this.getImage(nightImageURL);

	MediaTracker tracker = new MediaTracker(this);

	tracker.addImage(dayImage, 0);
	tracker.addImage(nightImage, 1);

	try {
	    this.showStatus("Loading images...");
	    tracker.waitForAll();
	    this.showStatus("Running...");
	}
	catch (InterruptedException e3) {
	    this.showStatus("Unable to load images!");
	}

	// Store the images

	this.n_x = dayImage.getWidth(this);
	this.n_y = dayImage.getHeight(this);

	this.dayImage = new BufferedImage (this.n_x, this.n_y, BufferedImage.TYPE_INT_RGB);
	this.nightImage = new BufferedImage (this.n_x,this.n_y, BufferedImage.TYPE_INT_RGB);

	this.dayImage.getGraphics().drawImage(dayImage, 0, 0, this);
	this.nightImage.getGraphics().drawImage(nightImage, 0, 0, this);

	// Create space for the blended image

	this.blendedImage = new BufferedImage (this.n_x, this.n_y, BufferedImage.TYPE_INT_RGB);

	// Set up the update timer task

	TimerTask updateTask = new TimerTask() {
		public void run() { 
		    setBlendedImage(); 
		    repaint();
		}
	    };

	// Start the updates

	Timer timer = new Timer();
	timer.schedule(updateTask, 0, 60000);

    }

    public void paint(Graphics g) {
	g.drawImage(blendedImage, 0, 0, Color.black, this);
    }

    public void update(Graphics g) {
	g.drawImage(blendedImage, 0, 0, Color.black, this);
    }

    void setBlendedImage () {

	// Get a calendar for the current GMT

	Calendar calendar = Calendar.getInstance(TimeZone.getTimeZone("GMT"));

	// Calculate the Greenwich Sidereal Time

	double GST = getGreenwichSiderealTime(calendar);

	// Calculate the solar right ascension and declination

	double alpha = getSolarRightAscension(calendar);
	double delta = getSolarDeclination(calendar);

	// Loop over the y-pixels of the night/day images

	int i_x;
	int i_y;

	for(i_x = 0; i_x < this.n_x; i_x++) {

	    // Calculate the longitude

	    double longitude = 180. - (i_x+0.5)/this.n_x*360.;

	    // Calculate the solar hour angle 

	    double HA = GST*360./24. - longitude - alpha;

	    for(i_y = 0; i_y < this.n_y; i_y++) {

		// Calculate the latitude

		double latitude = 90. - (i_y+0.5)/this.n_y*180.;

		// Calculate the altitude of the sun

		double alt = getAltitude(HA, delta, latitude, longitude);

		// Work out the interpolation factor for drawing the pixel

		double intFactor;

		if(alt > dayAltMin) {
		    intFactor = 1.;
		}
		else if(alt < nightAltMax) {
		    intFactor = 0.;
		}
		else {
		    intFactor = (alt - nightAltMax)/(dayAltMin - nightAltMax);
		}

		// Get the RGB pixels of the day and night images

		int dayRGB = this.dayImage.getRGB(i_x, i_y);
		int nightRGB = this.nightImage.getRGB(i_x, i_y);

		int dayR = dayRGB >>> 16 & 0xFF;
		int dayG = dayRGB >>> 8 & 0xFF;
		int dayB = dayRGB & 0xFF;

		int nightR = nightRGB >>> 16 & 0xFF;
		int nightG = nightRGB >>> 8 & 0xFF;
		int nightB = nightRGB & 0xFF;

		// Calculate the interpolated value of the blended pixel

		int blendedRGB = (int) (dayR*intFactor + nightR*(1.-intFactor)) << 16 |
		    (int) (dayG*intFactor + nightG*(1.-intFactor)) << 8 |
		    (int) (dayB*intFactor + nightB*(1.-intFactor));

		this.blendedImage.setRGB(i_x, i_y, blendedRGB);

	    }
	}

    }

    private double getDaysSinceJ2000(Calendar calendar) {

	// Calculate the number of days from the epoch J2000.0
        // the specified date

	int year = calendar.get(Calendar.YEAR);
	int month = calendar.get(Calendar.MONTH)+1;
	int day = calendar.get(Calendar.DAY_OF_MONTH);

	double D0 = 367*year - 7*(year+(month+9)/12)/4 + 275*month/9 + day - 730531.5;

	double D = D0 + getGreenwichMeanTime(calendar)/24.;

	return D;

    }


    private double getGreenwichMeanTime(Calendar calendar) {

	// Get the Greenwich Mean Time, in hours

	int hour = calendar.get(Calendar.HOUR_OF_DAY);
	int minute = calendar.get(Calendar.MINUTE);
	int second = calendar.get(Calendar.SECOND);

	double GMT = hour + minute/60. + second/3600.;

	return GMT;

    }


    private double getGreenwichSiderealTime(Calendar calendar) {

	// Get the number of days since J2000.0

	double D = getDaysSinceJ2000(calendar);

	// Calculate the GST

	double T = D/36525.;

	double GST = (280.46061837 + 360.98564736629*D + 0.000388*T*T)*24./360.;

	// Phase it to within 24 hours

	while(GST < 0.) { GST += 24.; }
	while(GST >= 24.) { GST -= 24.; }

	return GST;

    }


    private double getSolarRightAscension(Calendar calendar) {

	// Calculate the number of days from the epoch J2000.0

	double D = getDaysSinceJ2000(calendar);

	// Convert this into centuries

	double T = D/36525.;

	// Calculate the mean longitude and anomaly

	double L = 279.697 + 36000.769*T;
	double M = 358.476 + 35999.050*T;

	// Calculate the true longitude

	double lambda = L + (1.919 - 0.005*T)*Math.sin(M*DtoR) + 0.020*Math.sin(2*M*DtoR);

	// Calculate the obliquity

	double epsilon = 23.452 - 0.013*T;

	// Calculate the right ascension, in degrees

	double alpha = Math.atan2(Math.sin(lambda*DtoR)*Math.cos(epsilon*DtoR),Math.cos(lambda*DtoR)) * RtoD;

	return alpha;

    }


    private double getSolarDeclination(Calendar calendar) {

	// Calculate the number of days from the epoch J2000.0

	double D = getDaysSinceJ2000(calendar) + getGreenwichMeanTime(calendar)/24.;

	// Convert this into centuries

	double T = D/36525.;

	// Calculate the obliquity

	double epsilon = 23.452 - 0.013*T;

	// Calculate the declination, in degrees

	double delta = Math.asin(Math.sin(getSolarRightAscension(calendar)*DtoR)*Math.sin(epsilon*DtoR)) * RtoD;

	return delta;

    }

    private double getAltitude(double HA, double delta, double latitude, double longitude) {

	// Calculate the altitude, in degrees

	double alt = Math.asin(Math.sin(latitude*DtoR)*Math.sin(delta*DtoR) + 
                               Math.cos(latitude*DtoR)*Math.cos(delta*DtoR)*Math.cos(HA*DtoR))*RtoD;

	return alt;

    }

}