Calculate and Draw Moon Phase

it is very enlightening for my personal project...

diligent hands rule....

In the MoonAge function ip is calculated using:

int j = JulianDate(d, m, y);
//Calculate the approximate phase of the moon
ip = (j + 4.867) / 29.53059;
ip = ip - Math.Floor(ip);

I have searched the internet and have been unable to find reference to j + 4.867

Internet consensus appears to be:
In vb.net..

Dim DaySince As Double = j - 2451550.1
Dim NewMoons As Double = DaySince / 29.53059
ip = NewMoons - Math.Floor(NewMoons)
ag = ip * 29.53059

What is the theory behind using j + 4.867 ?
Is it a restructured formula to allow more direct use in sub DrawMoon?
If yes, how was it derived?

Thanks for your work,
David

Hi Mostafa,

Thanks for a great program. I needed some code to calculate the phase of the moon and your program seemed to be just what I wanted. However, when I downloaded the VB code and ran it I found a problem. The graphic displaying the moon didn’t look quite right. I realised what the problem was. I live in New Zealand and our view of things such as the moon is somewhat different to you people in the northern hemisphere (we are upside down remember). It looked perfect if I stood on my head (not very practical). Could you please tell me what changes to your code I would need to make so that the moon graphic displays correctly for us people down here in the antipodes.

Many thanks, Terry.

I think you need to adjust time zone difference because you are in Australia...

diligent hands rule....

Hi
Calculation JulianDate could be done a lot easier:

C#

private double DateTimeToJulianDate(DateTime date)
{
    return date.ToOADate() + 2415018.5;
}

private double DateTimeToJulianDate(int d, int m, int y)
{
    return DateTimeToJulianDate(new DateTime(y, m, d));
}

And the long explanation:

C#

/*if (Month< 3) ...: To make the magic numbers work our right, they're putting February at the 'end' of the year.
      (153 * Month - 457) / 5: Wow, that's some serious magic numbers.

      Normally, the number of days in each month is 31 28 31 30 31 30 31 31 30 31 30 31, but after that adjustment in the if statement,
      it becomes 31 30 31 30 31 31 30 31 30 31 31 28. Or, subtract 30 and you end up with 1 0 1 0 1 1 0 1 0 1 1 -2.
      They're creating that pattern of 1s and 0s by doing that division in integer space.

      Re-written to floating point, it would be(int)(30.6 * Month - 91.4).
      30.6 is the average number of days per month, excluding February(30.63 repeating, to be exact).
      91.4 is almost the number of days in 3 average non-February months. (30.6 * 3 is 91.8).

      So, let's remove the 30, and just focus on that 0.6 days. If we multiply it by the number of months, and then truncate to an integer, we'll get a pattern of 0s and 1s.
      •  0.6 * 0 = 0.0 -> 0
      •  0.6 * 1 = 0.6 -> 0 (difference of 0)
      •  0.6 * 2 = 1.2 -> 1 (difference of 1)
      •  0.6 * 3 = 1.8 -> 1 (difference of 0)
      •  0.6 * 4 = 2.4 -> 2 (difference of 1)
      •  0.6 * 5 = 3.0 -> 3 (difference of 1)
      •  0.6 * 6 = 3.6 -> 3 (difference of 0)
      •  0.6 * 7 = 4.2 -> 4 (difference of 1)
      •  0.6 * 8 = 4.8 -> 4 (difference of 0)

      See that pattern of differences in the right? That's the same pattern in the list above, the number of days in each month minus 30.
      The subtraction of 91.8 would compensate for the number of days in the first three months, that were moved to the 'end' of the year,
      and adjusting it by 0.4 moves the successive differences of 1 (0.6 * 4 and 0.6 * 5 in the above table) to align with the adjacent months that are 31 days.

      Since February is now at the 'end' of the year, we don't need to deal with its length.
      It could be 45 days long (46 on a leap year), and the only thing that would have to change is the constant for the number of days in a year, 365.
      Note that this relies on the pattern of 30 and 31 month days.If we had two months in a row that were 30 days, this would not be possible.
      365 * Year: Days per year (Year / 4) - (Year / 100) + (Year / 400): Plus one leap day every 4 years, minus one every 100, plus one every 400.

      + 1721119: This is the Julian Date of March 2nd, 1 BC.
      Since we moved the 'start' of the calendar from January to March, we use this as our offset, rather than January 1st.Since there is no year zero,
      1 BC gets the integer value 0. As for why March 2nd instead of March 1st, I'm guessing that's because that whole month calculation was still a little off at the end.
      If the original writer had used - 462 instead of - 457 (- 92.4 instead of - 91.4 in floating point math), then the offset would have been to March 1st.*/

Could you add the hilal counting please? Because the islamic can be helped much with your program.

I am sorry, I haven't equations for this option.

Kindly add more controls and event handlers for moonrise, moonset, next crescent, etc...

I haven't equations for moonrise, moonset and other, but I shall search for this.
Thank you very much.

Some good C# code for sunrise/moonrise can be found @ http://wiki.crowe.co.nz/Calculate Sunrise%2fSunset.ashx.

Nice job. Keep it up...

Excellent!

Thank you my friend.

Nice

First of all, I must say that I don't like the term "Julian date" as it's ambiguous; it could mean a date in the Julian calendar. a date in the Julian period (which is how your code uses it), or the number of days since the start of the year (which is how the US Naval Observatory website, among others, uses the term). I prefer to say "Julian-calendar date", "Julian-period number" and "day of year" respectively.

More seriously, your routine to calculate the Julian-period number has a major error; it assumes that only dates up to 2299160 (1582-10-04 Julian) are in the Julian calendar and all others are Gregorian. This is only true of some countries; Britain and its (then) territories was Julian up to 1752-09-02 Julian (1752-09-13 Gregorian), and some countries were still using the Julian calendar up to the early 20th century. See the Calendar FAQ section 2.2.4 for further information. (The Calendar FAQ itself has an error here; it refers to the Western calendar as the "Christian" calendar, despite the fact that it was introduced in 45BC, over 70 years before the establishment of Christianity. I prefer to call it the "Roman" or "solar" calendar.)

My main reason for this post is that I think I have a simpler (and probably more efficient) way of calculating the age of the moon for dates 1900-2099 inclusive; calculate the number of days since 1899-12-31, subtract 1, and take the result modulo 29.530588853 (the latter figure of course being the length in days of a lunar month). The following is my program for this in TI-89 programming language (© lines are comments):

© Age of the moon
moonage(y,m,d)
Func
Return mod(schn(y,m,d)-1,29.530588853)
EndFunc

The above routine depends on the following:
mod(x,y): Remainder of x/y. Equivalent to x-y*floor(x/y); on the TI-89, can take real arguments.

© Days since 1899-12-31 (Schecter's Number)
© Valid for range 1900-2099 (1901-2099 if only truncation downward, not toward zero, is available)
schn(y,m,d)
Func
Return iPart((1461*(y-1900)-1)/4)+doy(y,m,d)
EndFunc

iPart(x): Integer part of x, truncated toward zero

© Days from start of year to given date (inclusive)
doy(y,m,d)
Func
Return iPart(275*m/9)-iPart((m+9)/12)*(notleap(y)+1)+d-30
EndFunc

© 0 if leap year, 1 otherwise (assumes Gregorian calendar)
© For Julian calendar, omit the second and third parts
notleap(y)
Func
Return iPart((mod(y,4)+2)/3)-iPart((mod(y,100)+99)/100)+iPart((mod(y,400)+399)/400)
EndFunc

Oops, in the above I said "calculate the number of days since 1899-12-31 [and] subtract 1" without explaining why. The reason being, of course, that 1900-1-1 (day number 1 in the applicable range) was a new moon.

Dear Robert,
Thank you very much for your important and good knowledge.
This article and its code from my older papers!

Hello,

I am programming in VB6 and i want to translate DrawMoon procedure to VB6.
any idea ?

Regards.

1- Let your Form has following controls:
MonthView, its Name = MyCalendar
PictureBox, its Name = PicMoon, change its BackColor to Dark Blue, and its ForeColor to White
Label, its Name = lblAge
CommandButton, its Name = btnClose, its Caption = "Close"
CommandButton, its Name = btnToDay, its Caption = "Moon today"

2- In Declarations section, write following code, or copy and paste it:
Option Explicit
Dim ip As Double
Dim ag As Double

Private Sub DrawMoon()
Dim x As Single, y As Single
Dim X1 As Single, X2 As Single
Dim R As Single
Dim Pixels As Single
Dim Phase As Double
Phase = ip
Pixels = ScaleHeight
PicMoon.Cls
PicMoon.Scale (-1.5, -1.5)-(1.5, 1.5)
For y = 0 To 1 Step 4 / Pixels
x = Sqr(1 - y * y)
PicMoon.Line (-x, y)-(x, y), 0
PicMoon.Line (-x, -y)-(x, -y), 0
R = 2 * x
If Phase < 0.5 Then
X1 = -x
X2 = R - 2 * Phase * R - x
Else
X1 = x
X2 = x - 2 * Phase * R + R
End If
PicMoon.Line (X1, y)-(X2, y)
PicMoon.Line (X1, -y)-(X2, -y)
Next y
PrintAge
End Sub

Private Function JulianDate(d As Integer, m As Integer, y As Integer) As Long
Dim mm As Integer, yy As Integer
Dim k1 As Long, k2 As Long, k3 As Long
Dim j As Long
yy = y - Int((12 - m) / 10)
mm = m + 9
If (mm >= 12) Then
mm = mm - 12
End If
k1 = Int(365.25 * (yy + 4712))
k2 = Int(30.6001 * mm + 0.5)
k3 = Int(Int((yy / 100) + 49) * 0.75) - 38
j = k1 + k2 + d + 59
If (j > 2299160) Then
j = j - k3
End If
JulianDate = j
End Function

Private Sub MoonAge(ByVal d As Integer, ByVal m As Integer, ByVal y As Integer)
Dim j As Long
j = JulianDate(d, m, y)
ip = (j + 4.867) / 29.53059
ip = ip - Int(ip)
If (ip < 0.5) Then
ag = ip * 29.53059 + 29.53059 / 2
Else
ag = ip * 29.53059 - 29.53059 / 2
End If
ag = Int(ag) + 1
End Sub

Private Sub PrintAge()
Dim theAge As String
theAge = "Moon age"
theAge = theAge + " " + ":" + " " + Str(ag)
If (ag = 1) Then
theAge = theAge + " " + "day"
Else
theAge = theAge + " " + "days"
End If
lblAge.Caption = theAge
End Sub

Private Sub ShowMoon()
YourChoice
DrawMoon
PrintAge
End Sub

Private Sub YourChoice()
Dim Aday, Amonth, Ayear As Integer
Aday = MyCalendar.Day
Amonth = MyCalendar.Month
Ayear = MyCalendar.Year
MoonAge Aday, Amonth, Ayear
End Sub

3- Write following events:
Private Sub btnClose_Click()
Unload Me
End Sub

Private Sub btnToDay_Click()
MyCalendar.Value = Now
YourChoice
Form_Paint
End Sub

Private Sub Form_Load()
PicMoon.ScaleMode = 3
MyCalendar.Value = Now
YourChoice
End Sub

Private Sub Form_Paint()
DrawMoon
End Sub

Private Sub MyCalendar_DateClick(ByVal DateClicked As Date)
YourChoice
Form_Paint
End Sub
--------------------------------------------------------------------------------------------

Please let me know about the result.

Thanks,
Mostafa Kaisoun

It is amazing, but the quality isn't good at all

Very nice article!

A fine job - thank you for sharing.

Thank you very much.

I like how you draw the moon.
Of course the lunar phase can be figured easily using the algorithm for the hebrew calendar, and I assume the muslam calendar as well? On hebrew calendar, the first day of the month is always new moon and 14th is full moon. I don't know if muslam calendar has the month starting on new moon or full moon, but it's lunar-based.

Of course, seeing as you're from Egypt, you probably know that Smile | :)

Thank you very much,
Muslam calendar has the month starting on new moon.
Thank you again.

Calculate and Draw Moon Phase

Introduction

Background

About the Code

Final Words

License

Comments and Discussions