Moon times - code

Hello,

Does anyone want to share the code that calculates Moon Times? (rise and set)

Thanks,

  • function calculateMoon(moment, lat, lng) 
    	{
    		var rise, set;
    		var nextInfo = Gregorian.info(moment, Time.FORMAT_SHORT);
    		var options = 
    		{
        		:year   => nextInfo.year,
        		:month  => nextInfo.month, 
        		:day    => nextInfo.day,
        		:hour   => 0,
        		:minute => 0,
        		:second => 0
    		};
            		
    		var date = Gregorian.moment(options);
    		var hc = 0.133 * RAD;
    		
    		
    		var h0 = getMoonPosition(date, lat, lng) - hc;
    		var h1, h2, a, b, xe, ye, d, roots, x1, x2, dx;
    			
    		for(var i=1; i<=24; i+=2)
    	    {
    			var duration1 = new Time.Duration(i * Gregorian.SECONDS_PER_HOUR);
    			var date1 = date.add(duration1);
    				
    			var duration2 = new Time.Duration((i+1) * Gregorian.SECONDS_PER_HOUR);
    			var date2 = date.add(duration2);
    				
    			h1 = getMoonPosition(date1, lat, lng) - hc;
    			h2 = getMoonPosition(date2, lat, lng) - hc;
    				
    			a = (h0 + h2) / 2 - h1;
            	b = (h2 - h0) / 2;
            		
            	xe = -b / (2 * a);
            		
            	ye = (a * xe + b) * xe + h1;
            	d = b * b - 4 * a * h1;
            	roots = 0;
    				
    			if (d >= 0) 
    			{
    				dx = Math.sqrt(d) / (a.abs() * 2);
                		
                	x1 = xe - dx;
                	x2 = xe + dx;
                	if (x1.abs() <= 1) 
                	{
                		roots++;
                	}
                	if (x2.abs() <= 1) 
                	{
                		roots++;
                	}
                	if (x1 < -1) 
                	{
                		x1 = x2;
                	}
            	}
    
           	 	if (roots == 1) 
           	 	{
                	if (h0 < 0) 
                	{
                		rise = i + x1;
                	}
               		else 
               		{
               			set = i + x1;
               		}
    
            	} 
            	else if (roots == 2) 
            	{
                	rise = i + (ye < 0 ? x2 : x1);
                	set = i + (ye < 0 ? x1 : x2);
            	}
    
            	if (rise != null && set != null)
            	{
            		break;
            	}
    
            	h0 = h2;
    	    }
            
    		var durationRise = new Time.Duration(rise * Gregorian.SECONDS_PER_HOUR);
    		var dateRise = date.add(durationRise);
    		
    		var nextRise = Gregorian.info(dateRise, Time.FORMAT_SHORT);
    		System.println(nextRise.hour + ":" + nextRise.min); 
    		
    		var durationSet = new Time.Duration(set * Gregorian.SECONDS_PER_HOUR);
    		var dateSet = date.add(durationSet);
    		
    		var nextSet = Gregorian.info(dateSet, Time.FORMAT_SHORT);
    		System.println(nextSet.hour + ":" + nextSet.min); 
    
    		return dateRise;
    	}
    	
    	function getMoonPosition(moment1, lat, lng)
    	{
    		var lw  = RAD * -lng;
            var phi = RAD * lat;
            var dd = moment1.value().toDouble() / DAYS - 0.5 + J1970 - J2000;
            
            var L = RAD * (218.316 + 13.176396 * dd); // ecliptic longitude
            var M = RAD * (134.963 + 13.064993 * dd); // mean anomaly
            var F = RAD * (93.272 + 13.229350 * dd);  // mean distance
    
            var l  = L + RAD * 6.289 * Math.sin(M); // longitude
            var b  = RAD * 5.128 * Math.sin(F);     // latitude
            
            var e = RAD * 23.4397; 
            
            var ra = Math.atan2(Math.sin(l) * Math.cos(e) - Math.tan(b) * Math.sin(e), Math.cos(l)); 
            var dec = Math.asin(Math.sin(b) * Math.cos(e) + Math.cos(b) * Math.sin(e) * Math.sin(l)); 
            
            var H = RAD * (280.16 + 360.9856235 * dd) - lw - ra;
            var h = Math.asin(Math.sin(phi) * Math.sin(dec) + Math.cos(phi) * Math.cos(dec) * Math.cos(H));
            
    		if (h < 0) 
    		{
    			// the following formula works for positive altitudes only.
            	return h; // if h = -0.08901179 a div/0 would occur.
            }
    		
        	return h + 0.0002967 / Math.tan(h + 0.00312536 / (h + 0.08901179)); // altitude correction for refraction
    	}

    I used that code but it doesn't display correct times. Here is the code translate in MonkeyC

  • I've just tested you code and I am good at around 2min.

  • Didn't you chnage anything in the code?

    For me isn't working. For rise I have a difference of about 30 min and for set aroung 4 hours

  • no, I haven't change anything, except the way you deal with date, but even with your way it is ok,

    did you use double for the coordinates?

  • Thank you for your time.

    I have found where the error was. I was seding the LAN and LNG in radians which wasn't correct.