var c1 = new Array(
0,0,1,0,6288774,
2,0,-1,0,1274027,
2,0,0,0,658314,
0,0,2,0,213618,
0,1,0,0,-185116,
0,0,0,2,-114332,
2,0,-2,0,58793,
2,-1,-1,0,57066,
2,0,1,0,53322,
2,-1,0,0,45758,
0,1,-1,0,-40923,
1,0,0,0,-34720,
0,1,1,0,-30383,
2,0,0,-2,15327,
0,0,1,2,-12528,
0,0,1,-2,10980,
4,0,-1,0,10675,
0,0,3,0,10034,
4,0,-2,0,8548,
2,1,-1,0,-7888,
2,1,0,0,-6766,
1,0,-1,0,-5163,
1,1,0,0,4987,
2,-1,1,0,4036,
2,0,2,0,3994,
4,0,0,0,3861,
2,0,-3,0,3665,
0,1,-2,0,-2689,
2,0,-1,2,-2602,
2,-1,-2,0,2390,
1,0,1,0,-2348,
2,-2,0,0,2236,
0,1,2,0,-2120,
0,2,0,0,-2069,
2,-2,-1,0,2048,
2,0,1,-2,-1773,
2,0,0,2,-1595,
4,-1,-1,0,1215,
0,0,2,2,-1110,
3,0,-1,0,-892,
2,1,1,0,-810,
4,-1,-2,0,759,
0,2,-1,0,-713,
2,2,-1,0,-700,
2,1,-2,0,691,
2,-1,0,-2,596,
4,0,1,0,549,
0,0,4,0,537,
4,-1,0,0,520,
1,0,-2,0,-487,
2,1,0,-2,-399,
0,0,2,-2,-381,
1,1,1,0,351,
3,0,-2,0,-340,
4,0,-3,0,330,
2,-1,2,0,327,
0,2,1,0,-323,
1,1,-1,0,299,
2,0,3,0,294);

var c2 = new Array(
0,0,1,0,-20905355,
2,0,-1,0,-3699111,
2,0,0,0,-2955968,
0,0,2,0,-569925,
0,1,0,0,48888,
0,0,0,2,-3149,
2,0,-2,0,246158,
2,-1,-1,0,-152138,
2,0,1,0,-170733,
2,-1,0,0,-204586,
0,1,-1,0,-129620,
1,0,0,0,108743,
0,1,1,0,104755,
2,0,0,-2,10321,
0,0,1,-2,79661,
4,0,-1,0,-34782,
0,0,3,0,-23210,
4,0,-2,0,-21636,
2,1,-1,0,24208,
2,1,0,0,30824,
1,0,-1,0,-8379,
1,1,0,0,-16675,
2,-1,1,0,-12831,
2,0,2,0,-10445,
4,0,0,0,-11650,
2,0,-3,0,14403,
0,1,-2,0,-7003,
2,-1,-2,0,10056,
1,0,1,0,6322,
2,-2,0,0,-9884,
0,1,2,0,5751,
2,-2,-1,0,-4950,
2,0,1,-2,4130,
4,-1,-1,0,-3958,
3,0,-1,0,3258,
2,1,1,0,2616,
4,-1,-2,0,-1897,
0,2,-1,0,-2117,
2,2,-1,0,2354,
4,0,1,0,-1423,
0,0,4,0,-1117,
4,-1,0,0,-1571,
1,0,-2,0,-1739,
0,0,2,-2,-4421,
0,2,1,0,1165,
2,0,-1,-2,8752);

var c3 = new Array(
0,0,0,1,5128122,
0,0,1,1,280602,
0,0,1,-1,277693,
2,0,0,-1,173237,
2,0,-1,1,55413,
2,0,-1,-1,46271,
2,0,0,1,32573,
0,0,2,1,17198,
2,0,1,-1,9266,
0,0,2,-1,8822,
2,-1,0,-1,8216,
2,0,-2,-1,4324,
2,0,1,1,4200,
2,1,0,-1,-3359,
2,-1,-1,1,2463,
2,-1,0,1,2211,
2,-1,-1,-1,2065,
0,1,-1,-1,-1870,
4,0,-1,-1,1828,
0,1,0,1,-1794,
0,0,0,3,-1749,
0,1,-1,1,-1565,
1,0,0,1,-1491,
0,1,1,1,-1475,
0,1,1,-1,-1410,
0,1,0,-1,-1344,
1,0,0,-1,-1335,
0,0,3,1,1107,
4,0,0,-1,1021,
4,0,-1,1,833,
0,0,1,-3,777,
4,0,-2,1,671,
2,0,0,-3,607,
2,0,2,-1,596,
2,-1,1,-1,491,
2,0,-2,1,-451,
0,0,3,-1,439,
2,0,2,1,422,
2,0,-3,-1,421,
2,1,-1,1,-366,
2,1,0,1,-351,
4,0,0,1,331,
2,-1,1,1,315,
2,-2,0,-1,302,
0,0,1,3,-283,
2,1,1,-1,-229,
1,1,0,-1,223,
1,1,0,1,223,
0,1,-2,-1,-220,
2,1,-1,-1,-220,
1,0,1,1,-185,
2,-1,-2,-1,181,
0,1,2,1,-177,
4,0,-2,-1,176,
4,-1,-1,-1,166,
1,0,1,-1,-164,
4,0,1,-1,132,
1,0,-1,-1,-119,
4,-1,0,-1,115,
2,-2,0,1,107);
		
var n1 = c1.length/5;		
var n2 = c2.length/5;		
var n3 = c3.length/5;		

// Alternative version of Sun position based on Schlyter's method

// 'Meeus' means "Astronomical Algorithms", 2nd ed. by Jean Meeus

// ecliptic position of the Sun relative to Earth (basically simplified version of planetxyz calc)
function sunxyz(jday) {
	// return x,y,z ecliptic coordinates, distance, true longitude
	// days counted from 1999 Dec 31.0 UT
	var d=jday-2451543.5;
	var w = 282.9404 + 4.70935E-5 * d;		// argument of perihelion
	var e = 0.016709 - 1.151E-9 * d; 
	var M = rev(356.0470 + 0.9856002585 * d); // mean anomaly
	var E = M + e*RAD2DEG * sind(M) * ( 1.0 + e * cosd(M) );
	var xv = cosd(E) - e;
	var yv = Math.sqrt(1.0 - e*e) * sind(E);
	var v = atan2d( yv, xv );		// true anomaly
	var r = Math.sqrt( xv*xv + yv*yv );	// distance
	var lonsun = rev(v + w);	// true longitude
	var xs = r * cosd(lonsun);		// rectangular coordinates, zs = 0 for sun 
	var ys = r * sind(lonsun);
	return new Array(xs,ys,0,r,lonsun,0);
}

function SunAlt(jday,obs) {
// return alt, az, time angle, ra, dec, ecl. long. and lat=0, illum=1, 0, dist, brightness 
	var sdat=sunxyz(jday);
	var ecl = 23.4393 - 3.563E-7 * (jday-2451543.5); 
	var xe = sdat[0];
	var ye = sdat[1]*cosd(ecl);
	var ze = sdat[1]*sind(ecl);
	var ra = rev(atan2d(ye,xe));
	var dec = atan2d(ze, Math.sqrt(xe*xe+ye*ye));
	var topo=radec2aa(ra,dec,jday,obs);
	return new Array( topo[0] , topo[1] , topo[2] , ra , dec , sdat[4], 0, 1, 0, sdat[3], -26.74 );
}

// Sun rise and set times (if twilight==-0.833) or desired twilight time. Return julian days
function sunrise(obs,twilight) {
	// obs is a reference variable make a copy
	var obscopy=new Object();
	for (var i in obs) obscopy[i] = obs[i];
	obscopy.hours=12;
	obscopy.minutes=0;
	obscopy.seconds=0;
	var riseset=new Array(0.0, 0.0, false, 0.0, 0.0);
	var lst=local_sidereal(obscopy);
	var jday=jdnoDelta(obscopy);
	var radec = SunAlt(jday, obscopy);
	var ra = radec[3]; var dec = radec[4];
	var UTsun=12.0+ra/15.0-lst;
	if (UTsun < 0.0) UTsun+=24.0;
	if (UTsun > 24.0) UTsun-=24.0;
	var cosLHA = (sind(twilight)-sind(obs.latitude)*sind(dec)) / (cosd(obs.latitude)*cosd(dec));
	// Check for midnight sun and midday night. "riseset[2]" false if no rise and set found
	riseset[2]=false;
	if (cosLHA <= 1.0 && cosLHA >= -1.0) {
		// rise/set times allowing for not today.
		riseset[2]=true;
		var lha=acosd(cosLHA)/15.0;
		if ((UTsun-lha) < 0.0) {
			var rtime=24.0+UTsun-lha;
		} else {
			var rtime=UTsun-lha;
		}
		riseset[0]=jday+rtime/24.0-0.5;
		if ((UTsun+lha) > 24.0) {
			var stime=UTsun+lha-24.0;
		} else {
			var stime=UTsun+lha;
			riseset[4]=stime;
		}
		riseset[1]=jday+stime/24.0-0.5;
		// riseset[3] and [4] are times in UT hours
		riseset[3]=rtime;
		riseset[4]=stime;
	}
	return(riseset);
}

function calc_nutazione(JDE){
//"D","M","M1","F","omega"
var c= new Array( 
0,0,0,0,1,-171996,-174.2,92025,8.9,
-2,0,0,2,2,-13187,-1.6,5736,-3.1,
0,0,0,2,2,-2274,-0.2,977,-0.5,
0,0,0,0,2,2062,0.2,-895,0.5,
0,1,0,0,0,1426,-3.4,54,-0.1,
0,0,1,0,0,712,0.1,-7,0,
-2,1,0,2,2,-517,1.2,224,-0.6,
0,0,0,2,1,-386,-0.4,200,0,
0,0,1,2,2,-301,0,129,-0.1,
-2,-1,0,2,2,217,-0.5,-95,0.3,
-2,0,1,0,0,-158,0,0,0,
-2,0,0,2,1,129,0.1,-70,0,
0,0,-1,2,2,123,0,-53,0,
2,0,0,0,0,63,0,0,0,
0,0,1,0,1,63,0.1,-33,0,
2,0,-1,2,2,-59,0,26,0,
0,0,-1,0,1,-58,-0.1,32,0,
0,0,1,2,1,-51,0,27,0,
-2,0,2,0,0,48,0,0,0,
0,0,-2,2,1,46,0,-24,0,
2,0,0,2,2,-38,0,16,0,
0,0,2,2,2,-31,0,13,0,
0,0,2,0,0,29,0,0,0,
-2,0,1,2,2,29,0,-12,0,
0,0,0,2,0,26,0,0,0,
-2,0,0,2,0,-22,0,0,0,
0,0,-1,2,1,21,0,-10,0,
0,2,0,0,0,17,-0.1,0,0,
2,0,-1,0,1,16,0,-8,0,
-2,2,0,2,2,-16,0.1,7,0,
0,1,0,0,1,-15,0,9,0,
-2,0,1,0,1,-13,0,7,0,
0,-1,0,0,1,-12,0,6,0,
0,0,2,-2,0,11,0,0,0,
2,0,-1,2,1,-10,0,5,0,
2,0,1,2,2,-8,0,3,0,
0,1,0,2,2,7,0,-3,0,
-2,1,1,0,0,-7,0,0,0,
0,-1,0,2,2,-7,0,3,0,
2,0,0,2,1,-7,0,3,0,
2,0,1,0,0,6,0,0,0,
-2,0,2,2,2,6,0,-3,0,
-2,0,1,2,1,6,0,-3,0,
2,0,-2,0,1,-6,0,3,0,
2,0,0,0,1,-6,0,3,0,
0,-1,1,0,0,5,0,0,0,
-2,-1,0,2,1,-5,0,3,0,
-2,0,0,0,1,-5,0,3,0,
0,0,2,2,1,-5,0,3,0,
-2,0,2,0,1,4,0,0,0,
-2,1,0,2,1,4,0,0,0,
0,0,1,-2,0,4,0,0,0,
-1,0,1,0,0,-4,0,0,0,
-2,1,0,0,0,-4,0,0,0,
1,0,0,0,0,-4,0,0,0,
0,0,1,2,0,3,0,0,0,
0,0,-2,2,2,-3,0,0,0,
-1,-1,1,0,0,-3,0,0,0,
0,1,1,0,0,-3,0,0,0,
0,-1,1,2,2,-3,0,0,0,
2,-1,-1,2,2,-3,0,0,0,
0,0,3,2,2,-3,0,0,0,
2,-1,0,2,2,-3,0,0,0);

var T = (JDE-2451545.0)/36525.0;
var T2 = T*T;
var T3 = T2*T;
var T4 = T3*T;
var D = rev(297.85036 + 445267.111480*T - 0.0019142*T2 + T3/189474.0);
var M = rev(357.52772 + 35999.050340*T - 0.0001603*T2 - T3/300000.0);
var M1 = rev(134.96298 + 477198.867398*T + 0.0086972*T2 + T3/56250.0);
var F = rev(93.27191 + 483202.017538*T - 0.0036825*T2 + T3/327270.0);
var omega = rev(125.04452 - 1934.136261*T + 0.0020708*T2 + T3/450000.0);
	
var n = c.length/9;
var k = 0;
var DeltaPsi = 0;
var DeltaEpsilon = 0;
for (var i = 0; i < n; i++){
	var arg = rev(c[k]*D+c[k+1]*M+c[k+2]*M1+c[k+3]*F+c[k+4]*omega);
	DeltaPsi += sind(arg)*(c[k+5]+c[k+6]*T);
	DeltaEpsilon += cosd(arg)*(c[k+7]+c[k+8]*T);
	k+=9;
}
DeltaPsi /= 10000.0; //in secondi
DeltaEpsilon /= 10000.0; //in secondi

return Array(DeltaPsi,DeltaEpsilon);
}

function MoonPos(jday,obs) {
// MoonPos calculates the Moon position and distance, based on Meeus chapter 47
// and the illuminated percentage from Meeus equations 48.4 and 48.1
// OPN: This version of MoonPos calculates the position to a precision of about 2' or so
// NB: Time is not taken from obs but from jday (julian day)
// Returns alt, az, hour angle, ra, dec (geocentr!), eclip. long and lat (geocentr!), 
// illumination, distance, brightness and phase angle 
//jday=2448724.5; //per test
	var T=(jday-2451545.0)/36525.0;
	// Moons mean longitude L'
	var LP=rev(218.3164591 + 481267.88134236*T - 0.0013268*T*T + T*T*T/538841.0 - T*T*T*T/65194000.0);
	// Moons mean elongation
	var D=rev(297.8502042 + 445267.1115168*T - 0.00163*T*T + T*T*T/545868.0 - T*T*T*T/113065000.0);
	// Suns mean anomaly
	var M=rev(357.5291092 + 35999.0502909*T - 0.0001536*T*T + T*T*T/24490000.0);
	// Moons mean anomaly M'
	var MP=rev(134.9634114 + 477198.8676313*T + 0.008997*T*T + T*T*T/69699.0 - T*T*T*T/14712000.0);
	// Moons argument of latitude
	var F=rev(93.2720993 + 483202.0175273*T - 0.0034029*T*T - T*T*T/3526000.0 + T*T*T*T/863310000.0);
	// The "further arguments" A1, A2 and A3, and e
	var A1 = rev(119.75 + 131.849*T);
	var A2 = rev(53.09 + 479264.29*T);
	var A3 = rev(313.45 + 481266.484*T);
	var e = 1 - 0.002516*T - 0.000074*T*T;
/*
	// Sum of most significant terms from table 45.A and 45.B (terms less than 0.004 deg / 40 km dropped)
	var Sl = 6288774*sind(MP) + 1274027*sind(2*D-MP) + 658314*sind(2*D) + 213618*sind(2*MP) -
		185116*sind(M)*e - 114332*sind(2*F) + 58793*sind(2*D-2*MP) + 57066*sind(2*D-M-MP)*e +
		53322*sind(2*D+MP) + 45758*sind(2*D-M)*e - 40923*sind(M-MP)*e - 34720*sind(D) - 
		30383*sind(M+MP)*e + 15327*sind(2*D-2*F) - 12528*sind(MP+2*F) + 10980*sind(MP-2*F) +
		10675*sind(4*D-MP) + 10034*sind(3*MP) + 8548*sind(4*D-2*MP) - 7888*sind(2*D+M-MP)*e -
		6766*sind(2*D+M)*e - 5163*sind(D-MP) + 4987*sind(D+M)*e + 4036*sind(2*D-M+MP)*e +
		3994*sind(2*D+2*MP) + 3871*sind(4*D) + 3665*sind(2*D-3*MP) - 2689*sind(M-2*MP)*e -
		2602*sind(2*D-MP+2*F) + 2390*sind(2*D-M-2*MP)*e - 2348*sind(D+MP) + 2236*sind(2*D-2*M)*e*e;
	var Sb = 5128122*sind(F) + 280602*sind(MP+F) + 277602*sind(MP-F) + 173237*sind(2*D-F) +
		55413*sind(2*D-MP+F) + 46271*sind(2*D-MP-F) + 32573*sind(2*D+F) + 17198*sind(2*MP+F) +
		9266*sind(2*D+MP-F) + 8822*sind(2*MP-F) + 8216*sind(2*D-M-F)*e + 4324*sind(2*D-2*MP-F) +
		4200*sind(2*D+MP+F);
	var Sr = (-20905355)*cosd(MP) - 3699111*cosd(2*D-MP) - 2955968*cosd(2*D) - 569925*cosd(2*MP) +
		246158*cosd(2*D-2*MP) - 152138*cosd(2*D-M-MP)*e - 170733*cosd(2*D+MP) - 204586*cosd(2*D-M)*e -
		129620*cosd(M-MP)*e + 108743*cosd(D) + 104755*cosd(M+MP)*e + 79661*cosd(MP-2*F) + 48888*cosd(M)*e;
*/		
	var k = 0;
	var Sl = 0;		
	for (var i = 0; i < n1; i++){
		Sl += c1[k+4]*Math.pow(e,Math.abs(c1[k+1]))*sind(c1[k]*D+c1[k+1]*M+c1[k+2]*MP+c1[k+3]*F);
		k+=5;
	}
	Sl=Math.round(Sl);
			
	k = 0;
	var Sb = 0;		
	for (var i = 0; i < n3; i++){
		Sb += c3[k+4]*Math.pow(e,Math.abs(c3[k+1]))*sind(c3[k]*D+c3[k+1]*M+c3[k+2]*MP+c3[k+3]*F);
		k+=5;
	}
	Sb=Math.round(Sb);
	
	k = 0;
	var Sr = 0;		
	for (var i = 0; i < n2; i++){
		Sr += c2[k+4]*Math.pow(e,Math.abs(c2[k+1]))*cosd(c2[k]*D+c2[k+1]*M+c2[k+2]*MP+c2[k+3]*F);
		k+=5;
	}
	Sr=Math.round(Sr);

  //additive terms
  Sl += 3958*sind(A1) + 1962*sind(LP-F) + 318*sind(A2);
  Sb += -2235*sind(LP) + 382*sind(A3) + 175*sind(A1-F) + 175*sind(A1+F) + 127*sind(LP-MP) - 115*sind(LP+MP);
  // geocentric longitude, latitude
	var mglong = rev(LP+Sl/1000000.0);
	var mglat = Sb/1000000.0;
	// nutation
	var nut = calc_nutazione(jday);
	mglat += nut[0]/3600.0;
	// Obliquity of Ecliptic
	var obl = 23.4392911 + (-46.815*T - 0.00059*T*T + 0.001813*T*T*T)/3600.0;
	obl += nut[1]/3600.0;
	var ra = rev(atan2d(sind(mglong)*cosd(obl)-tand(mglat)*sind(obl),cosd(mglong)));
	var dec = asind(sind(mglat)*cosd(obl)+cosd(mglat)*sind(obl)*sind(mglong));
	var moondat=radec2aa(ra,dec,jday,obs);
	// phase angle (48.4)
	var pa = Math.abs(180.0 - D - 6.289*sind(MP) + 2.100*sind(M) - 1.274*sind(2*D-MP) - 
		0.658*sind(2*D) - 0.214*sind(2*MP) - 0.11*sind(D));
	var k = (1+cosd(pa))/2;
	var mr = Math.round(385000.56+Sr/1000.0);
	var h = moondat[0];
	// correct for parallax, equatorial horizontal parallax, see Meeus p. 337
	h -= asind(6378.14/mr)*cosd(h);
	// brightness, use Paul Schlyter's formula (based on common phase law for Moon)
	var sdat=sunxyz(jday);
	var r = sdat[3];	// Earth (= Moon) distance to Sun in AU
	var R = mr/149598000;	// Moon distance to Earth in AU
	var mag = 0.23 + 5*log10(r*R) + 0.026 * pa + 4.0E-9 * pa*pa*pa*pa
	return new Array(h,moondat[1],moondat[2],ra,dec,mglong, mglat, k, r, mr, mag);
}	// Moonpos()