/*******************************************************************************
 * Copyright (c) 2010 Attila Kovacs <attila_kovacs@post.harvard.edu>.
 * All rights reserved. 
 * 
 * This file is part of crush.
 * 
 *     crush 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 3 of the License, or
 *     (at your option) any later version.
 * 
 *     crush 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 crush.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * Contributors:
 *     Attila Kovacs <attila_kovacs@post.harvard.edu> - initial API and implementation
 ******************************************************************************/

	
// TODO precess with proper motion...
// TODO nutation correction...


package util;

public class Precession {
	private CoordinateEpoch fromEpoch, toEpoch;
	private double[][] P;
	private double[] v = new double[3], l = new double[3];
	
	public Precession(double fromJulianEpoch, double toJulianEpoch) {
		this(new JulianEpoch(fromJulianEpoch), new JulianEpoch(toJulianEpoch));
	}
	
	public Precession(CoordinateEpoch from, CoordinateEpoch to) {
		fromEpoch = from;
		toEpoch = to;
		if(fromEpoch.equals(toEpoch)) P = null;
		else getMatrix();
	}

	//  Precession from Lederle & Schwan, Astronomy and Astrophysics, 134, 1-6 (1984)
	private void getMatrix() {
		double fromJulianYear = fromEpoch.getJulianYear();
		double toJulianYear = toEpoch.getJulianYear();
		
		final double tau = (fromJulianYear - 2000.0) * Unit.year / Unit.julianCentury;
		final double t = (toJulianYear - fromJulianYear) * Unit.year / Unit.julianCentury;


		double eta = (2305.6997 + (1.39744 + 0.000060 * tau) * tau 
				+ (0.30201 - 0.000270 * tau + 0.017996 * t) * t) * t * Unit.arcsec;

		double z = (2305.6997 + (1.39744 + 0.000060 * tau) * tau 
				+ (1.09543 + 0.000390 * tau + 0.018326 * t) * t) * t * Unit.arcsec;

		double theta = (2003.8746 - (0.85405 + 0.000370 * tau) * tau
				- (0.42707 + 0.000370 * tau + 0.041803 * t) * t) * t * Unit.arcsec;	

		P = R3(-z).dot(R2(theta)).dot(R3(-eta)).value;	
	}
	
	
	public void precess(EquatorialCoordinates equatorial) {		
		if(P == null) return;
		
		v[0] = equatorial.cosLat * Math.cos(equatorial.RA()); 
		v[1] = equatorial.cosLat * Math.sin(equatorial.RA());
		v[2] = equatorial.sinLat;

		for(int i=0; i<l.length; i++) {
			double sum = 0.0;
			final double[] Pi = P[i];
			for(int j=0; j<v.length; j++) sum += Pi[j] * v[j];
			l[i] = sum;
		}
	
		double ml = Math.hypot(l[0], l[1]);

		equatorial.setRA(Math.atan2(l[1], l[0]));
		equatorial.setDEC(Math.atan2(l[2], ml));

		equatorial.epoch = toEpoch;
	}

	/*
	private Matrix R1(double phi) {
		double c = Math.cos(phi);
		double s = Math.sin(phi);

		double[][] R = { {  1, 0, 0 }, 
				{  0, c, s },
				{  0, -s, c } };	
		return new Matrix(R);
	}
	*/
	
	private Matrix R2(double phi) {
		double c = Math.cos(phi);
		double s = Math.sin(phi);

		double[][] R = { {  c, 0,-s }, 
				{  0, 1, 0 },
				{  s, 0, c } };	
		return new Matrix(R);
	}

	private Matrix R3(double phi) {
		double c = Math.cos(phi);
		double s = Math.sin(phi);

		double[][] R = { {  c, s, 0 }, 
				{ -s, c, 0 },
				{  0, 0, 1 } };	
		return new Matrix(R);		
	}
	
	
}
