// Functions for comets

// Copyright Ole Nielsen 2002-2007
// Please read copyright notice in astrotools2.html source

// Formulae from Paul Schlyter's article "Computing planetary positions" available at 
// http://hem.passagen.se/pausch/comp/ppcomp.html

// Orbital elements in Skymap format, retrieve fresh elements from 
// http://cfa-www.harvard.edu/cfa/ps/Ephemerides/Comets/SoftwareComets.html
// All elements must be referred to J2000.0 
// Paste desired lines from the Skymap file into the variable assignment below, 
// start line with a " then the pasted line and append ", to end of line.
// use the template line just above the var declaration for aiding alignment of columns
// Fields: name: 0-46; Ty: 47-51; Tm: 52-53; Td: 55-61; q: 64-71; e: 79-86; wL 88-95; N: 97-104; i: 106-113; G: 116-119; H: 122-125

// namenamenamenamenamenamenamename             yyyy mm dd.dddd  q.qqqqqq       e.eeeeee www.wwww NNN.NNNN iii.iiii  GG.G  HH.H
var c_elements = [
"29P Schwassmann-Wachmann                       2004 06 25.9245  5.720326       0.044383  47.8992 312.6523   9.3916   4.0   4.0",
"C/2006 M4 SWAN                                 2006 09 28.7276  0.783144       1.000439  62.5838 148.7273 111.8193   8.5   4.0",
"C/2006 L1 Garradd                              2006 10 17.9998  1.462068       0.997421 338.4078 101.7622 143.2432  12.0   4.0",
"4P Faye                                        2006 11 15.4569  1.667345       0.566680 205.0172 199.3082   9.0316   8.0   6.0",
"C/2006 P1 McNaught                             2007 01 12.7983  0.170741       1.000011 155.9780 267.4146  77.8371  10.0   4.0",
"P/2001 Q2 Petriew                              2007 02 24.6172  0.937604       0.698138 181.9227 214.1023  13.9745  11.0   4.0",
"96P Machholz                                   2007 04 04.6194  0.124618       0.958684  14.6181  94.5507  59.9553  13.0   4.8",
"2P Encke                                       2007 04 19.3117  0.339269       0.847040 186.5230 334.5714  11.7543  11.5   6.0",
""];

// Help for elements retrieved from JPL (http://ssd.jpl.nasa.gov/sbdb.cgi):
// q, e and i are the same
// T = tp (in ymd format)
// w = peri
// N = node
// G = M1
// H = K1/2.5

var c_index = -1;	// currently selected comet

function comet(name,Ty,Tm,Td,q,e,w,N,i,G,H) {
// The comet object
	this.name=name;	// IAU designation
	this.Ty=Ty;	// year of perihelion
	this.Tm=Tm;	// month of perihelion
	this.Td=Td;	// date of perihelion
	this.q=q;	// perihelion distance
	this.e=e;	// eccentricity
	this.w=w;	// argument of perihelion (= lowercase omega)
	this.N=N; 	// longitude of ascending node (= uppercase omega)
	this.i=i;	// inclination
	this.G=G; 	// absolute magnitude
	this.H=H;	// slope parameter
}

var comets=new Array();	// orbital elements in internal format


function add_comets() {
// convert element text format to internal format
	var elm = c_elements;
	for (var i=0; i < elm.length-1; i++) {
		comets[i] = new comet(elm[i].substring(0,31),
			parseInt(elm[i].substring(47,51),10),
			parseInt(elm[i].substring(52,54),10),
			parseFloat(elm[i].substring((elm[i].charAt(55)=="0"?56:55),62)),
			parseFloat(elm[i].substring(64,72)),
			parseFloat(elm[i].substring(79,87)),
			parseFloat(elm[i].substring(88,96)),
			parseFloat(elm[i].substring(97,105)),
			parseFloat(elm[i].substring(106,114)),
			parseFloat(elm[i].substring(116,120)),
			parseFloat(elm[i].substring(122,126)));
	}
}	// add_comets()


function comet_xyz(cn,jday) {
// heliocentric xyz for comet (cn is index to comets)
// based on Paul Schlyter's page http://www.stjarnhimlen.se/comp/ppcomp.html
// returns heliocentric x, y, z, distance, longitude and latitude of object
	var d = jday-2451543.5;
	var cm = comets[cn];
	var q = cm.q;
	var e = cm.e; 
	var Tj = jd0(cm.Ty,cm.Tm,cm.Td);	// get julian day of perihelion time
	if (e > 0.98)	{ 
		// treat as near parabolic (approx. method valid inside orbit of Pluto)
		var k = 0.01720209895;	// Gaussian gravitational constant
		var a = 0.75 * (jday-Tj) * k * Math.sqrt( (1 + e) / (q*q*q) ); 
		var b = Math.sqrt( 1 + a*a );
		var W = cbrt(b + a) - cbrt(b - a);
		var c = 1 + 1/(W*W);
		var f = (1 - e) / (1 + e);
		var g = f / (c*c);
		var a1 = (2/3) + (2/5) * W*W;
		var a2 = (7/5) + (33/35) * W*W + (37/175) * W*W*W*W;
		var a3 = W*W * ( (432/175) + (956/1125) * W*W + (84/1575) * W*W*W*W);
		var w = W * ( 1 + g * c * ( a1 + a2*g + a3*g*g ) );
		var v = 2 * atand(w);
		var r = q * ( 1 + w*w ) / ( 1 + w*w * f );
	}
	else {		// treat as elliptic
		var a = q / (1.0 - e);
		var P = 365.2568984 * Math.sqrt(a*a*a);	// period in days
		var M = 360.0 * (jday-Tj)/P; 	// mean anomaly
		// eccentric anomaly E
		var E0 = M + RAD2DEG*e*sind(M) * ( 1.0+e*cosd(M) );
		var E1 = E0 - ( E0-RAD2DEG*e*sind(E0)-M ) / ( 1.0-e*cosd(E0) );
		while (Math.abs(E0-E1) > 0.0005) {
			E0 = E1;
			E1 = E0 - ( E0 - RAD2DEG*e*sind(E0)-M ) / ( 1.0-e*cosd(E0) );
		}
		var xv = a*(cosd(E1) - e);
		var yv = a*Math.sqrt(1.0 - e*e) * sind(E1);
		var v = rev(atan2d( yv, xv ));		// true anomaly
		var r = Math.sqrt( xv*xv + yv*yv );	// distance
	}	// from here common for all orbits
	var N = cm.N + 3.82394E-5 * d;	// precess from J2000.0 to now;
	var w = cm.w;	// why not precess this value?
	var i = cm.i;
	var xh = r * ( cosd(N)*cosd(v+w) - sind(N)*sind(v+w)*cosd(i) );
	var yh = r * ( sind(N)*cosd(v+w) + cosd(N)*sind(v+w)*cosd(i) );
	var zh = r * ( sind(v+w)*sind(i) );
	var lonecl = atan2d(yh, xh);
	var latecl = atan2d(zh, Math.sqrt(xh*xh + yh*yh + zh*zh));
	return new Array(xh,yh,zh,r,lonecl,latecl);
}	// comet_xyz()


function CometAlt(jday,obs) {
// Alt/Az, hour angle, ra/dec, ecliptic long. and lat, illuminated fraction (=1.0), dist(Sun), dist(Earth), brightness of planet p
	var cn = c_index;
	var sun_xyz = sunxyz(jday);
	var planet_xyz = comet_xyz(cn,jday);
	var dx = planet_xyz[0]+sun_xyz[0];
	var dy = planet_xyz[1]+sun_xyz[1];
	var dz = planet_xyz[2]+sun_xyz[2];
	var lon = rev( atan2d(dy, dx) );
	var lat = atan2d(dz, Math.sqrt(dx*dx+dy*dy));
	var radec = radecr(planet_xyz, sun_xyz, jday, obs);
	var ra = radec[0]; 
	var dec = radec[1];
	var altaz = radec2aa(ra, dec, jday, obs);
	var dist = radec[2];
	var r = planet_xyz[3];
	var mag = comets[cn].G  + 5*log10(dist) + 2.5*comets[cn].H*log10(r);	// Schlyter's formula is wrong!
	return new Array (altaz[0], altaz[1], altaz[2], ra, dec, lon, lat, 1.0, r, dist, mag); 
}	// CometAlt()