Attribute VB_Name = "Module4"
Option Explicit

' This module contains solar and lunar position related code, segregated into
' a separate module for programming purposes.  It is dependent on the helper functions
' in Module 1.
'    K. Fisher 5/2006 fisherka@csolutions.net

Public Function PositionMoonCoords(ByVal Year_tmp As Integer, ByVal Month_tmp As Integer, ByVal Day_tmp As Integer, ByVal UTC_Hour_tmp As Integer, _
 ByVal UTC_Min_tmp As Integer, ByVal UTC_Sec_tmp As Double, ByVal IsGregorian_tmp As Integer) As String
' Returns the RA and Dec of the Moon from a UTC date in J2000.0
' Returns other Moon ephemeris characteristics

' Source: Per 1999 _Astronomical Almanac_
' Source accuracy: 0.°11 (6.6 arcmins) with a maximum error of 0.°35 (27 arcmins)

' Accuracy: Against the JPL/NASA ephemeris << http://ssd.jpl.nasa.gov/horizons.cgi >>,
' this alogrithm returns positions within ~approx. 10 arcminutes of the the JPL webapplet
' These errors where not statistically evaluated.

' ACCURACY WARNING: Distance and dependant variable apparent diameter are approx. 0.7%
' different from NASA/JPL values

' Uses TimeJulianDay, TimeJulianCenturiesPeriod, Convert_DecDeg2RADMS and minor trig utilities

' Output:
'   1) RA of Moon without atmospheric refract  ' alf in decimal degrees
'   2) Dec of Moon w/o atmos refraction        ' delta in decimal degrees
'   3) RA of Moon without atmospheric refract  ' alf_hms as hms string in "HHMMSS.SS" format
'   4) Dec of Moon w/o atmos refraction        ' delta_dms as hms string in "sDDMMSS.S" format
'   5) Eclipitic Longitude of Moon             ' gamma in decimal degrees
'   6) Eclipitic Latitude of Moon              ' beta in decimal degrees
'   7) Moon-Earth distance in kilometers - geocentric ' R_kms in kilometer
'   8) Moon apparent diameter surface in arcsecs       ' Moon_dia_arcs_surface
'   9) Moon apparent diameter geocentric in arcsecs       ' Moon_dia_arcs_geocentric
'  10) X ecliptic geocentric coords of Moon in kilometers      ' s_x_km
'  11) Y ecliptic geocentric of Moon in km      ' s_y_km
'  12) Z ecliptic geocentric of Moon in km      ' s_z_km
'  13) Date time string                       ' dtBD
'   13)  Date
'   14)  Hour part
'   15)  AM or PM

' Use
' PositionMoonCoords(2006,5,1,0,0,0,1) give the partial string
' 80.0126573316393 28.2389036096722 052003.04 +281420.1 . . .
' or a JCoords of 052003.04 +281420.1
'    or
' GetCSWord(PositionMoonCoords(2006,5,1,0, 0, 0, 1), 3)
' 052003.04
' GetCSWord(PositionMoonCoords(2006,5,1,0, 0, 0, 1), 4)
' +281420.1

' Testing
'               Duffett-Smith, P.  1988.  Practical Astronomy with
'                    your calculator. 3rd ed. Cambridge Univ.
'                    Press. p. 88 gives for Epoch 1990, position of the Sun
'                    on July 27, 1980 at 00:00:00 as 8h25m46s +19d13m48"


'   Against the JPL/NASA Horizons Webapplet
' PositionMoonCoords(2006,5,1,0,0,0,1) give the partial string
' 80.0126573316393 28.2389036096722 052003.04 +281420.1 . . .
' or a JCoords of 052003.04 +281420.1
' JPL/NASA applet gives 052135.88 +281326.9

' Initialize variables
    Dim A As Double                         ' intermediate variable
    Dim B As Double                         ' intermediate variable
    Dim intY, intM, intD, intH, intMin As Integer
    Dim intY2, intM2, intD2, intH2, intMin2 As Integer  ' aliased working variables
    Dim IsGregorian As Integer
    Dim bolGregorian As Boolean
    Dim dblSec, dblDay As Double     ' Working variables
    Dim dblSec2 As Double                   ' aliased working variable
    Dim S As String                         ' Output string
    
' Constants for Moon code
    Const Earth_radius_km As Double = 6378.14  ' Earth radius (1976 IAU value)
    Const Moon_diameter_km As Double = 3475.6  ' Moon's diameter
    
' Working variables for Moon code
    Dim dblJD As Double                     ' Julian Date
    Dim dtBD As Date                        ' Besslian Day
    Dim intBDoY As Integer                  ' Besslian Day of Year
    Dim t As Double                         ' Julian centuries epoch to J2000.0
    Dim n As Double                         ' number of Julian days between J2000 and now
    Dim g As Double                         ' mean anomaly of the Moon
    Dim gamma As Double                     ' ecliptic's longitude
    Dim beta As Double                      ' ecliptic latitude
    Dim eps As Double                       ' obliquity of the ecliptic
    Dim f, t1 As Double                     ' working variables
    Dim alf As Double                       ' Right Ascension in degrees
    Dim delta As Double                     ' Declination in degrees
    Dim alf_hms As String                   ' Right Ascension as string in HHMMSS.SS format
    Dim delta_dms As String                 ' Declination as string in sDDMMSS.S format
    Dim R_km As Double                      ' Distance from the Moon in kilometers geocentric
    Dim s_x_km, s_y_km, s_z_km As Double    ' equatorial (geocentric) rectangular coordinates of the Moon in kilometers
    Dim Moon_dia_arcs_surface As Double     ' Apparent diameter of the Moon in arcsecs from the surface
    Dim Moon_dia_arcs_geocentric As Double  ' Apparent diameter of the Moon in arcsecs geocentric
    Dim pi_hp As Double                     ' Horizontal parallax of Moon
    
' Working variables for polynomial evaluation
    Dim Ag0, Ag1, Ag2, Ag3, Ag4, Ag5, Ag6   ' accumulators for polymonials for ecliptic longitude gamma
    Dim Ab0, Ab1, Ab2, Ab3                  ' accumulators for polymonials for ecliptic latitude  beta
    Dim Ap0, Ap1, Ap2, Ap3, Ap4             ' accumulators for polymonials for horizontal parallax pi_hp
    Dim Eg1, Eg2, Eg3, Eg4, Eg5, Eg6        ' Evection polymonials for ecliptic longitude gamma
    Dim Eb0, Eb1, Eb2, Eb3                  ' Evection for polymonials for ecliptic latitude  beta
    Dim Ep1, Ep2, Ep3, Ep4                  ' Evection for polymonials for horizontal parallax pi_hp

' Transfer interface input
' Get M/D/Y UTC-Hour/UTC-Min/UTC-Sec

intY = Year_tmp
intM = Month_tmp
intD = Day_tmp
intH = UTC_Hour_tmp
intMin = UTC_Min_tmp
dblSec = UTC_Sec_tmp
IsGregorian = IsGregorian_tmp
Select Case IsGregorian
    Case 0
        bolGregorian = False
    Case 1
        bolGregorian = True
    Case Else
        bolGregorian = True
End Select

' Compute the Julian Date
    ' Alias names to reduce runtime conflict
intY2 = intY
intM2 = intM
intD2 = intD
intH2 = intH
intMin2 = intMin
dblSec2 = dblSec

    ' Get the current Julian date
S = TimeJulianDay(intY2, intM2, intD2, intH2, intMin2, dblSec2, 1)
dblJD = CDbl(GetCSWord(S, 1))

' Compute n, number of days between J2000.0 at noon and now
'      Programming note: Use direct computation of n instead of
'      algorithm method
' Besselian Day of Year
S = CStr(intM) & "/" & CStr(intD) & "/" & CStr(intY)
dtBD = DateValue(S)
'  Note seconds are truncated
dtBD = dtBD + TimeValue(CStr(intH) & ":" & CStr(intMin) & ":" & CStr(dblSec))
intBDoY = DatePart("y", dtBD)
'  Clear working variable
S = ""

'      Get J2000.0 Julian day, which starts at Julian day noon, not midnight
S = TimeJulianDay(2000, 1, 1, 12, 0, 0, 1)
A = CDbl(GetCSWord(S, 1))
n = dblJD - A
'  Clear working variables
S = ""
A = 0

' Compute t, the Julian Centuries between epoch to date and J2000.0
t = TimeJulianCenturiesPeriod(intY2, intM2, intD2, intH2, intMin2, dblSec2, 2000, 1, 1, 12, 0, 0, 1)

' Find the lunar ecliptic longitude gamma
'      Compute the E polynomials bringing them within a 360 deg range
Eg1 = Adj360Rotation((477198.85 * t) + 134.9)
Eg2 = Adj360Rotation((-413335.38 * t) + 259.2)
Eg3 = Adj360Rotation((890534.23 * t) + 235.7)
Eg4 = Adj360Rotation((954397.7 * t) + 269.9)
Eg5 = Adj360Rotation(((35999.05 * t) + 357.5))
Eg6 = Adj360Rotation(((966404.05 * t) + 186.6))

'      Compute the A polynomials
Ag0 = ((481267.883 * t) + 218.32)
Ag1 = 6.29 * Sin(Deg2Rad(Eg1))
Ag2 = -1.27 * Sin(Deg2Rad(Eg2))
Ag3 = 0.66 * Sin(Deg2Rad(Eg3))
Ag4 = 0.21 * Sin(Deg2Rad(Eg4))
Ag5 = -0.19 * Sin(Deg2Rad(Eg5))
Ag6 = -0.11 * Sin(Deg2Rad(Eg6))

'      Sum the A polynomials, bringing the result within a 360 deg range
gamma = Adj360Rotation(Ag0 + Ag1 + Ag2 + Ag3 + Ag4 + Ag5 + Ag6)

'      If the raw gamma is less than 0 degs, add 360 degs
If gamma < 0 Then
   gamma = gamma + 360#
End If

' Find the lunar ecliptic latitude beta
'      Compute the E polynomials bringing them within a 360 deg range
Eb0 = Adj360Rotation((483202.03 * t) + 93.3)
Eb1 = Adj360Rotation((960400.87 * t) + 228.2)
Eb2 = Adj360Rotation((6003.18 * t) + 318.3)
Eb3 = Adj360Rotation((-407332.2 * t) + 217.6)

'      Compute the A polynomials
Ab0 = 5.13 * Sin(Deg2Rad(Eb0))
Ab1 = 0.28 * Sin(Deg2Rad(Eb1))
Ab2 = -0.28 * Sin(Deg2Rad(Eb2))
Ab3 = -0.17 * Sin(Deg2Rad(Ep3))

'      Sum the A polynomials, bringing the result within a 360 deg range
beta = Ab0 + Ab1 + Ab2 + Ab3

' Find the lunar horizontal parallax pi_hp
'      Compute the E polynomials bringing them within a 360 deg range
Ep1 = Adj360Rotation((477198.85 * t) + 134.9)
Ep2 = Adj360Rotation((-413335.38 * t) + 259.2)
Ep3 = Adj360Rotation((890534.23 * t) + 235.7)
Ep4 = Adj360Rotation((954397.7 * t) + 269.9)

'      Compute the A polynomials
Ap0 = 0.9508
Ap1 = 0.0518 * Cos(Deg2Rad(Ep1))
Ap2 = 0.0095 * Cos(Deg2Rad(Ep2))
Ap3 = 0.0078 * Cos(Deg2Rad(Ep3))
Ap4 = 0.0028 * Cos(Deg2Rad(Ep4))

'      Sum the A polynomials, bringing the result within a 360 deg range
pi_hp = Adj360Rotation(Ap0 + Ap1 + Ap2 + Ap3 + Ap4)

' ACCURACY WARNING: Distance and dependant variable apparent diameter are approx. 0.7%
' different from NASA/JPL values

' Find the distance to the Moon - geocentric
R_km = Earth_radius_km / Sin(Deg2Rad(pi_hp))

' Find the apparent size of the Moon in arcsecs at the surface and at the geocenter
'   Convert geocentric to surface of Earth distance
A = R_km - Earth_radius_km
Moon_dia_arcs_surface = ScaleHalfAngleFormula(Moon_diameter_km, A)
Moon_dia_arcs_geocentric = ScaleHalfAngleFormula(Moon_diameter_km, R_km)
' Clear working variables
A = 0

' Find the Right Ascension and Declination at Epoch to Date
'   Do the preliminary computations of factors and the obliquity of the ecliptic
'      Compute eps, obliquity of the ecliptic
eps = 23.439 - (0.0000004 * n)

' Programming note: First version ecliptic to equatorial adjusted
' gamma to alf was empirically adjusted from
' Allen's Astrophysical Quantities (2000) at 668-669 algorithm to convert
' instead of implementing the usual ecliptic to equatorial spherical
' triangle solution
' Decided to not use in favor of traditional Duffet Smith conversion utility
'      Compute f, a working variable to find RA
' f = 180# / PI()
'      Compute t1, a working variable to find RA
' t1 = (Tan(Deg2Rad(eps / 2#))) ^ 2
'      Compute alf, Right Ascension, sin method
' A = (f * t1 * Sin(Deg2Rad(2# * gamma)))
' B = (f / 2#) * (t1 ^ 2) * Sin(Deg2Rad(4# * gamma))
' - Cos(beta) adjusted
' alf = gamma - A + B - Cos(Deg2Rad(beta))
'      Put alf in the range 0 =< alf < 360
' alf = Adj360Rotation(alf)
'      Clear working variables
' A = 0
' B = 0

' Compute alf, Right Ascension, basic Duffet-Smith method at p. 40, sec. 27
A = ((Sin(Deg2Rad(gamma)) * Cos(Deg2Rad(eps))) - (Tan(Deg2Rad(beta)) * Sin(Deg2Rad(eps))))
B = Cos(Deg2Rad(gamma))
alf = Adj360Rotation(Rad2Deg(Atan2(A, B)))
'      Clear working variables
A = 0
B = 0

'      Compute delta, eclipitc beta to Declination
A = Sin(Deg2Rad(eps)) * Sin(Deg2Rad(gamma)) * Cos(Deg2Rad(beta))
B = Sin(Deg2Rad(beta)) * Cos(Deg2Rad(eps))
delta = Rad2Deg(ArcSin(A + B))
'      Clear working variables
A = 0
B = 0

'     Prepare string reformatted RA and Dec
alf_hms = Convert_DecDeg2RADMS(alf, "HHMMSS.SS")
delta_dms = Convert_DecDeg2RADMS(delta, "sDDMMSS.S")

' Convert s_x, s_y, s_z, the rectangular coordinates of the Moon in kilometers
s_x_km = R_km * Cos(Deg2Rad(gamma)) * Cos(Deg2Rad(beta))
'      Clear working variables
A = 0
B = 0

' Find y coord
A = Cos(Deg2Rad(eps)) * Sin(Deg2Rad(gamma)) * Cos(Deg2Rad(beta))
B = Sin(Deg2Rad(beta)) * Sin(Deg2Rad(eps))
s_y_km = R_km * (A - B)
'      Clear working variables
A = 0
B = 0

' Find z coord
A = Sin(Deg2Rad(eps)) * Sin(Deg2Rad(gamma)) * Cos(Deg2Rad(beta))
B = Sin(Deg2Rad(beta)) * Cos(Deg2Rad(eps))
s_z_km = R_km * (A + B)
'      Clear working variables
A = 0
B = 0

' Assemble result string
'   Clear working "S" string
S = ""

'   1) RA of Moon without atmospheric refract  ' alf in decimal degrees
S = Trim(CStr(alf))

'   2) Dec of Moon w/o atmos refraction        ' delta in decimal degrees
S = S & " " & Trim(CStr(delta))

'   3) RA of Moon without atmospheric refract  ' alf_hms as hms string in "HHMMSS.SS" format
S = S & " " & Trim(alf_hms)

'   4) Dec of Moon w/o atmos refraction        ' delta_dms as hms string in "sDDMMSS.S" format
S = S & " " & Trim(delta_dms)

'   5) Eclipitic Longitude of Moon             ' gamma in decimal degrees
S = S & " " & Trim(CStr(gamma))

'   6) Eclipitic Latitude of Moon              ' beta in decimal degrees
S = S & " " & Trim(CStr(beta))

'   7) Moon-Earth distance in kilometers - geocentric ' R_kms in kilometer
S = S & " " & Trim(CStr(R_km))

'   8) Moon apparent diameter surface in arcsecs       ' Moon_dia_arcs_surface
S = S & " " & Trim(CStr(Moon_dia_arcs_surface))

'   9) Moon apparent diameter geocentric in arcsecs       ' Moon_dia_arcs_geocentric
S = S & " " & Trim(CStr(Moon_dia_arcs_geocentric))

'  10) X ecliptic geocentric coords of Moon in kilometers      ' s_x_km
S = S & " " & Trim(CStr(s_x_km))

'  11) Y ecliptic geocentric of Moon in km      ' s_y_km
S = S & " " & Trim(CStr(s_y_km))

'  12) Z ecliptic geocentric of Moon in km      ' s_z_km
S = S & " " & Trim(CStr(s_z_km))

'  13) Date time string                       ' dtBD
'   13)  Date
'   14)  Hour part
'   15)  AM or PM

'   1) Date string                            ' dtBD
'       Programming note: Correct runtime error
If intH = 0 And intMin = 0 And dblSec = 0# Then
    S = S & " " & CStr(dtBD) & " 00:00:00 AM"
Else
    S = S & " " & CStr(dtBD)
End If

' Return result
PositionMoonCoords = S

' Debug.Print "gamma", gamma
' Debug.Print "beta", beta
' Debug.Print "pi_hp", pi_hp
' Debug.Print "R_km", R_km
' Debug.Print "Moon_dia_arcs_surface", Moon_dia_arcs_surface
' Debug.Print "Moon_dia_arcs_geocentric ", Moon_dia_arcs_geocentric
' Debug.Print "alf", alf
' Debug.Print "delta", delta
' Debug.Print "RA-Dec", alf_hms, delta_dms
' Debug.Print "s_x_km", s_x_km
' Debug.Print "s_y_km", s_y_km
' Debug.Print "s_z_km", s_z_km

End Function

Public Function PositionSolarCoords(ByVal Year_tmp As Integer, ByVal Month_tmp As Integer, ByVal Day_tmp As Integer, ByVal UTC_Hour_tmp As Integer, _
 ByVal UTC_Min_tmp As Integer, ByVal UTC_Sec_tmp As Double, ByVal IsGregorian_tmp As Integer) As String
' Returns the RA and Dec of the Sun from a UTC date and imte between 1950 and 2050

' Source: Allen's Astrophysical Quantities (2000) at 668-669 after _Astronomical Almanac_
' Source accuracy: "The following formulas . . . give the apparent coordinates of the Sun to a precision of
' 0.°01 (6 arcmins) and the equation of time to a precision of 0m.1 between 1950 and 2050 . . ."

' Accuracy: Against the JPL/NASA ephemeris << http://ssd.jpl.nasa.gov/horizons.cgi >>,
' this alogrithm returns positions within 1 arcminute or 60" of the the JPL webapplet

' WARNING: Except for the Ecliptic longitude and its dependent variables like geocentric
' coordinates of the Sun, this implementation appears reasonably accurate.  For a higher
' precision ecliptic longitude, add one hour to the time.

' Uses TimeJulianDay, Convert_DecDeg2RADMS and minor trig utilities

' Output:
'   1) RA of Sun without atmospheric refract  ' alf in decimal degrees
'   2) Dec of Sun w/o atmos refraction        ' delta in decimal degrees
'   3) RA of Sun without atmospheric refract  ' alf_hms as hms string in "HHMMSS.SS" format
'   4) Dec of Sun w/o atmos refraction        ' delta_dms as hms string in "sDDMMSS.S" format
'   5) Eclipitic Longitude of Sun             ' gamma in decimal degrees
'   6) Eclipitic Latitude of Sun              ' beta in decimal degrees
'   7) Mean longitude of the Sun              ' L in decimal degrees
'   8) Mean anomaly of the Sun                ' g in decimal degrees
'   9) Sun-Earth distance in AU               ' R_au in AU
'  10) Sun-Earth distance in kilometers       ' R_au in kilometers
'  11) Sun-Earth distance in light-minutes    ' R_au in light-minutes
'  12) Sun apparent diameter in arcsecs       ' Sun_dia_arcsecs
'  13) Equation of Time value in minutes      ' E in minutes
'  14) X Geocentric coords of Sun in AU       ' s_x
'  15) Y Geocentric coords of Sun in AU       ' s_y
'  16) Z Geocentric coords of Sun in AU       ' s_z
'  17) X Geocentric coords of Sun in km       ' s_x_km
'  18) Y Geocentric coords of Sun in km       ' s_y_km
'  19) Z Geocentric coords of Sun in km       ' s_z_km
'  20) Date time string                       ' dtBD
'   20)  Date
'   21)  Hour part
'   22)  AM or PM

' Use
' ? PositionSolarCoords(1980,7,27,0,0,0,1) give the partial string
' 126.573814005204 19.1995317049426 082617.72 +191158.3 . . .
' or a JCoords of 082617.72 +191158.3
'    or
' GetCSWord(PositionSolarCoords(1980, 7, 27, 0, 0, 0, 1), 3)
' 082617.72
' GetCSWord(PositionSolarCoords(1980, 7, 27, 0, 0, 0, 1), 4)
' +191158.3

' Testing
'               Duffett-Smith, P.  1988.  Practical Astronomy with
'                    your calculator. 3rd ed. Cambridge Univ.
'                    Press. p. 88 gives for Epoch 1990, position of the Sun
'                    on July 27, 1980 at 00:00:00 as 8h25m46s +19d13m48"
' ? PositionSolarCoords(1980,7,27,0,0,0,1) give the partial string
' 126.573814005204 19.1995317049426 082617.72 +191158.3 . . .
'                    or a JCoords of 082617.72 +191158.3
'                    which is within 30" of Duffet-Smith's result
' JPL/NASA ephemeris << http://ssd.jpl.nasa.gov/horizons.cgi >> for
' 2006-04-13 00:00:00 UTC gives:
'    RA/DEC   01 24 32.92 +08 53 06.7
'    Sun diameter  1914.394
'    Sun dist AU 1.00255724667839
'    Ecliptic Long/Lat 22.9311258  -0.0000783
' This algorithm code gives:
' ? GetCSWord(PositionSolarCoords(2006, 4, 13, 0, 0, 0, 1), 3)
'    012440.81  ' or ~ 8" RA difference
' ? GetCSWord(PositionSolarCoords(2006, 4, 13, 0, 0, 0, 1), 4)
'    +085351.9  ' or ~ 46" Dec difference
' ? GetCSWord(PositionSolarCoords(2006,4,13,0,0,0,1),5)
'  22.883521781257 ' or ~ 0.05 degs or 2 arcmin Eclipitic longitude difference
' Adding one hour to the time index gives a more accurate Ecliptic longitude
' GetCSWord(PositionSolarCoords(2006, 4, 13, 1, 0, 0, 1), 5)
'  22.924123329028 ' or ~ 0.007 degs or ~ 25"
' ? GetCSWord(PositionSolarCoords(2006,4,13,0,0,0,1),9)
'  1.00260081687329 ' or ~ 0.000043 AU difference in RA_distance
' ? GetCSWord(PositionSolarCoords(2006,4,13,0,0,0,1),12)
'  1914.54063042379 ' or ~ 0.146 arcsec difference in Sun diameter
' An altitude and az check for Salt Lake City, Utah,
' JPL/NASA gives Az/Alt 262.2359  22.4714
' ? Equ2Horizon_coord_convert(GetCSWord(PositionSolarCoords(2006,4,13,1,0,0,1),1),GetCSWord(PositionSolarCoords(2006,4,13,1,0,0,1),2),"00:00:00.00", #4/13/2006#, -111.890833, 40.756389)
' 262.236324389033 22.4676321504436
'   or an az difference of ~ 0.0004 degs or ~ 1.44"
'   and an alt difference of ~ 0.004 degs or ~14.4"

' Initialize variables
    Dim A As Double                         ' intermediate variable
    Dim B As Double                         ' intermediate variable
    Dim intY, intM, intD, intH, intMin As Integer
    Dim intY2, intM2, intD2, intH2, intMin2 As Integer  ' aliased working variables
    Dim IsGregorian As Integer
    Dim bolGregorian As Boolean
    Dim dblSec, dblDay As Double     ' Working variables
    Dim dblSec2 As Double                   ' aliased working variable
    Dim S As String                         ' Output string
    
' Constants
    Const kAU2km As Double = 149597870.691   ' 1 AU= 149597870.691 km
    Const c_light_kms As Double = 299792.458 ' speed of light km/s c= 299792.458 km/s
    Const c_light_kmm = (299792.458 * 60#)   ' speed of light km/s c= 299792.458 km/s * 60
'   Not used stored
    Const parallax_horiz_deg As Double = 0.00244281888  ' Horizontal parallax in degrees
    Const parallax_horiz_arcsec As Double = 8.794148    ' Horizontal parallax in arcsecs
    Const Light_time_days As Double = 0.0058            ' Light time distance in light days

' Working variables for Astronomical Almanac algorithm code
    Dim dblJD As Double                     ' Julian Date
    Dim dtBD As Date                        ' Besslian Day
    Dim intBDoY As Integer                  ' Besslian Day of Year
    Dim dtfracBD As Double                  ' fraction of day from 0 UTC
    Dim n As Double                         ' number of Julian days between J2000 and now
    Dim L As Double                         ' mean longitude of the Sun
    Dim g As Double                         ' mean anomaly of the Sun
    Dim gamma As Double                     ' ecliptic's longitude
    Dim beta As Double                      ' ecliptic latitude
    Dim eps As Double                       ' obliquity of the ecliptic
    Dim f, t As Double                      ' working variables
    Dim alf As Double                       ' Right Ascension in degrees
    Dim delta As Double                     ' Declination in degrees
    Dim alf_hms As String                   ' Right Ascension as string in HHMMSS.SS format
    Dim delta_dms As String                 ' Declination as string in sDDMMSS.S format
    Dim R_au As Double                      ' Distance from the Sun in AU
    Dim R_km As Double                      ' Distance from the Sun in kilometers
    Dim R_lmin As Double                    ' Distance from the Sun in light minutes
    Dim s_x, s_y, s_z                       ' equatorial (geocentric) rectangular coordinates of the Sun in AU
    Dim s_x_km, s_y_km, s_z_km              ' equatorial (geocentric) rectangular coordinates of the Sun in AU
    Dim E As Double                         ' Value of the Equation of Time in minutes
    Dim Sun_dia_deg                         ' Apparent diameter of the Sun in degrees
    Dim Sun_dia_arcsecs                     ' Apparent diameter of the Sun in arcsecs

' Transfer interface input
' Get M/D/Y UTC-Hour/UTC-Min/UTC-Sec

intY = Year_tmp
intM = Month_tmp
intD = Day_tmp
intH = UTC_Hour_tmp
intMin = UTC_Min_tmp
dblSec = UTC_Sec_tmp
IsGregorian = IsGregorian_tmp
Select Case IsGregorian
    Case 0
        bolGregorian = False
    Case 1
        bolGregorian = True
    Case Else
        bolGregorian = True
End Select

' Compute the Julian Date
    ' Alias names to reduce runtime conflict
intY2 = intY
intM2 = intM
intD2 = intD
intH2 = intH
intMin2 = intMin
dblSec2 = dblSec

    ' Get the current Julian date
S = TimeJulianDay(intY2, intM2, intD2, intH2, intMin2, dblSec2, 1)
dblJD = CDbl(GetCSWord(S, 1))

'      Programming note: skip - use direct computation of n
' Besselian Day of Year
S = CStr(intM) & "/" & CStr(intD) & "/" & CStr(intY)
dtBD = DateValue(S)
'  Note seconds are truncated
dtBD = dtBD + TimeValue(CStr(intH) & ":" & CStr(intMin) & ":" & CStr(dblSec))
intBDoY = DatePart("y", dtBD)
'  Clear working variable
S = ""

' Compute the fraction of the day
'      Programming note: skip - use direct computation of n

' Compute n, number of days between J2000.0 at noon and now
'      Programming note: Use direct computation of n instead of
'      algorithm method

' Programming runtine accuracy note:
'   When comparing this to the JPL/NASA ephemeris << http://ssd.jpl.nasa.gov/horizons.cgi >>
'   I noticed that the ecliptic longitude off by 1 hour.
'   RA and Dec are almost exactly on.  If the time is advanced one hour in this algorithm
'   it is almost equal to the JPL/NASA ephemeris

'      Get J2000.0 Julian day, which starts at Julian day noon, not midnight
S = TimeJulianDay(2000, 1, 1, 12, 0, 0, 1)
A = CDbl(GetCSWord(S, 1))
n = dblJD - A
'  Clear working variables
S = ""
A = 0

' Compute L, mean longitude of the Sun corrected for aberration
L = 280.46 + (0.9856474 * n)

' Compute g, mean anomaly of the Sun
g = 357.528 + (0.9856003 * n)

' Put L and g in the range 0 =< L,g < 360
L = Adj360Rotation(L)
g = Adj360Rotation(g)

' Compute gamma, ecliptic longitude
gamma = L + (1.915 * Sin(Deg2Rad(g))) + (0.2 * Sin(Deg2Rad(2 * g)))
' Put gamma in the range 0 =< gamma < 360
gamma = Adj360Rotation(gamma)

' Compute beta, ecliptic latitude
beta = 0

' Compute eps, obliquity of the ecliptic
eps = 23.439 - (0.0000004 * n)

' Compute f, a working variable to find RA
f = 180# / PI()

' Compute t, a working variable to find RA
t = (Tan(Deg2Rad(eps / 2#))) ^ 2

' Compute alf, Right Ascension, sin method
A = (f * t * Sin(Deg2Rad(2# * gamma)))
B = (f / 2#) * (t ^ 2) * Sin(Deg2Rad(4# * gamma))
alf = gamma - A + B
' Put alf in the range 0 =< alf < 360
alf = Adj360Rotation(alf)
' Clear working variables
A = 0
B = 0

' Compute delta, Declination
A = Sin(Deg2Rad(eps)) * Sin(Deg2Rad(gamma))
delta = Rad2Deg(ArcSin(A))
'  Clear working variables
A = 0

' Prepare string reformatted RA and Dec
alf_hms = Convert_DecDeg2RADMS(alf, "HHMMSS.SS")
delta_dms = Convert_DecDeg2RADMS(delta, "sDDMMSS.S")

' Compute R_au, the distance to the Sun
A = (0.01671 * Cos(Deg2Rad(g)))
B = (0.00014 * Cos(Deg2Rad(2# * g)))
R_au = 1.00014 - A - B

' Convert R_au, the distance to the Sun, to one-way kilometers
R_km = R_au * kAU2km

' Convert R_au, the distance to the Sun, to one-way light minutes
R_lmin = R_km / c_light_kmm

' Compute s_x, s_y, s_z, the rectangular coordinates of the Sun in AU
s_x = R_au * Cos(Deg2Rad(gamma))
s_y = R_au * Cos(Deg2Rad(eps)) * Sin(Deg2Rad(gamma))
s_z = R_au * Sin(Deg2Rad(eps)) * Sin(Deg2Rad(gamma))

' Convert s_x, s_y, s_z, the rectangular coordinates of the Sun in AU to kilometers
s_x_km = s_x * kAU2km
s_y_km = s_y * kAU2km
s_z_km = s_z * kAU2km

' Compute E, value of the Equation of Time in minutes
E = 4 * (L - alf)

' Set a constant, parallax_horiz_deg, horizontal parallax of the Sun 1 AU = 8.79 arcsecs
'   See "Constants" above

' Compute Sun_dia_deg, Sun's apparent diameter
'    Allen's alogrithm computes the semidiameter or radius
Sun_dia_deg = 2 * 0.2666 / R_au
Sun_dia_arcsecs = Sun_dia_deg * 3600#

' Compute Light_time_days, the distance to the Sun in light-days
'    Not implemented - see "Convert to Light-minutes" above

' Assemble result string
'   Clear working "S" string
S = ""

'   1) RA of Sun without atmospheric refract  ' alf in decimal degrees
S = Trim(CStr(alf))

'   2) Dec of Sun w/o atmos refraction        ' delta in decimal degrees
S = S & " " & Trim(CStr(delta))

'   3) RA of Sun without atmospheric refract  ' alf_hms as hms string in "HHMMSS.SS" format
S = S & " " & Trim(alf_hms)

'   4) Dec of Sun w/o atmos refraction        ' delta_dms as hms string in "sDDMMSS.S" format
S = S & " " & Trim(delta_dms)

'   5) Eclipitic Longitude of Sun             ' gamma in decimal degrees
S = S & " " & Trim(CStr(gamma))

'   6) Eclipitic Latitude of Sun              ' beta in decimal degrees
S = S & " " & Trim(CStr(beta))

'   7) Mean longitude of the Sun              ' L in decimal degrees
S = S & " " & Trim(CStr(L))

'   8) Mean anomaly of the Sun                ' g in decimal degrees
S = S & " " & Trim(CStr(g))

'   9) Sun-Earth distance in AU               ' R_au in AU
S = S & " " & Trim(CStr(R_au))

'  10) Sun-Earth distance in kilometers       ' R_au in kilometers
S = S & " " & Trim(CStr(R_km))

'  11) Sun-Earth distance in light-minutes    ' R_au in light-minutes
S = S & " " & Trim(CStr(R_lmin))

'  12) Sun apparent diameter in arcsecs       ' Sun_dia_arcsecs
S = S & " " & Trim(CStr(Sun_dia_arcsecs))

'  13) Equation of Time value in minutes      ' E in minutes
S = S & " " & Trim(CStr(E))

'  14) X equatorial geocentric coords of Sun in AU       ' s_x
S = S & " " & Trim(CStr(s_x))

'  15) Y equatorial geocentric of Sun in AU       ' s_y
S = S & " " & Trim(CStr(s_y))

'  16) Z equatorial geocentric of Sun in AU       ' s_z
S = S & " " & Trim(CStr(s_z))

'  17) X equatorial geocentric of Sun in km       ' s_x_km
S = S & " " & Trim(CStr(s_x_km))

'  18) Y equatorial geocentric of Sun in km       ' s_y_km
S = S & " " & Trim(CStr(s_y_km))

'  19) Z equatorial geocentric of Sun in km       ' s_z_km
S = S & " " & Trim(CStr(s_z_km))

'  20) Date time string                       ' dtBD
'   20)  Date
'   21)  Hour part
'   22)  AM or PM

'   1) Date string                            ' dtBD
'       Programming note: Correct runtime error
If intH = 0 And intMin = 0 And dblSec = 0# Then
    S = S & " " & CStr(dtBD) & " 00:00:00 AM"
Else
    S = S & " " & CStr(dtBD)
End If

' Return result
PositionSolarCoords = S

End Function