B4A Library Astro library

Sub SunPosition(  Lat As Double,  Lon As Double,  Timenow As Long ) As Double()
 
        ' Main variables
         Dim  twopi As Double = 2*cPI
        Dim rad  = cPI/180     As Double
        Dim dEarthMeanRadius  =   6371.01    As Double' In km
        Dim dAstronomicalUnit  =  149597890.0 As Double' In km
        Dim dElapsedJulianDays As Double
        Dim dDecimalHours As Double
        Dim dEclipticLongitude As Double
        Dim dEclipticObliquity As Double
        Dim dRightAscension As Double
        Dim dDeclination As Double

        ' Auxiliary variables
        Dim dY As Double
        Dim dX As Double

        ' Calculate difference in days between the current Julian Day
        ' And JD 2451545.0, which Is noon 1 January 2000 Universal Time
     
                ' Calculate time of the day in UT decimal hours
            Dim offsetTicks As Long = DateTime.TimeZoneOffset * DateTime.TicksPerHour
            Dim UTCTimenow As Long = Timenow - offsetTicks
            dDecimalHours = DateTime.GetHour(UTCTimenow ) + DateTime.GetMinute(UTCTimenow )/60 + DateTime.GetSecond(UTCTimenow )/3600.0
            ' Calculate current Julian Day
        Dim    liAux1 As Double = (DateTime.GetMonth(Timenow) - 14) / 12
        Dim    liAux2 As Double = (1461 * (DateTime.GetYear(Timenow) + 4800 + liAux1)) / 4 + _
                                (367*(DateTime.GetMonth(Timenow) _
                                - 2 -12 * liAux1)) / 12- (3 * ((DateTime.GetYear(Timenow) + 4900  _     
                                 + liAux1) / 100)) / 4 + DateTime.GetDayOfMonth(Timenow) - 32075
        Dim    dJulianDate As Double = liAux2 - 0.5 + dDecimalHours / 24.0
            ' Calculate difference between current Julian Day And JD 2451545.0
            dElapsedJulianDays = dJulianDate - 2451545.0
        '    dElapsedJulianDays = Julianday2(UTCTimenow , offsetTicks ) - 2001545.0
        '   dElapsedJulianDays = UTCTimenow  / DateTime.TicksPerDay - 2001545.0  + 2440587.5

        ' Calculate ecliptic coordinates (ecliptic longitude And obliquity of the
        ' ecliptic in radians but without limiting the angle To be less than 2*Pi
        ' (i.e., the result may be greater than 2*Pi)
        '
            Dim  dMeanLongitude As Double
            Dim dMeanAnomaly As Double
            Dim dOmega As Double
            dOmega = 2.1429 - 0.0010394594 * dElapsedJulianDays
            dMeanLongitude = 4.8950630 + 0.017202791698 * dElapsedJulianDays '  Radians
            dMeanAnomaly = 6.2400600 + 0.0172019699 * dElapsedJulianDays
            dEclipticLongitude = dMeanLongitude + 0.03341607 * Sin(( dMeanAnomaly ) _
                    + 0.00034894 * Sin( 2*dMeanAnomaly ) - 0.0001134 - 0.0000203 * Sin(dOmega))
            dEclipticObliquity = 0.4090928 - 6.2140e-9 * dElapsedJulianDays    + 0.0000396 * Cos(dOmega)
     

        ' Calculate celestial coordinates ( right ascension And declination ) in radians
        ' but without limiting the angle To be less than 2*Pi (i.e., the result may be
        ' greater than 2*Pi)
        '
            Dim dSin_EclipticLongitude As Double
            dSin_EclipticLongitude = Sin( dEclipticLongitude )
            dY = Cos( dEclipticObliquity ) * dSin_EclipticLongitude
            dX = Cos( dEclipticLongitude )
            dRightAscension = ATan2( dY,dX )
            If( dRightAscension < 0.0 ) Then dRightAscension = dRightAscension + twopi
            dDeclination = ASin( Sin( dEclipticObliquity ) * dSin_EclipticLongitude )
     

        ' Calculate local coordinates ( azimuth And zenith angle ) in degrees
     
            Dim dGreenwichMeanSiderealTime As Double
            Dim dLocalMeanSiderealTime As Double
            Dim dLatitudeInRadians As Double
            Dim dHourAngle As Double
            Dim dCos_Latitude As Double
            Dim dSin_Latitude As Double
            Dim dCos_HourAngle As Double
            Dim dParallax As Double
            Dim dZenithAngle As Double
            Dim dAzimuth As Double
            dGreenwichMeanSiderealTime = 6.6974243242 + 0.0657098283 * dElapsedJulianDays + dDecimalHours
            dLocalMeanSiderealTime = (dGreenwichMeanSiderealTime * 15 + Lon) * rad
            dHourAngle = dLocalMeanSiderealTime - dRightAscension
            dLatitudeInRadians = Lat * rad
            dCos_Latitude = Cos( dLatitudeInRadians )
            dSin_Latitude = Sin( dLatitudeInRadians )
            dCos_HourAngle= Cos( dHourAngle )
            dZenithAngle = (ACos( dCos_Latitude * dCos_HourAngle *Cos(dDeclination) + Sin( dDeclination ) * dSin_Latitude))
            dY = -Sin( dHourAngle )
            dX = Tan( dDeclination ) * dCos_Latitude - dSin_Latitude * dCos_HourAngle
            dAzimuth = ATan2( dY, dX )
            If ( dAzimuth < 0.0 ) Then     dAzimuth = dAzimuth + twopi
            dAzimuth = dAzimuth / rad
            ' Parallax Correction
            dParallax = (dEarthMeanRadius / dAstronomicalUnit) *Sin(dZenithAngle)
            dZenithAngle = (dZenithAngle + dParallax) / rad
            Dim dbl(2)  As Double
            dbl(0) = dAzimuth
            dbl(1) = 90 - dZenithAngle
            Return dbl

End Sub

Sub GetNoonTime(h() As String) As String
    Dim sf As StringFunctions
    sf.Initialize
    Dim h1 As Int = sf.Left(h(0),h(0).IndexOf(":"))*60+sf.Right(h(0),h(0).IndexOf(":"))
    Dim h2 As Int = sf.Left(h(1),h(1).IndexOf(":"))*60+sf.Right(h(1),h(1).IndexOf(":"))
    Dim h3 As Int = (h1 + h2)/2
    Return NumberFormat(h3/60,2,0) & ":" & NumberFormat(h3 Mod 60,2,0)
End Sub