/*
  MDJ DK 2011
  Morten Dam Jørgensen

  Orbital mechanics for top part

	Use as you please but cite as:
		"Morten Dam Joergensen, 2011, http://mdj.dk"


*/

dk = {};
dk.mdj = {};

skyTracerUpdate = 40; // tracer point draw rate
skyUpdateTime = 10000; // how often shall we redraw the sky
objects = {}; //objects

skyTracerCounter = -1;

dk.mdj.orbital = function(){

  //Local functions

   radeg = 180/Math.PI;
   degrad = Math.PI / 180.0;
   
   skyRange = {
     "x" : {"start" : 90, "stop" : 270}, // 0 to 360 full circle
     "y" : {"start": 0, "stop" : 60} // 0 to 90 (horizon to zenith)
   };
   // 
   // start = 90; //deg (East)
   // stop = 270; //deg (West)
   // 

   function skyTransform (a,z) { //simple transform
  
	var canvas_height = $("#sky").height();
	var canvas_width = $("#sky").width();
	
     var xfactor = canvas_width /Math.abs(skyRange.x.stop - skyRange.x.start);
     var yfactor = canvas_height / Math.abs(skyRange.y.stop - skyRange.y.start);
     
     return {"x" : (a - skyRange.x.start)* xfactor, "y" :  (canvas_height - ((z - skyRange.y.start) * yfactor))};
   }

  function sind (ang) {
    return Math.sin( ang * degrad );
  }
  function cosd (ang) {
    return Math.cos( ang * degrad );
  }
  function tand (ang) {
    return Math.tan(ang * degrad);
  }
  
  function asind (arg) {
    return Math.asin(arg * radeg);
  }
  function acosd (arg) {
    return Math.acos(arg* radeg);
  }
  function atand (arg) {
    return Math.atan(arg* radeg);
  }
  function atan2d (y,x) {
    return (radeg * Math.atan2(y,x));
  }
  function trimDeg (deg) {
    deg = deg - Math.floor(deg/360.0)*360.0; //aligns
    if (deg < 0.0) {deg = deg + 360.0;}
    return deg;
  }
  function toRad (r) {
    return (180.0/Math.PI) * r;
  }
  function toDeg (r) {
    return radeg * r;
  }
  function timeToDec (time) { //convert HH:MM:SS time ot decimal hours
    return time.hh + time.mm/60.0 + time.ss /3600; 
  }
  function localTimeToDec (time) { //convert HH:MM:SS time ot decimal hours
    return (time.hh + time.mm/60.0 + time.ss /3600); 
  }
  function degToHr (deg) {
    return deg/15;
  }
  function degToHrCorrected (deg) {
    var hr = deg/15;
    if (hr < 0) {hr = 24 + hr;}
    return hr;
  }
  function HrToDeg (hr) {
    return hr*15;
  }
  function getDay (time) { //day number
    //{"y": 1990, "m" : 4, "d" : 19, "hh" : 0, "mm" : 0, "ss" : 0}
    var d =  parseInt(367 * parseInt(time.y) - 7 * ( parseInt(time.y) + (parseInt(time.m) + 9)/12 ) / 4 + 275* parseInt(time.m)/9 + parseInt(time.d) - 730530);
  
    return d + (time.hh / 24.0);
  }
  
  //eccentric anomaly
  function findE (obj) {
    var E0 = obj.M + (180/Math.PI) * obj.e * sind(obj.M) * (1.0 + obj.e * cosd(obj.M));
    var E1 = E0 - (E0 - (180/Math.PI) * obj.e * sind(E0) - obj.M) / (1 - obj.e * cosd(E0));

    while ((E0 - E1) > 0.005) {
      E0 = E1;
      E1 = E0 - (E0 - (180/Math.PI) * obj.e * sind(E0) - obj.M) / (1 - obj.e * cosd(E0));
    }
    return E1;
  }
  function ecl (d) { //obliquity of the ecliptic
    return 23.4393 - 3.563E-7 * d;
  }
  
  function orbitalObj (obj, d) { //orbital element
    var orbital = {
      
      "name" : obj.name,
      "type" : obj.type,
      "N" : trimDeg(eval(obj.N)),  // N = longitude of the ascending node
      "i" : trimDeg(eval(obj.i)), //  i = inclination to the ecliptic (plane of the Earth's orbit)
      "w" : trimDeg(eval(obj.w)), //  w = argument of perihelion
      "a" : eval(obj.a), //  a = semi-major axis, or mean distance from Sun
      "e" : eval(obj.e),  //  e = eccentricity (0=circle, 0-1=ellipse, 1=parabola)
      "M" : trimDeg(eval(obj.M))   //  M = mean anomaly (0 at perihelion; increases uniformly with time)
    };
    
    orbital.w1 = orbital.N + orbital.w; //longitude of perihelion
    orbital.L = trimDeg(orbital.M + orbital.w1); //mean longitude
    orbital.q = orbital.a*(1.0-orbital.e); //perihelion distance
    orbital.Q = orbital.a*(1.0+orbital.e); //aphelion distance
    orbital.P = Math.pow(orbital.a,1.5); //orbital period if a is in AU
    orbital.oblecl = ecl(d); //Obiquity of ecliptic
    orbital.E = findE(orbital); // Eccentric anomaly
    //console.log(orbital);
    return orbital;
  }

  
  function rv (obj,d) { //distance and true anomaly (seems to work with all bodies)
    //console.log("rv()");
    //eccentric anomaly E from the mean anomaly M and from the eccentricity e (E and M in degrees):
    //E = findE(obj) //obj.M + obj.e*(180/Math.PI) * sind(obj.M) * ( 1.0 + obj.e * cosd(obj.M) );

    //Then compute the bodys distance r and its true anomaly v from:
    var xv = obj.a *  (cosd(obj.E) - obj.e); // r * Math.cos(obj.v) = 
    var yv = obj.a *  (Math.sqrt(1.0 - obj.e*obj.e) * sind((obj.E))); // r * Math.sin(obj.v) 

    var v = trimDeg(Math.atan2( yv, xv ) * radeg);
    var r = Math.sqrt( xv*xv + yv*yv );
    var lon = trimDeg(obj.w + v);
    obj.RV = {"R" : r, "V" : v, "lon" : lon};
    //console.log(obj);
    return {"v" : v, "r" : r};
}

  function sunEclipticRectCoord (obj) {
    //console.log("sunEclipticRectCoord()");
    var x = obj.RV.R * cosd(obj.RV.lon);
    var y = obj.RV.R * sind(obj.RV.lon);
    var z = 0.0;

    
    var xequat = x;
    var yequat = y * cosd(obj.oblecl) + z * sind(obj.oblecl);
    var zequat = y * sind(obj.oblecl) + z * cosd(obj.oblecl);
    obj.ecliptic = {"x" : xequat, "y" : yequat, "z" : zequat};

    var r = obj.RV.R; //AU
    var RA = degToHr(atan2d(yequat, xequat)); //HOUR
    var Decl = atan2d(zequat , Math.sqrt(xequat*xequat + yequat*yequat)); // DEGREE
    obj.RADecl = {"r" : r, "RA" : RA, "Decl" : Decl};
    obj.geoRA = obj.RADecl;
  }
  
  function sunSiderealTime (obj, time) {
    var GMST0 = degToHr(trimDeg( obj.L + 180 ));
    var UT = localTimeToDec(time);
   //console.log(GMST0 + "+" + UT + " + "  + degToHr(time.lon));
    var siderealTime = GMST0 + UT + degToHr(time.lon); //HOUR
    obj.sidtime = siderealTime; //HOUR
    return siderealTime;
  }

  function sunAtAz (obj, time) {
    var HA  = HrToDeg(obj.sidtime - obj.RADecl.RA); //HOUR --> DEG
    
    var x = cosd(HA) * cosd(obj.RADecl.Decl);
    var y = sind(HA) * cosd(obj.RADecl.Decl);
    var z = sind(obj.RADecl.Decl);

    var xhor = x * sind(time.lat) - z * cosd(time.lat);
    var yhor = y;
    var zhor = x * cosd(time.lat) + z * sind(time.lat);
    
    var azimuth = atan2d( yhor, xhor) + 180;
    var altitude = atan2d(zhor, Math.sqrt(xhor*xhor + yhor*yhor));
    obj.ataz = {"azimuth" : azimuth, "altitude" : altitude};

  }
  
  function planetRV (obj, time) {
    var x = obj.a * (cosd(obj.E) - obj.e);
    var y = obj.a * Math.sqrt(1.0 - obj.e*obj.e) * sind(obj.E);
    
    //console.log(x,y);
    
    var r = Math.sqrt(x*x + y*y);
    var v = trimDeg(atan2d(y,x));
    obj.RV = {"R" : r, "V" : v};
    //console.log(r, v);
  }
  function planetEcliptic (obj) {
    var xeclip = obj.RV.R * ( cosd(obj.N) * cosd(obj.RV.V + obj.w) - sind(obj.N) * sind(obj.RV.V + obj.w) * cosd(obj.i));
    var yeclip = obj.RV.R * ( sind(obj.N) * cosd(obj.RV.V + obj.w) + cosd(obj.N) * sind(obj.RV.V + obj.w) * cosd(obj.i));
    var zeclip = obj.RV.R * sind(obj.RV.V + obj.w) * sind(obj.i);
    obj.ecliptic = {"x" : xeclip, "y" : yeclip, "z" : zeclip};

    var r = obj.RV.R; //ER
    var lon = trimDeg(atan2d(yeclip, xeclip)); //Degree
    var lat = atan2d(zeclip , Math.sqrt(xeclip*xeclip + yeclip*yeclip)); // DEGREE
    //console.log(r,lon,lat);
    obj.eclip = {"r": r, "lon" : lon, "lat" : lat}; //heliocentric
  }
  
  function moonPertub (moon, sun) {
    var Ms = sun.M;
    var Mm = moon.M;
    var Ls = sun.L;
    var Lm = moon.L;
    var D = moon.L - sun.L;
    var F = moon.L - moon.N;
    
    var pLon = -1.274 * sind(Mm - 2*D); // (Evection)
    pLon += 0.658 * sind(2*D); // (Variation)
    pLon += -0.186 * sind(Ms); //(Yearly equation)
    pLon += -0.059 * sind(2*Mm - 2*D);
    pLon += -0.057 * sind(Mm - 2*D + Ms);
    pLon += 0.053 * sind(Mm + 2*D);
    pLon += 0.046 * sind(2*D - Ms);
    pLon += 0.041 * sind(Mm - Ms);
    pLon += -0.035 * sind(D);
    pLon += -0.031 * sind(Mm + Ms);
    pLon += -0.015 * sind(2*F - 2*D);
    pLon += 0.011 * sind(Mm - 4*D);
    //console.log(pLon);
    
    var pLat = -0.173 * sind(F - 2*D);
    pLat += -0.055 * sind(Mm - F - 2*D);
    pLat += -0.046 * sind(Mm + F - 2*D);
    pLat += 0.033 * sind(F + 2*D);
    pLat += 0.017 * sind(2*Mm + F);

    var pRadius = -0.58 * cosd(Mm - 2*D);
    pRadius += - 0.46 * cosd(2*D);
    
    var r = moon.eclip.r + pRadius;
    var lon = moon.eclip.lon + pLon;
    var lat = moon.eclip.lat + pLat;

    moon.eclip.lon = lon;
    moon.eclip.lat = lat;
    moon.eclip.r = r;    
  }
  
  function planetRADec (obj) {
    // ecliptic to rectangular
    var x = obj.eclip.r * cosd(obj.eclip.lon) * cosd(obj.eclip.lat);
    var y = obj.eclip.r * sind(obj.eclip.lon) * cosd(obj.eclip.lat);
    var z = obj.eclip.r * sind(obj.eclip.lat);
    
    //rotatre to angle coorsponding to obliquity of ecliptic
    var xequat = x;
    var yequat = y * cosd(obj.oblecl) - z * sind(obj.oblecl);
    var zequat = y * sind(obj.oblecl) + z * cosd(obj.oblecl);
    obj.ecliptic = {"x" : xequat, "y" : yequat, "z" : zequat};    
    
    //convert to spherical
    
    var r = Math.sqrt(xequat*xequat + yequat*yequat + zequat*zequat);
    var RA = trimDeg(atan2d(yequat, xequat)); //HOUR
    var Decl = atan2d(zequat , Math.sqrt(xequat*xequat + yequat*yequat)); // DEGREE
    obj.RADecl = {"r" : r, "RA" : RA, "Decl" : Decl};
    //console.log(obj.RADecl);

  }
  function moonTopocentric(moon, time, sun) {
    var mpar = asind(1.0 / moon.RADecl.r); //do ppar = (8.794/3600) / r for planets
    var gclat = time.lat - 0.1924 * sind(2 * time.lat);
    var rho = 0.99833 + 0.00167 * cosd(2*time.lat);
    
    var HA = trimDeg(HrToDeg(sun.sidtime) - moon.RADecl.RA); //Hour angle - sidtime in hr
    //console.log(mpar, gclat, HA);
    var g = toDeg(atand(tand(gclat) / cosd(HA)));
    
    var topRA = moon.RADecl.RA - mpar * rho * cosd(gclat) * sind(HA) / cosd(moon.RADecl.Decl);
    var topDecl = moon.RADecl.Decl - mpar * rho * sind(gclat) * sind(g - moon.RADecl.Decl) / sind(g);
    //console.log(topRA, topDecl);
    moon.RADecl.RA = topRA;
    moon.RADecl.Decl = topDecl;
    moon.geoRA = {"RA" :  topRA, "Decl" : topDecl, "r" : moon.RADecl.r};
    
  }
  function geoCentric (planet, sun) {
    var x = planet.ecliptic.x + sun.ecliptic.x;
    var y = planet.ecliptic.y + sun.ecliptic.y;
    var z = planet.ecliptic.z + sun.ecliptic.z;
    planet.geoCentric = {"x" : x, "y" : y, "z": z};
    
    var r = Math.sqrt(x*x + y*y + z*z);
    var RA = (atan2d(y, x));
    var Decl = atan2d(z , Math.sqrt(x*x + y*y));
    
    planet.geoRA = {"RA" : RA, "Decl" : Decl, "r" : r};
  }
  function planetAtAz (obj, time, sun) {
    var HA  = HrToDeg(sun.sidtime - degToHr(obj.geoRA.RA)); //HOUR --> DEG
    //console.log(HA, sun.sidtime);
    var x = cosd(HA) * cosd(obj.geoRA.Decl);
    var y = sind(HA) * cosd(obj.geoRA.Decl);
    var z = sind(obj.geoRA.Decl);

    var xhor = x * sind(time.lat) - z * cosd(time.lat);
    var yhor = y;
    var zhor = x * cosd(time.lat) + z * sind(time.lat);
    
    var azimuth = atan2d( yhor, xhor) + 180 ;
    var altitude = atan2d(zhor, Math.sqrt(xhor*xhor + yhor*yhor));
    obj.ataz = {"azimuth" : azimuth, "altitude" : altitude};

  }

  return {
  //External functions
  
  getElements : function(url) {
    $.getJSON(url, {async:false},
     function(json){

		raphaelCan = Raphael($("#sky")[0]);
        dk.mdj.orbital.saveObjects(json);
		


        for (var i=0; i < json.orbitals.length; i++) {
          $("#sky").append("<div class='planet text planetimg' id='p"+ i +"'>" + json.orbitals[i].name + "</div>");
          
            if (objects.orbitals.name = "Moon") {
                // var now = new Date();
                // var time = {"y": now.getUTCFullYear(), "m" : (now.getUTCMonth()+1), "d" : now.getUTCDate(), "hh" : now.getUTCHours(), "mm" : now.getUTCMinutes(), "ss" : now.getUTCSeconds()};
                // 
                // var phase = dk.mdj.orbital.moonphase(time.y, time.m, time.d);
                var whi = objects.orbitals[i].rep = raphaelCan.circle(10 + i, 10 + i, objects.orbitals[i].size).attr({"fill" : objects.orbitals[i].color, "stroke-width" : 0}).hide();
    			objects.orbitals[i].rep = whi
            } else {
                objects.orbitals[i].rep = raphaelCan.circle(10 + i, 10 + i, objects.orbitals[i].size).attr({"fill" : objects.orbitals[i].color, "stroke-width" : 0}).hide();
            }
          

			objects.orbitals[i].trace = raphaelCan.set();

			$("#p" + i).hover(
		  		function () {
					var j = parseInt($(this).attr("id").substring(1));
					// console.log(planets[j]);
					// console.log(objects.orbitals);
					// console.log(objects[j]);
					// console.log($(this).position());
					var hoverbox = $("<div/>");
					hoverbox.addClass("hoverbox");

					var topbox = $("<div/>");					
					topbox.css("min-height" ,"50px");

					
					if (planets[j].name == "Moon"){
					    var img = $("<div />");
					    var phase = dk.mdj.orbital.moonphase(objects.time.y, objects.time.m, objects.time.d);
					    var row = Math.floor(phase / 4);
					    var column = phase - row *4;
                        // console.log(phase, row, column);
					    img.css({"width" : 30, "height" : 31, "background-image" : "url(http://mdj.dk/ndxz-studio/site/mdjb/phases.png)", "background-position" :  (-column * 24 +4)+ "px " + (-row * 31 - 5 ) + "px"})
					    img.addClass("planet_img");				        
    					
				    } else {
    					var img = $("<img/>");
    					img.attr("src", objects.orbitals[j].img);
    					img.addClass("planet_img");				        
				    }
					
					topbox.append(img)
					// hoverbox.append(img);
					topbox.append($("<b>" + planets[j].name + "</b> <br/><i>" + planets[j].type + "</i>"));
					hoverbox.append(topbox);
					// Add cool info here (rise/set time, position, magnitude, etc...)
					hoverbox.append("<br/>RA/DE: " + degToHrCorrected(planets[j].geoRA.RA).toPrecision(4) + "h/" + planets[j].geoRA.Decl.toPrecision(4) + "&deg;");
					hoverbox.append("<br/>Az/Alt: " + planets[j].ataz.azimuth.toPrecision(4) + "&deg;/" + planets[j].ataz.altitude.toPrecision(4) + "&deg;");
					
					if (planets[j].name == "Moon"){
						hoverbox.append("<br/>Distance: " + (planets[j].geoRA.r * 6367.5).toPrecision(6) + " km");
                        // hoverbox.append("<br/>Phase: "+ dk.mdj.orbital.moonphase(objects.time.y, objects.time.m, objects.time.d));
						
					} else {
						hoverbox.append("<br/>Distance: " + planets[j].geoRA.r.toPrecision(4) + " AU");
						
					}
					
					// hoverbox.append("<br/> " + planets[j].geoRA.rv + "");

					var l = $(this).offset().left + 10;
					var t = $(this).offset().top  + 10;
					hoverbox.css({"left" : l + "px", "top" : t + "px"});

		    		$("body").append(hoverbox);

		  }, 
		  	function () {
		    	$("body").find("div.hoverbox").remove();

		  }
		).click(function() { // Link to documentation
			var j = parseInt($(this).attr("id").substring(1));
			document.location = objects.orbitals[j].link;
		});
        }

//Uncomment to add local positioning
		// $.getJSON('/ajax/getPosition',{},
		//  function(json){
		//     //10.8, "lat" : 56
		
		
		// USe HTML 5 navigator.geolocation.getCurrentPosition(callbackfun) instead
		
			pos = {
				"lat" : 56,
				"lon" : 10.8
			};
			// if (json.lat != "from") {
			// pos.lon = parseInt(json.lon);
			// pos.lat = parseInt(json.lat);
			// }

	        dk.mdj.orbital.drawSky();
	        sky = setInterval('dk.mdj.orbital.drawSky();',skyUpdateTime);
	       // visits = setInterval('dk.mdj.orbital.renderVisits()', 10000)			
		//});
		

    });

    
  },
  
  saveObjects : function(json) {
    objects = json;
    return "saved";
  },
  calcPositions : function() {
    planets = []; //global objects
    var now = new Date();
//    var time = {"y": now.getUTCFullYear(), "m" : (now.getUTCMonth()+1), "d" : now.getUTCDate(), "hh" : now.getUTCHours(), "mm" : now.getUTCMinutes(), "ss" : now.getUTCSeconds(), "lon" : 10.8, "lat" : 56};
//cern 
   var time = {"y": now.getUTCFullYear(), "m" : (now.getUTCMonth()+1), "d" : now.getUTCDate(), "hh" : now.getUTCHours(), "mm" : now.getUTCMinutes(), "ss" : now.getUTCSeconds(), "lon" : pos.lon, "lat" : pos.lat};
    objects.time = time;
    objects.updated = time;
    var d = getDay(time);
    objects.time.day = d;
    
    var sun = orbitalObj(objects.orbitals[0], d);
    var RV = rv(sun,d);
    sunEclipticRectCoord(sun);
    sunSiderealTime(sun, time);
    sunAtAz(sun, time);
    planets[0] = sun;
    //Moon
    var moon = orbitalObj(objects.orbitals[1], d);
    planetRV(moon, time);
    planetEcliptic(moon);
    moonPertub (moon, sun); 
    planetRADec( moon);
    moonTopocentric(moon, time, sun);
    planets[1] = moon;
    planetAtAz(moon, time, sun);

    for (var i = objects.orbitals.length - 1; i >= 2; i--){
      var planet = orbitalObj(objects.orbitals[i], d); // orbital properties
      planetRV(planet, time); // planet distance and true anomaly
      planetEcliptic(planet); // planet longitude and latitude 
      planetRADec(planet);
      geoCentric(planet, sun); //RA and Decl in geocentric values
      planets[i] = planet;

      planetAtAz(planet, time, sun);
      
    }
    
  },
  drawSky : function() {
    dk.mdj.orbital.calcPositions();  //Calculate current positions
    dk.mdj.orbital.renderSky(); //Draw the objects
    
  },
  
  renderVisits: function() {
    d = new Date();
    var date =  d.getUTCFullYear() + "-" + (d.getUTCMonth()+1) + "-" + d.getUTCDate() + " " + d.getUTCHours() + ":" + d.getUTCMinutes() +  ":" + d.getUTCSeconds();

    $.getJSON('/ajax/orbital/getUsersSince',{
     t: date, y : d.getUTCFullYear(), m: (d.getUTCMonth()+1), d: d.getUTCDate(), h: d.getUTCHours(), mm : d.getUTCMinutes(),s : d.getUTCSeconds()},
     function(json){


        //$("#data").empty().append("Visits since " + date + ": " + json.c);
        
    });
    
    
  },
  
  renderSky : function() {


 
 var dim = {"w" : $("#skyC").width(), "h" : $("#skyC").height()};

  skyTracerCounter += 1;
//$("#data").empty().append("Julian day: " + objects.time.day + " time: " + objects.time.y + "/" + objects.time.m + "/" + objects.time.hh + "/" + objects.time.mm+ "/" + objects.time.ss + "  At longitude: " +  objects.time.lon + " Latitude: " + objects.time.lat + "<table><tr><td>Name</td><td>RA</td><td>Decl</td><td>Distance</td></tr></table>");
for (var i = planets.length - 1; i >= 0; i--){
//$("#data table").append("<tr><td>" + planets[i].name + "</td><td>" + degToHrCorrected(planets[i].geoRA.RA) + "</td><td>" + planets[i].geoRA.Decl + "</td><td>" + planets[i].geoRA.r + "</td></tr>");

 if ((planets[i].ataz.azimuth > skyRange.x.start && planets[i].ataz.azimuth < skyRange.x.stop) && (planets[i].ataz.altitude > skyRange.y.start && planets[i].ataz.altitude <  skyRange.y.stop)) {
      $("#p"+ i).show("medium"); //show text
   		objects.orbitals[i].rep.show();


      var cord = skyTransform(planets[i].ataz.azimuth,planets[i].ataz.altitude);
		
       $("#p"+ i).css("left",(cord.x-2) + "px").css("top",(cord.y-2) + "px");
		objects.orbitals[i].rep.animate({"cx" : cord.x, "cy" : cord.y});
	

       if (skyTracerCounter % skyTracerUpdate == 0) {
		// Add tracer
		objects.orbitals[i].trace.push(raphaelCan.circle(cord.x,cord.y, 1.2).attr({"fill" : "#333", "stroke-width" : 0}).toBack());		
       }
     }
  else {
      $("#p"+ i).hide(); //hide text    
    	objects.orbitals[i].rep.hide();

  }
  // ctx.restore();

}

  },
  flare : function() {
    // add a flare if new visitors enter site

    
  },
  
  julday: function(year, month, day) {
  	if (year < 0) { year ++; }
  	var jy = parseInt(year);
  	var jm = parseInt(month) +1;
  	if (month <= 2) {jy--;	jm += 12;	} 
  	var jul = Math.floor(365.25 *jy) + Math.floor(30.6001 * jm) + parseInt(day) + 1720995;
  	if (day+31*(month+12*year) >= (15+31*(10+12*1582))) {
  		ja = Math.floor(0.01 * jy);
  		jul = jul + 2 - ja + Math.floor(0.25 * ja);
  	}
  	return jul;
  },
  
  moonphase : function(year,month,day) {
      //  (0 to 29, where 0=new moon, 15=full etc.) 
  	n = Math.floor(12.37 * (year -1900 + ((1.0 * month - 0.5)/12.0)));
  	RAD = 3.14159265/180.0;
  	t = n / 1236.85;
  	t2 = t * t;
  	as = 359.2242 + 29.105356 * n;
  	am = 306.0253 + 385.816918 * n + 0.010730 * t2;
  	xtra = 0.75933 + 1.53058868 * n + ((1.178e-4) - (1.55e-7) * t) * t2;
  	xtra += (0.1734 - 3.93e-4 * t) * Math.sin(RAD * as) - 0.4068 * Math.sin(RAD * am);
  	i = (xtra > 0.0 ? Math.floor(xtra) :  Math.ceil(xtra - 1.0));
  	j1 = dk.mdj.orbital.julday(year,month,day);
  	jd = (2415020 + 28 * n) + i;
  	return (j1-jd + 30)%30;
  }
    
};
}();