Want to avoid the missteps to gaining all the benefits of the cloud? Learn more about the different assessment options from our Cloud Advisory team.
'Position of the Sun at 11:00 UT on 1997 August 7th
'1. Find the days before J2000.0 (d) from the table
Dim test_d As Double
'11 = hour
'24 = day
test_d = 11 / 24 + 212 + 7 - 1096.5 '= -877.04167
Text1.Text = test_d
Dim y, m, d, hours, minutes, seconde As Integer
Dim h As Double
y = 1997 'Year(Now)
m = 8 'Month(Now)
d = 7 'Day(Now)
hours = 11 'Hour(Now)
minutes = 0 'Minute(Now)
seconds = 0 'Second(Now)
h = CDbl(hours + (minutes / 60#) + (seconds / 3600#))
dd = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 + 275 * m \ 9 + d - 730531.5 + h / 24
Text2.Text = dd
'2. Find the Mean Longitude (L) of the Sun
L = 280.461 + 0.9856474 * dd
= -583.99284 + 720
(add multiples of 360 to bring in range 0 to 360)
= 136.00716
Dim y, x, a, alpha, delta As Double
y = Cos(epsilon) * Sin(lambda)
x = Cos(lambda)
a = ArcTan(y / x)
If x < 0 Then
alpha = a + 180
End If
If y < 0 And x > 0 Then
alpha = a + 360
Else
alpha = a
End If
'y = 0.6489924
Text7.Text = y
'x = -0.7068507
Text8.Text = x
'a = -42.556485
Text9.Text = a
'alpha = -42.556485 + 180
'alpha = a + 180
'= 137.44352 (degrees)
Text10.Text = alpha
delta = ArcSin(Sin(epsilon) * Sin(lambda))
'= 16.342193 degrees
Text11.Text = delta
Public Function ArcSin(x As Double) As Double
'Inverse Sine
On Error Resume Next
ArcSin = Atn(x / Sqr(-x * x + 1))
On Error GoTo 0
End Function
Public Function ArcTan(x As Double) As Double
'Inverse Tangent
On Error Resume Next
ArcTan = Atn(x) * (180 / pi)
On Error GoTo 0
End Function
'Position of the Sun at 11:00 UT on 1997 August 7th
'1. Find the days before J2000.0 (d) from the table
Dim test_d As Double
'11 = hour
'24 = day
test_d = 11 / 24 + 212 + 7 - 1096.5 '= -877.04167
text1.Text = test_d
Dim y, m, d, hours, minutes, seconds As Integer
Dim h As Double
Dim dd As Double
y = 1997 'Year(Now)
m = 8 'Month(Now)
d = 7 'Day(Now)
hours = 11 'Hour(Now)
minutes = 0 'Minute(Now)
seconds = 0 'Second(Now)
h = CDbl(hours + (minutes / 60.0#) + (seconds / 3600.0#))
dd = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 + 275 * m \ 9 + d - 730531.5 + h / 24
text2.Text = dd
'2. Find the Mean Longitude (L) of the Sun
Dim L As Double
' according to the order of operations, * or / before + or -, this is the order it would calculate
L = ((0.9856474 * dd) + 280.461)
' assuming L is a negative number... if not, I'm not sure if you'd subtract if L > 360
' but, you are trying to get L in the range of 0 to 360
While L < 0
L = L + 360
End While
text3.Text = L
'3. Find the Mean anomaly (g) of the Sun
Dim g As Double
g = 357.528 + 0.9856003 * dd
'= -506.88453 + 720
g = g + 720
'= 213.11547
text4.Text = g
'4. Find the ecliptic longitude (lambda) of the sun
Dim lambda As Double
'lambda = L + 1.915 * Sin(g) + 0.02 * Sin(2 * g)
'lambda = L + 1.915 * RadiansToDegrees(Sin(DegreesToRadians(g))) + 0.02 * RadiansToDegrees(Sin(DegreesToRadians(2 * g)))
lambda = L + 2.060165 * Sin(g) + 0.02 * Sin(2 * g)
'= 134.97925
text5.Text = lambda ' = 135.05047 should be 134.97925 !!!
'5. Find the obliquity of the ecliptic plane (epsilon)
Dim epsilon As Double
epsilon = 23.439 - (0.0000004 * dd)
'= 23.439351
text6.Text = epsilon
'6. Find the Right Ascension (alpha) and Declination (delta) of the(Sun)
Dim XX As Double
Dim YY As Double
Dim AA As Double
Dim alpha As Double
Dim delta As Double
YY = Cos(DegreesToRadians(epsilon)) * Sin(DegreesToRadians(lambda))
XX = Cos(DegreesToRadians(lambda))
'YY = 0.6489924
'XX = -0.7068507
'AA = arctan(YY / XX)
AA = RadiansToDegrees(Atan(YY / XX))
'AA = -42.556485
If XX < 0 Then
alpha = AA + 180
ElseIf (YY < 0 And XX > 0) Then
alpha = AA + 360
Else
alpha = AA
End If
'alpha = -42.556485 + 180
'= 137.44352 '(degrees)
text7.Text = alpha
'delta = arcsin(Sin(epsilon) * Sin(lambda))
'delta = Asin(Sin(epsilon) * Sin(lambda))
delta = RadiansToDegrees(Asin(Sin(DegreesToRadians(epsilon)) * Sin(DegreesToRadians(lambda))))
'= 16.342193 degrees
text8.Text = delta
End Sub
Private Function DegreesToRadians(ByVal X As Double) As Double
Const pi = 3.14159265358979
'To convert degrees to radians, multiply degrees by pi/180.
DegreesToRadians = X * pi / 180
End Function
Private Function RadiansToDegrees(ByVal X As Double) As Double
Const pi = 3.14159265358979
'To convert radians to degrees, multiply radians by 180/pi.
RadiansToDegrees = X * (180 / pi)
End Function
'Position of the Sun at 11:00 UT on 1997 August 7th
'1. Find the days before J2000.0 (d) from the table
Dim test_d As Double
'11 = hour
'24 = day
test_d = 11 / 24 + 212 + 7 - 1096.5 '= -877.04167
Text1.Text = test_d
Dim y, m, d, hours, minutes, seconds As Integer
Dim h As Double
Dim dd As Double
y = 1997 'Year(Now)
m = 8 'Month(Now)
d = 7 'Day(Now)
hours = 11 'Hour(Now)
minutes = 0 'Minute(Now)
seconds = 0 'Second(Now)
' y = Year(Now)
' m = Month(Now)
' d = Day(Now)
' hours = Hour(Now)
' minutes = Minute(Now)
' seconds = Second(Now)
h = CDbl(hours + (minutes / 60#) + (seconds / 3600#))
dd = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 + 275 * m \ 9 + d - 730531.5 + h / 24
Text2.Text = dd
'2. Find the Mean Longitude (L) of the Sun
Dim L As Double
' according to the order of operations, * or / before + or -, this is the order it would calculate
L = ((0.9856474 * dd) + 280.461)
' assuming L is a negative number... if not, I'm not sure if you'd subtract if L > 360
' but, you are trying to get L in the range of 0 to 360
While L < 0
L = L + 360
Wend
Text3.Text = L
'3. Find the Mean anomaly (g) of the Sun
Dim g As Double
g = 357.528 + 0.9856003 * dd
'= -506.88453 + 720
g = g + 720
'= 213.11547
Text4.Text = g
'4. Find the ecliptic longitude (lambda) of the sun
Dim lambda As Double
'lambda = L + 1.915 * Sin(g) + 0.02 * Sin(2 * g)
'lambda = L + 1.915 * RadiansToDegrees(Sin(DegreesToRadians(g))) + 0.02 * RadiansToDegrees(Sin(DegreesToRadians(2 * g)))
lambda = L + 2.060165 * Sin(g) + 0.02 * Sin(2 * g)
'= 134.97925
Text5.Text = lambda ' = 135.05047 should be 134.97925 !!!
'5. Find the obliquity of the ecliptic plane (epsilon)
Dim epsilon As Double
epsilon = 23.439 - (0.0000004 * dd)
'= 23.439351
Text6.Text = epsilon
'6. Find the Right Ascension (alpha) and Declination (delta) of the(Sun)
Dim XX As Double
Dim YY As Double
Dim AA As Double
Dim alpha As Double
Dim delta As Double
YY = Cos(DegreesToRadians(epsilon)) * Sin(DegreesToRadians(lambda))
XX = Cos(DegreesToRadians(lambda))
'YY = 0.6489924
Text7.Text = YY
'XX = -0.7068507
Text8.Text = XX
AA = ArcTan(YY / XX)
'AA = RadiansToDegrees(Atan(YY / XX))
'AA = RadiansToDegrees(ArcTan(YY / XX))
'AA = -42.556485
Text9.Text = AA
If XX < 0 Then
alpha = AA + 180
ElseIf (YY < 0 And XX > 0) Then
alpha = AA + 360
Else
alpha = AA
End If
'alpha = -42.556485 + 180
'= 137.44352 '(degrees)
Text10.Text = alpha
'delta = arcsin(Sin(epsilon) * Sin(lambda))
'delta = Asin(Sin(epsilon) * Sin(lambda))
'delta = RadiansToDegrees(Asin(Sin(DegreesToRadians(epsilon)) * Sin(DegreesToRadians(lambda))))
delta = RadiansToDegrees(ArcSin(Sin(DegreesToRadians(epsilon)) * Sin(DegreesToRadians(lambda))))
'= 16.342193 degrees
Text11.Text = delta
'7. Find the Aziumuth (Z) and altitude (A) of the Sun, we first
'find the local sidereal time, then the hour angle of the Sun.
'Long is the longitude (west negative) of the place, and lat is
'the latitude. LST and RA are taken in degrees. Note that the
'calculation of LST needs about 12 places of precision.
' I tried with two different variable types Decimal and Double
' I ended up with the same result, however, I think it might
' need to be this decimal type, and the other Double values
' typecast to decimals to get greater accuracy, as the value
' rounded slightly
Dim LST_DEC As Decimal
LST_DEC = 280.46061837D + (360.98564736629D * dd) + (-1.9166702079166957)
While LST_DEC < 0
LST_DEC += 360
End While
text12.Text = LST_DEC.ToString("###,##0.000000000000")
Dim LST As Double
LST = CDbl(280.46061837D + (360.98564736629D * dd) + (-1.9166702079166957))
While LST < 0
LST += 360
End While
text13.Text = LST.ToString("###,##0.000000000000")
'LST = -316320.911064 degs
'= 119.088936 degs
Dim ha As Double
ha = LST - alpha
While ha < 0
ha += 360
End While
text14.Text = LST.ToString("###,##0.000000000000")
'ha = 119.088936 - 137.44352degs
'= -18.354584 degs
'= 341.64662711177084 (they forgot to 0-360 the value)
Dim DAY0 As Double
Dim M0 As Double
Dim L0 As Double
Dim EQC1 As Double
Dim EQC2 As Double
Dim OBLQ As Double
Dim TheYear As Double
Dim TheDay As Double
Dim TheHour As Double
Dim TheMinute As Double
Dim TheSecond As Double
Dim YRDAY As Double
Dim FRDAY As Double
Dim TIM As Double
DAY0 = 726482 ' Day number of Epoch (origin approx. year BC 1)
M0 = 356.6349 ' Mean anomaly at Epoch
L0 = 279.4033 ' Mean longitude along ecliptic
EQC1 = 1.9151 ' Equation of centre harmonic coefficient
EQC2 = 0.02 ' Equation of centre harmonic coefficient
OBLQ = 23.4406 ' Obliquity of ecliptic (inclination)
'Position of the Sun at 11:00 UT on 1997 August 7th
'TheYear = 1997
'TheDay = DatePart("y", "1997-08-7")
'TheHour = 11
'TheMinute = 0
'TheSecond = 0
'LAT 15.0302 North
'LON 204.9828 West
'TheYear = 1985
'TheDay = DatePart("y", "1985-08-12")
'TheHour = 1
'TheMinute = 45
'TheSecond = 0
'LAT 15.8453 North
'LON 057.7893 West
TheYear = 2009
TheDay = DatePart("y", "2009-05-03")
TheHour = 15
TheMinute = 48
TheSecond = 0
' Step 1. Calculate number of days since epoch
'Enter Year (e.g.1985), Day number (Jan 1 = 1), hours, minutes, seconds:
YRDAY = Int((TheYear - 1) * 365.25) - DAY0 ' Days from year difference
FRDAY = (TheHour + ((TheMinute + (TheSecond / 60)) / 60)) / 24 ' Fraction of day
TIM = YRDAY + TheDay + FRDAY ' Time in days since epoch
'Note: These formulae only valid 1900 Mar 01 - 2100 Feb 28.
'Step 2. Compute GHAA and True Longitude
Dim MA As Double
Dim ML As Double
Dim GHAA As Double
Dim TL As Double
MA = M0 + 360 * TIM / 365.2596413 ' Sun's Mean Anomaly
While MA < 0
MA = MA + 360
End While
While MA > 360
MA = MA - 360
End While
ML = L0 + 360 * TIM / 365.2421988 ' Mean Longitude
While ML < 0
ML = ML + 360
End While
While ML > 360
ML = ML - 360
End While
GHAA = ML - 180 + 360 * FRDAY ' Greenwich Hour Angle Aries
While GHAA < 0
GHAA = GHAA + 360
End While
While GHAA > 360
GHAA = GHAA - 360
End While
' True Longitude
TL = ML + EQC1 * Sin(DegreesToRadians(MA)) + EQC2 * Sin(DegreesToRadians(2 * MA))
'Note: 1. 'Longitude' means angle from ascending node up the ecliptic.
' 2. The(Sun) 's argument of perigee (AP) is given by AP = ML - MA
'Step 3a. Compute Unit Vector, Geocentric Rectangular Coordinates
Dim Xsun As Double
Dim Ysun As Double
Dim Zsun As Double
Xsun = Cos(DegreesToRadians(TL))
Ysun = Sin(DegreesToRadians(TL)) * Cos(DegreesToRadians(OBLQ))
Zsun = Sin(DegreesToRadians(TL)) * Sin(DegreesToRadians(OBLQ))
'Step 3b. Unit Vector, Geocentric, Greenwich Coordinates
Dim C As Double
Dim S As Double
Dim Xg As Double
Dim Yg As Double
Dim Zg As Double
C = Cos(DegreesToRadians(GHAA))
S = Sin(DegreesToRadians(GHAA))
Xg = (C * Xsun) + (S * Ysun)
Yg = (-S * Xsun) + (C * Ysun)
Zg = Zsun
'Step 4a. Astronomic Coordinates
Dim DEC As Double
Dim RA As Double
DEC = RadiansToDegrees(ArcSin(Zsun)) ' Declination
RA = ArcTan(DegreesToRadians(Ysun) / DegreesToRadians(Xsun)) ' Right Ascension
If Xsun < 0 Then RA = RA + 180
If RA < 0 Then RA = RA + 360
'ATN is assumed to yield a value in the range +- 90 degrees. RA is in the
'same quadrant as true longitude, TL.
'Step 4b. Geographic Coordinates
Dim LAT As Double
Dim LON As Double
LAT = DEC ' Sun's Latitude +North
LON = GHAA - RA ' Sun's Longitude, +West
If LAT >= 0 Then
'text15.Text = Format(LAT, "##0.0000") & "N"
text15.Text = LAT.ToString("##0.0000") & "N"
Else
'text15.Text = Format(Abs(LAT),"##0.0000") & "S"
text15.Text = Abs(LAT).ToString("##0.0000") & "S"
End If
If LON >= 0 Then
'text16.Text = Format(LON,"000.0000") & "W"
text16.Text = LON.ToString("000.0000") & "W"
Else
'text16.Text = Format(Abs(LON),"000.0000") & "E"
text16.Text = Abs(LON).ToString("000.0000") & "E"
End If
Dim DAY0 As Double
Dim M0 As Double
Dim L0 As Double
Dim EQC1 As Double
Dim EQC2 As Double
Dim OBLQ As Double
Dim TheDate As Date
Dim TheYear As Double
Dim TheDay As Double
Dim TheHour As Double
Dim TheMinute As Double
Dim TheSecond As Double
Dim DaysBefore As Double
Dim YRDAY As Double
Dim FRDAY As Double
Dim TIM As Double
DAY0 = 726482 ' Day number of Epoch (origin approx. year BC 1)
M0 = 356.6349 ' Mean anomaly at Epoch
L0 = 279.4033 ' Mean longitude along ecliptic
EQC1 = 1.9151 ' Equation of centre harmonic coefficient
EQC2 = 0.02 ' Equation of centre harmonic coefficient
OBLQ = 23.4406 ' Obliquity of ecliptic (inclination)
'Position of the Sun at 11:00 UT on 1997 August 7th
TheDate = "1997-08-07 11:00:00"
TheYear = 1997
TheDay = DatePart("y", "1997-08-7")
TheHour = 11
TheMinute = 0
TheSecond = 0
'LAT 15.0302 North
'LON 204.9828 West
'TheDate = "1985-08-12 01:45:00"
'TheYear = 1985
'TheDay = DatePart("y", "1985-08-12")
'TheHour = 1
'TheMinute = 45
'TheSecond = 0
'LAT 15.8453 North
'LON 057.7893 West
'TheDate = "2009-05-03 15:48:00"
'TheYear = 2009
'TheDay = DatePart("y", "2009-05-03")
'TheHour = 15
'TheMinute = 48
'TheSecond = 0
' Step 1. Calculate number of days since epoch
'Enter Year (e.g.1985), Day number (Jan 1 = 1), hours, minutes, seconds:
YRDAY = Int((TheYear - 1) * 365.25) - DAY0 ' Days from year difference
FRDAY = (TheHour + ((TheMinute + (TheSecond / 60)) / 60)) / 24 ' Fraction of day
TIM = YRDAY + TheDay + FRDAY ' Time in days since epoch
DaysBefore = 367 * Year(TheDate) - 7 * (Year(TheDate) + (Month(TheDate) + 9) \ 12) \ 4 + 275 * Month(TheDate) \ 9 + Microsoft.VisualBasic.Day(TheDate) - 730531.5 + Hour(TheDate) / 24
'Note: These formulae only valid 1900 Mar 01 - 2100 Feb 28.
'Step 2. Compute GHAA and True Longitude
Dim MA As Double
Dim ML As Double
Dim GHAA As Double
Dim TL As Double
MA = M0 + 360 * TIM / 365.2596413 ' Sun's Mean Anomaly
While MA < 0
MA = MA + 360
End While
While MA > 360
MA = MA - 360
End While
ML = L0 + 360 * TIM / 365.2421988 ' Mean Longitude
While ML < 0
ML = ML + 360
End While
While ML > 360
ML = ML - 360
End While
GHAA = ML - 180 + 360 * FRDAY ' Greenwich Hour Angle Aries
While GHAA < 0
GHAA = GHAA + 360
End While
While GHAA > 360
GHAA = GHAA - 360
End While
' True Longitude
TL = ML + EQC1 * Sin(DegreesToRadians(MA)) + EQC2 * Sin(DegreesToRadians(2 * MA))
'Note: 1. 'Longitude' means angle from ascending node up the ecliptic.
' 2. The(Sun) 's argument of perigee (AP) is given by AP = ML - MA
'Step 3a. Compute Unit Vector, Geocentric Rectangular Coordinates
Dim Xsun As Double
Dim Ysun As Double
Dim Zsun As Double
Xsun = Cos(DegreesToRadians(TL))
Ysun = Sin(DegreesToRadians(TL)) * Cos(DegreesToRadians(OBLQ))
Zsun = Sin(DegreesToRadians(TL)) * Sin(DegreesToRadians(OBLQ))
'Step 3b. Unit Vector, Geocentric, Greenwich Coordinates
Dim C As Double
Dim S As Double
Dim Xg As Double
Dim Yg As Double
Dim Zg As Double
C = Cos(DegreesToRadians(GHAA))
S = Sin(DegreesToRadians(GHAA))
Xg = (C * Xsun) + (S * Ysun)
Yg = (-S * Xsun) + (C * Ysun)
Zg = Zsun
'Step 4a. Astronomic Coordinates
Dim DEC As Double
Dim RA As Double
DEC = RadiansToDegrees(ArcSin(Zsun)) ' Declination
RA = ArcTan(DegreesToRadians(Ysun) / DegreesToRadians(Xsun)) ' Right Ascension
If Xsun < 0 Then RA = RA + 180
If RA < 0 Then RA = RA + 360
'ATN is assumed to yield a value in the range +- 90 degrees. RA is in the
'same quadrant as true longitude, TL.
'Step 4b. Geographic Coordinates
Dim LAT As Double
Dim LON As Double
LAT = DEC ' Sun's Latitude +North
LON = GHAA - RA ' Sun's Longitude, +West
If LAT >= 0 Then
'text15.Text = Format(LAT, "##0.0000") & "N"
text15.Text = LAT.ToString("##0.0000") & "N"
Else
'text15.Text = Format(Abs(LAT),"##0.0000") & "S"
text15.Text = Abs(LAT).ToString("##0.0000") & "S"
End If
If LON >= 0 Then
'text16.Text = Format(LON,"000.0000") & "W"
text16.Text = LON.ToString("000.0000") & "W"
Else
'text16.Text = Format(Abs(LON),"000.0000") & "E"
text16.Text = Abs(LON).ToString("000.0000") & "E"
End If
Dim LST As Double
Dim HA As Double
Dim observer_lat As Double
Dim observer_lon As Double
' Birmingham UK, you might get this from user entered text fields or GPS
observer_lat = 52.5
observer_lon = -1.91667
'LST = 280.46061837 + 360.98564736629 * -877.04167 + observer_lon
LST = 280.46061837 + 360.98564736629 * DaysBefore + observer_lon
While LST < 0
LST = LST + 360
End While
While LST > 360
LST = LST - 360
End While
HA = LST - RA
' Instructions say to do this, but then they do not do it
'While HA < 0
' HA = HA + 360
'End While
'While HA > 360
' HA = HA - 360
'End While
Dim Alt As Double
Alt = Sin(DegreesToRadians(DEC)) * Sin(DegreesToRadians(observer_lat)) + Cos(DegreesToRadians(DEC)) * Cos(DegreesToRadians(observer_lat)) * Cos(DegreesToRadians(HA))
Alt = RadiansToDegrees(ArcSin(Alt))
Dim X1 As Double
Dim Y1 As Double
Dim AZ As Double
Y1 = -Cos(DegreesToRadians(DEC)) * Cos(DegreesToRadians(observer_lat)) * Sin(DegreesToRadians(HA))
X1 = Sin(DegreesToRadians(DEC)) - Sin(DegreesToRadians(observer_lat)) * Sin(DegreesToRadians(Alt))
AZ = ArcTan(Y1 / X1)
If X1 < 0 Then
AZ = AZ + 180
ElseIf Y1 < 0 And X1 > 0 Then
AZ = AZ + 360
Else
AZ = AZ
End If
http://www.usno.navy.mil/USNO/astronomical-applications
http://www.usno.navy.mil/USNO/astronomical-applications/software-products/mica
Software algorithms (in C and Fortran):
http://www.usno.navy.mil/USNO/astronomical-applications/software-products/novas
And check out the Standards Of Fundamental Astronomy (SOFA) Libraries:
http://www.iau-sofa.rl.ac.uk/2009_0201.html