// JavaScript by Peter Hayes http://www.aphayes.pwp.blueyonder.co.uk/
// Copyright 2001-2004
// This code is made freely available but please keep this notice.
// I accept no liability for any errors in my coding but please
// let me know of any errors you find. My address is on my home page.

// Functions for the planets

// The planet object
//
// Revised by CJL - Oct 2007
// Now includes solar data

function place(name,latitude,longitude) {
  this.name      = name;
  this.latitude  = latitude;
  this.longitude = longitude;
}

function observatory(place,year,month,day,hr,mins,sec,timez) {
  this.name = place.name;
  this.year = year;
  this.month = month;
  this.day = day;
  this.hours = hr;
  this.minutes = mins;
  this.seconds = sec;
  this.tz = timez;
  this.latitude = place.latitude;
  this.longitude = place.longitude;
}

function planet(name,L,a,e,i,N,P) {
  this.name=name;
  this.L=L;
  this.a=a;
  this.e=e;
  this.i=i;
  this.N=N;
  this.P=P;
}

var planets=new Array();
planets[0]=new planet("Mercury",
   new Array(252.250906, 149474.0722491,    0.00030397,    0.000000018),
   new Array(  0.387098310,   0.0,          0.0,           0.0),
   new Array(  0.20563175,    0.000020406, -0.0000000284, -0.00000000017),
   new Array(  7.004986,      0.0018215,   -0.00001809,    0.000000053),
   new Array( 48.330893,      1.1861890,    0.00017587,    0.000000211),
   new Array( 77.456119,      1.5564775,    0.00029589,    0.000000056));

planets[1]=new planet("Venus",
   new Array(181.979801,  58519.2130302,    0.00031060,    0.000000015),
   new Array(  0.723329820,   0.0,          0.0,           0.0),
   new Array(  0.00677188,   -0.000047766,  0.0000000975,  0.00000000044),
   new Array(  3.394662,      0.0010037,   -0.00000088,   -0.000000007),
   new Array( 76.679920,      0.9011190,    0.00040665,   -0.000000080),
   new Array(131.563707,      1.4022188,   -0.00107337,   -0.000005315));

planets[2]=new planet("Earth",
   new Array(100.466449,  36000.7698231,    0.00030368,    0.000000021),
   new Array(  1.000001018,   0.0,          0.0,           0.0),
   new Array(  0.01670862,   -0.000042037, -0.0000001236,  0.00000000004),
   new Array(  0.0,           0.0,          0.0,           0.0),
   new Array(  0.0,           0.0,          0.0,           0.0),
   new Array(102.937348,      1.7195269,    0.00045962,    0.000000499));

planets[3]=new planet("Mars",
   new Array(355.433275,  19141.6964746,    0.00031097,    0.000000015),
   new Array(  1.523679342,   0.0,          0.0,           0.0),
   new Array(  0.09340062,    0.000090483, -0.0000000806, -0.00000000035),
   new Array(  1.849726,     -0.0006010,    0.00001276,   -0.000000006),
   new Array( 49.558093,      0.7720923,    0.00001605,    0.000002325),
   new Array(336.060234,      1.8410331,    0.00013515,    0.000000318));

planets[4]=new planet("Jupiter",
   new Array( 34.351484,   3036.3027889,    0.00022374,    0.000000025),
   new Array(  5.202603191,   0.0000001913, 0.0,           0.0),
   new Array(  0.04849485,    0.000163244, -0.0000004719, -0.00000000197),
   new Array(  1.303270,     -0.0054966,    0.00000465,   -0.000000004),
   new Array(100.464441,      1.0209550,    0.00040117,    0.000000569),
   new Array( 14.331309,      1.6126668,    0.00103127,   -0.000004569));

planets[5]=new planet("Saturn",
   new Array( 50.077471,   1223.5110141,    0.00051952,   -0.000000003),
   new Array(  9.554909596,  -0.0000021389, 0.0,           0.0),
   new Array(  0.05550862,   -0.000346818, -0.0000006456,  0.00000000338),
   new Array(  2.488878,     -0.0037363,   -0.00001516,    0.000000089),
   new Array(113.665524,      0.8770979,   -0.00012067,   -0.000002380),
   new Array( 93.056787,      1.9637694,    0.00083757,    0.000004899));

planets[6]=new planet("Uranus",
   new Array(314.055005,    429.8640561,    0.00030434,    0.000000026),
   new Array( 19.218446062,  -0.0000000372, 0.00000000098, 0.0),
   new Array(  0.04629590,   -0.000027337,  0.0000000790,  0.00000000025),
   new Array(  0.773196,      0.0007744,    0.00003749,   -0.000000092),
   new Array( 74.005947,      0.5211258,    0.00133982,    0.000018516),
   new Array(173.005159,      1.4863784,    0.00021450,    0.000000433));

planets[7]=new planet("Neptune",
   new Array(304.348665,    219.8833092,    0.00030926,    0.000000018),
   new Array( 30.110386869,  -0.0000001663, 0.00000000069, 0.0),
   new Array(  0.00898809,    0.000006408, -0.0000000008, -0.00000000005),
   new Array(  1.769952,     -0.0093082,   -0.00000708,    0.000000028),
   new Array(131.784057,      1.1022057,    0.00026006,   -0.000000636),
   new Array( 48.123691,      1.4262677,    0.00037918,   -0.000000003));

planets[8]=new planet("Sun",
   new Array( 0.0, 0.0, 0.0, 0.0),
   new Array( 0.0, 0.0, 0.0, 0.0),
   new Array( 0.0, 0.0, 0.0, 0.0),
   new Array( 0.0, 0.0, 0.0, 0.0),
   new Array( 0.0, 0.0, 0.0, 0.0),
   new Array( 0.0, 0.0, 0.0, 0.0));

// heliocentric xyz for planet p and observer obs

function helios(p,obs) {
  var T=(jd(obs)-2451545.0)/36525;
  var T2=T*T;
  var T3=T2*T;
  // longitude of ascending node
  var N=rev(p.N[0]+p.N[1]*T+p.N[2]*T2+p.N[3]*T3);
  // inclination
  var i=p.i[0]+p.i[1]*T+p.i[2]*T2+p.i[3]*T3;
  // Mean longitude
  var L=rev(p.L[0]+p.L[1]*T+p.L[2]*T2+p.L[3]*T3);
  // semimajor axis
  var a=p.a[0]+p.a[1]*T+p.a[2]*T2+p.a[3]*T3;
  // eccentricity
  var e=p.e[0]+p.e[1]*T+p.e[2]*T2+p.e[3]*T3;
  // longitude of perihelion
  var P=rev(p.P[0]+p.P[1]*T+p.P[2]*T2+p.P[3]*T3);
  var M=rev(L-P);
  var w=rev(L-N-M);
  // Eccentric anomaly
  var E0=M+(180.0/Math.PI)*e*sind(M)*(1+e*cosd(M));
  var E=E0-(E0-(180.0/Math.PI)*e*sind(E0)-M)/(1-e*cosd(E0));
  while (Math.abs(E0-E) > 0.0005) {
    E0=E;
    E=E0-(E0-(180.0/Math.PI)*e*sind(E0)-M)/(1-e*cosd(E0));
  };
  var x=a*(cosd(E)-e);
  var y=a*Math.sqrt(1-e*e)*sind(E);
  var r=Math.sqrt(x*x+y*y);
  var v=rev(atan2d(y,x));
  // Heliocentric Ecliptic Rectangular Coordinates
  var xeclip=r*(cosd(N)*cosd(v+w)-sind(N)*sind(v+w)*cosd(i));
  var yeclip=r*(sind(N)*cosd(v+w)+cosd(N)*sind(v+w)*cosd(i));
  var zeclip=r*sind(v+w)*sind(i);
  return new Array(xeclip,yeclip,zeclip);
}

// radecr returns ra, dec and earth distance
// obj and base are Heliocentric Ecliptic Rectangular Coordinates
// for the object and earth and obs is the observer

function radecr(obj,base,obs) {
  // Equatorial co-ordinates
  var x=obj[0];
  var y=obj[1];
  var z=obj[2];
  // julian date
  var jdobs=jd(obs);
  // Obliquity of Ecliptic
  var obl=23.4393-3.563E-7*(jdobs-2451543.5);
  // Convert to Geocentric co-ordinates
  var x1=x-base[0];
  var y1=(y-base[1])*cosd(obl)-(z-base[2])*sind(obl);
  var z1=(y-base[1])*sind(obl)+(z-base[2])*cosd(obl);
  // RA and dec
  var ra=rev(atan2d(y1,x1))/15.0;
  var dec=atan2d(z1,Math.sqrt(x1*x1+y1*y1));
  // Earth distance
  var r=Math.sqrt(x1*x1+y1*y1+z1*z1);
  return new Array(ra,dec,r);
}

// Function "posPlan" calculates the current ra,dec,alt,az for a given planet
// [0] Mercury to [7] Neptune
//
// Modified by CJL - 12/10/07

function posPlan(obs,pnum) {
  var obscopy=new Object(); 
  for (var i in obs) obscopy[i] = obs[i];
// Get postions of Earth and given planet
  var earth_xyz = helios(planets[2],obscopy);
  if (pnum==8) {
	  var planet_xyz = new Array(0.0,0.0,0.0);
  } else {
	  var planet_xyz = helios(planets[pnum],obscopy);
  }
// Calculate Ra,Dec,Alt & Azimuth
  var radec=radecr(planet_xyz,earth_xyz,obscopy);
  var altaz=radtoaa(radec[0],radec[1],obscopy);

  return new Array(planets[pnum].name,radec[0],radec[1],altaz[0],altaz[1],radec[2]);
}

// Function "risetPlan" calculates the rise, transit and set times for a given planet on that day
// [0] Mercury to [7] Neptune
//
// Modified by CJL - 12/10/07

function risetPlan(obs,pnum) {
  var obscopy=new Object(); 
  for (var i in obs) obscopy[i] = obs[i];
  // Get postions of Earth and given planet
  var earth_xyz = helios(planets[2],obscopy);
  obscopy.hours=12;
  obscopy.minutes=0;
  obscopy.seconds=0;
  if (pnum==8) {
	  var planet_xyz = new Array(0.0,0.0,0.0);
  } else {
	  var planet_xyz = helios(planets[pnum],obscopy);
  }  
  var lst=p_local_sidereal(obscopy);
  var radec=radecr(planet_xyz,earth_xyz,obscopy);
  var ra=radec[0];
  var dec=radec[1];
  var UTplanet=12.0+ra-lst;

  if (UTplanet < 0.0) UTplanet+=24.0;
  if (UTplanet > 24.0) UTplanet-=24.0;
  // refraction correction 0.583
  var cosLHA=(sind(-0.583)-sind(obs.latitude)*sind(dec)) / (cosd(obs.latitude)*cosd(dec));
  lha=acosd(cosLHA)/15.04107;
  if ((UTplanet-lha) < 0.0) {
    var rise=(24.0+UTplanet-lha);
  } else {
	var rise=(UTplanet-lha);
  }
  var transit=UTplanet;
  if ((UTplanet+lha) > 24.0) {
    var sets=(UTplanet+lha-24.0);
  } else {
    var sets=(UTplanet+lha);
  }

  return new Array(planets[pnum].name,rise,transit,sets);
}

function sitename() {
  var sname=observer.name;
  var latd=Math.abs(observer.latitude);
  var latdi=Math.floor(latd);
  sname+=((latdi < 10) ? " 0" : " ") + latdi;
  latm=60*(latd-latdi); latmi=Math.floor(latm);
  sname+=((latmi < 10) ? ":0" : ":") + latmi;
  lats=60*(latm-latmi); latsi=Math.floor(lats);
  sname+=((latsi < 10) ? ":0" : ":") + latsi;
  sname+=((observer.latitude >= 0) ? "N " : "S ");
  var longd=Math.abs(observer.longitude);
  var longdi=Math.floor(longd);
  sname+=((longdi < 10) ? "0" : "") + longdi;
  longm=60*(longd-longdi); longmi=Math.floor(longm);
  sname+=((longmi < 10) ? ":0" : ":") + longmi;
  longs=60*(longm-longmi); longsi=Math.floor(longs);
  sname+=((longsi < 10) ? ":0" : ":") + longsi;
  sname+=((observer.longitude >= 0) ? "W" : "E");
  return sname;
}

// The Julian date at observer time

function jd(obs) {
  var j = jd0(obs.year,obs.month,obs.day);
  j+=(obs.hours+((obs.minutes+obs.tz)/60.0)+(obs.seconds/3600.0))/24;
  return j;
}

// sidereal time in hours for observer

function p_local_sidereal(obs) {
  var res=g_sidereal(obs.year,obs.month,obs.day);
  res+=1.00273790935*(obs.hours+(obs.minutes+obs.tz+(obs.seconds/60.0))/60.0);
  res-=obs.longitude/15.0;
  while (res < 0) res+=24.0;
  while (res > 24) res-=24.0;
  return res;
}

// sidereal time in hours for Greenwich

function g_sidereal(year,month,day) {
  var T=(jd0(year,month,day)-2451545.0)/36525;
  var res=100.46061837+T*(36000.770053608+T*(0.000387933-T/38710000.0));
  return rev(res)/15.0;
}

// radtoaa converts ra and dec to altitude and azimuth

function radtoaa(ra,dec,obs) {
  var lst=p_local_sidereal(obs);
  var x=cosd(15.0*(lst-ra))*cosd(dec);
  var y=sind(15.0*(lst-ra))*cosd(dec);
  var z=sind(dec);
  // rotate so z is the local zenith
  var xhor=x*sind(obs.latitude)-z*cosd(obs.latitude);
  var yhor=y;
  var zhor=x*cosd(obs.latitude)+z*sind(obs.latitude);
  var azimuth=rev(atan2d(yhor,xhor)+180.0); // so 0 degrees is north
  var altitude=atan2d(zhor,Math.sqrt(xhor*xhor+yhor*yhor));
  return new Array(altitude,azimuth);
}

// Function to calculate angular separation given 2 RA and Decs

function AngSep (ra1,dec1,ra2,dec2) {
  var a_ra1=(ra1/24)*360; var a_ra2=(ra2/24)*360;
  var angle=acosd((cosd(90-dec1)*cosd(90-dec2))+(sind(90-dec1)*sind(90-dec2)*cosd(a_ra1-a_ra2)));
  return angle;
}

// datestring returns the datestring in form yyyy:mm:dd

function datestring(obs) {
  var datestr = "";  datestr += obs.year;
  datestr += ((obs.month < 10) ? ":0" : ":") + obs.month;
  datestr += ((obs.day < 10) ? ":0" : ":") + obs.day;
  return datestr;
}

// Comma used to seperate large numbers out with commas

function comma(SS) { var T='', S=String(SS), L=S.length-1, C, j
  for (j=0; j<=L; j++) {
    T+=C=S.charAt(j)
    if ((j < L) && ((L-j)%3 == 0) && (C != '-')) T+=',' }
  return T }

// Compass returns the compass direction of given angle

function compass(angle) {
  if ((angle>348.75) && (angle<=11.25)) var heading="North";
  if ((angle>11.25) && (angle<=33.75)) var heading="NNE";
  if ((angle>33.75) && (angle<=56.25)) var heading="NE";
  if ((angle>56.25) && (angle<=78.75)) var heading="ENE";
  if ((angle>78.75) && (angle<=101.25)) var heading="East";
  if ((angle>101.25) && (angle<=123.75)) var heading="ESE";
  if ((angle>123.75) && (angle<=146.25)) var heading="SE";
  if ((angle>146.25) && (angle<=168.75)) var heading="SSE";
  if ((angle>168.75) && (angle<=191.25)) var heading="South";
  if ((angle>191.25) && (angle<=213.75)) var heading="SSW";
  if ((angle>213.75) && (angle<=236.25)) var heading="SW";
  if ((angle>236.25) && (angle<=258.75)) var heading="WSW";
  if ((angle>258.75) && (angle<=281.25)) var heading="West";
  if ((angle>281.25) && (angle<=303.75)) var heading="WNW";
  if ((angle>303.75) && (angle<=326.25)) var heading="NW";
  if ((angle>326.25) && (angle<=348.75)) var heading="NNW";
  return (heading);
}

function PlanEvents(year,month) {
  var atlas=new Array(new place("Liverpool",53.4,3.0));
  var obs=new observatory(atlas[0],year,month,1,12,0,0,0);
  var jstart=julian(year,month,1,0);
  var weekday=new Array(7);  weekday[1]="Monday" ; weekday[2]="Tuesday" ; weekday[3]="Wednesday" ; 
  weekday[4]="Thursday" ; weekday[5]="Friday" ; weekday[6]="Saturday" ; weekday[0]="Sunday" ;
  if ((month==4) || (month==6) || (month==9) || (month==11)) {
	  var days=30;
  } else if (month==2) {
	  var days=28;
  } else {
	  var days=31;
  }
  var jstop=julian(year,month,days,24);
  var ticktock=Math.round((jstop-jstart)*24);
  var ev=0; var events=new Array(10);
// Loop through the planets in turn
  for (pcount=0; pcount < 8; pcount++) {
    if (pcount==2) {
	  continue;
	} else {
// Now loop through the month in hourly steps
      var databin=new Array(ticktock); var minvalue=new Array(270,0,0,0); var maxvalue=new Array(0,0,0.0);
      for (count = 0; count < ticktock; count++) {
        var time=jstart+(count/24);
	    var dateres=jdtocd(time); 
        obs.day=dateres[2]; obs.hours=dateres[4];
        var plan_one=posPlan(obs,pcount); var plan_two=posPlan(obs,8);
// Fill up 'databin' array with hourly angular separations
        databin[count]=AngSep(plan_one[1],plan_one[2],plan_two[1],plan_two[2]);
        if (databin[count]<minvalue[0]) {
			minvalue[0]=databin[count]; minvalue[1]=obs.day; minvalue[2]=obs.hours;
		}
        if (databin[count]>maxvalue[0]) {
			maxvalue[0]=databin[count]; maxvalue[1]=obs.day; maxvalue[2]=obs.hours;
		}
	  }
// Test for Conjunction - Inferior or Superior
      if ((minvalue[0]<databin[0]) && (minvalue[0]<databin[count-1])) {
         if ((plan_one[5]>1) && (pcount<2)) var aspect="Superior Conjunction";
		 if ((plan_one[5]<1) && (pcount<2)) var aspect="Inferior Conjunction";
		 if ((plan_one[5]>1) && (pcount>2)) var aspect="Conjunction";
         events[ev]=new Array(plan_one[0],aspect,minvalue[1],minvalue[2],minvalue[0]);
		 ev++; // Increase event counter
	  }
// Test for Opposition and Elongation
      if ((maxvalue[0]>databin[0]) && (maxvalue[0]>databin[count-1])) {
         if (pcount<2) var aspect="Greatest Elongation";
		 if (pcount>2) var aspect="Opposition";
         events[ev]=new Array(plan_one[0],aspect,maxvalue[1],maxvalue[2],maxvalue[0]);
		 ev++; // Increase event counter
	  }
    }
  }
  return new Array(events,ev);
}