Android Question Calculating Solar Azimuth [Solved]

[Roger Daley]

Hi All,

I am attempting to calculate the Solar Azimuth in order to rotate a map towards the sun. For almost a week I have been researching the Web, found many, many references but still can't get any to give a valid result.

I have been using two online calculators that agree with each other as my bench mark none of formula I have found agree with these over the day.


1. Yes. I have corrected for day light savings
2. Most articles are written from a northern hemisphere perspective and are slipshop in mentioning if a southern hemisphere correction is needed.
3. Some formula assume an output of 0-360 degress clockwise from Nth. Some assume a +/- divergence from Nth or Sth depending on northern/southern hemisphere.

Has anyone actually used a method which consistently produces a valid result? Sorry if this sounds like I am asking for someone else to do the work but this "simple" task has turned in to something like bare handed eel fishing.

Below, not usable code but some of the formula:

        a9 = (360/365.242199)*(NDay + 10)
        b9 = a9 + 1.9136790357*SinD((360/365.242199)*(NDay-2))                     '1.91... = 360 / cPI * 0.0167... ecc of earth
        c9 = (a9 - ATanD(TanD(b9) / CosD(23.43755447))) / 180
        EOT = 720 * (c9 - Round(c9))
       
        Declination = -ASinD(SinD(23.44)*CosD(b9))
      
        LSTM = 15*(Round(Longitude2/15))
        SolarTime = DeciHour +(4*(Longitude1-LSTM))+ EOT                                         'Time of day in minutes
        HourAngle = (SolarTime-720)/4                                    'Hour Angle in degrees
        AltitudeAngle = ASinD(CosD(Abs(Latitude1))*CosD(Declination)*CosD(HourAngle)+SinD(Abs(Longitude1))*SinD(Declination))
        'Azimuth = ACosD((SinD(AltitudeAngle)*SinD(Latitude1)-SinD(Declination))/(CosD(AltitudeAngle)*CosD(Latitude1)))
        'Azimuth = ACosD(SinD(Declination)-(SinD(AltitudeAngle)*SinD(Latitude1))/(CosD(AltitudeAngle)*CosD(Latitude1)))

Private SZA As Double
'Now we can calculate the Sun Zenith Angle (SZA): 
'Cos(SZA) = Sin(Latitude)*Sin(D)+Cos(Latitude)*Cos(D)*Cos(SHA)
        'SZA = ACosD(SinD(Latitude1)*SinD(Declination)+CosD(Latitude1)*CosD(Declination)*Cos(HourAngle))
'        Azimuth = ACosD((Sin(Declination)-SinD(Latitude1)*CosD(SZA))/(CosD(Latitude1)*Sin(SZA)))
        SZA=90-AltitudeAngle
        Azimuth = ACosD((SinD(Declination)*CosD(Latitude1)-CosD(HourAngle)*CosD(Declination)*SinD(Latitude1))/SinD(SZA))

Regards [Exasperated]
Roger

 

[derez]

https://www.b4x.com/android/forum/threads/astro-library.15975/
Last version (1.3) includes sun position method.

 

[Roger Daley]

Hi Derez,

It looks good but I am having trouble: Below is a test Sub and the resulting error. Any idea what I did wrong?

Sub SolarBearing1
    Private AzimuthAstro, LongitudeAstro, LatitudeAstro As Double
    Private TimeNowAstro As Long
    Private Mylocation1 As LatLng
    Private Astro2 As Astro
   
    Mylocation1 = gmap.MyLocation
    LatitudeAstro = Mylocation1.Latitude
    LongitudeAstro = Mylocation1.Longitude
    TimeNowAstro=DateTime.Now
   
AzimuthAstro = Astro2.SunPosition(LatitudeAstro,LongitudeAstro,TimeNowAstro)   

Log("AzimuthAstro        "&AzimuthAstro)
   
End Sub

B4A version: 6.31
Parsing code.    (0.06s)
Compiling code.    (0.18s)
Compiling layouts code.    (0.02s)
Organizing libraries.    (0.00s)
Generating R file.    (0.28s)
Compiling generated Java code.    Error
B4A line: 2185
AzimuthAstro = Astro2.SunPosition(LatitudeAstro,Lo
javac 1.8.0_65
src\horsetrailer\B4A\AntennaBearingTool\main.java:6197: error: incompatible types: double[] cannot be converted to double
_azimuthastro = (double)(_astro2.SunPosition(_latitudeastro,_longitudeastro,_timenowastro));Debug.locals.put("AzimuthAstro", _azimuthastro);
                        ^
1 error

Regards Roger

 

[DonManfred]

The sub returns an array

AzimuthAstro = Astro2.SunPosition(...)

should be

Private AzimuthAstro() As Double

i guess

 

[derez]

DonManfred is right, the method returns an array of two doubles, azimuth and elevation.

 

[Roger Daley]

Hi Guys.

Still not producing a valid result. To which convention does azimuth figure relate? [0-360, deviation from Nth/Sth]

Sub SolarBearing1
    Private AzimuthAstro, LongitudeAstro, LatitudeAstro As Double
    Private TimeNowAstro As Long
    Private Mylocation1 As LatLng
    Private Astro2 As Astro
    Private SunPosition1() As Double
   
    Mylocation1 = gmap.MyLocation
    LatitudeAstro = Mylocation1.Latitude
    LongitudeAstro = Mylocation1.Longitude
    TimeNowAstro=DateTime.Now
   
    SunPosition1 = Astro2.SunPosition(LatitudeAstro,LongitudeAstro,TimeNowAstro)   
    AzimuthAstro = SunPosition1(0)
    'AzimuthAstro = SunPosition1(1)

Log("AzimuthAstro        "&AzimuthAstro)
   
End Sub

Regards Roger

 

[derez]

Hi Guys.

Still not producing a valid result. To which convention does azimuth figure relate? [0-360, deviation from Nth/Sth]

Azimuth : 0-360 degrees, north is 0, East 90 etc.
Elevation: 0 - 90, 0 is horizon, 90 is Zenith.

What do you get ?

 

[Roger Daley]

Azimuth : 0-360 degrees, north is 0, East 90 etc.
Elevation: 0 - 90, 0 is horizon, 90 is Zenith.

What do you get ?

Hi Derez,


Firstly I think I have being doing something fundementaslly wrong with your Astro library, but don't know what.
I was in debug, using real time [AM] and getting results between 0-180, but all results were different to the online calculator at http://www.jgiesen.de/astro/suncalc/ for the same inputs.

I have since adapted the code you pasted here [post#7 SunAzEL.zip]. It appears to work perfectly but I have only used it [PM] both in debug and Release mode out on the street. I will continue to test in the real world tomorrow morning. If this continues to work I will try to emulate eastern longitudes and northern latitudes.
I will let you know the outcome.

regards Roger

PS The link in #6 in the other thread no longer works. I think this is the new link.

 

[derez]

Hi Derez,


Firstly I think I have being doing something fundementaslly wrong

The algorithm from that link is the one implemented in Astro library.
If you use southern coordinates use negative values for lat.

Why don't you show the code you use ? It will be easier to help !

 

[Roger Daley]

It being the same algorithm is why I figured I was doing something wrong with Astro.

Using that algorithm appears to be working AM and PM. I am happy with the result but I am still curious why I couldn't use Astro.
Many thanks for the code.
I have attached a zip of the project, I have left in previous attempts of code commented out.
Any feedback welcome.
Note: The "Sun Dial" mode is a last ditch method for the App.

Roger

 

[Roger Daley]

To close out the question this is the code I finally used. Many thanks to Derez.

Sub Process_Globals
    Private DegreesSDRoation As Float
    Private SDFunction As Int
    Private Azimuth1, Elevation1, SolarOffSet, SolarABOffSet As Double
   
End Sub


#Region Sundial

Sub BtnSunDial_click
    VibBip
    If SDFunction = 0 Then
        SunDialOn
    Else
        SunDialOff
    End If
   
End Sub

Sub SunDialOn
        SDFunction = 1
        IVSun.Visible = True
        IVShadowLine.Visible = True
        IVShadowLine.BringToFront
        IVSun.BringToFront
        SolarBearing1   
       
'UpDate Displays       
        SolarOffSet = Abs(Azimuth1 - AntBearing) Mod 360
        SolarOffSet = Min(360-SolarOffSet, SolarOffSet)
       
        SolarABOffSet = Abs(Azimuth1 - ABearing) Mod 360
        SolarABOffSet = Min(360 - SolarABOffSet, SolarABOffSet)

        lblASolarOffsetNNN.Text = NumberFormat2 (SolarABOffSet, 1 , 1, 1, False)
        lblSolarOffsetNNN.Text = NumberFormat2 (SolarOffSet, 1 , 1, 1, False)
       
        lblSolarAzimuthNNN.Text = NumberFormat2 (Azimuth1, 1 , 1, 1, False)
        lblASolarAzimuthNNN.Text = NumberFormat2 (Azimuth1, 1 , 1, 1, False)
       
        lblSolarElevationNNN.Text = NumberFormat2 (Elevation1, 1 , 1, 1, False)
        lblASolarElevationNNN.Text = NumberFormat2 (Elevation1, 1 , 1, 1, False)
       
        If Mode = 0 Then
            lblSolarOffsetNNN.Visible = True
            lblSolarOffset.Visible = True
            lblSolarAzimuth.Visible = True
            lblSolarAzimuthNNN.Visible = True
            lblSolarElevation.Visible = True
            lblSolarElevationNNN.Visible = True
        Else If Mode = 1 Then
            lblASolarOffsetNNN.Visible = True
            lblASolarOffset.Visible = True
            lblASolarAzimuth.Visible = True
            lblASolarAzimuthNNN.Visible = True
            lblASolarElevation.Visible = True
            lblASolarElevationNNN.Visible = True
        End If   
End Sub

Sub SunDialOff
        SDFunction = 0
        IVSun.Visible = False
        IVShadowLine.Visible = False
       
        lblSolarOffsetNNN.Visible = False
        lblSolarOffset.Visible = False
        lblSolarAzimuthNNN.Visible = False
        lblSolarAzimuth.Visible = False
        lblSolarElevation.Visible = False
        lblSolarElevationNNN.Visible = False
       
        lblASolarOffsetNNN.Visible = False
        lblASolarOffset.Visible = False
        lblASolarAzimuth.Visible = False
        lblASolarAzimuthNNN.Visible = False
        lblASolarElevation.Visible = False
        lblASolarElevationNNN.Visible = False
       
        If Mode = 0 Then cp.Initialize2(MapLat, MapLng, MapZoom, 0, 0 )  'Zero map
        If Mode = 1 Then cp.Initialize2(ABMapLat, ABMapLng, ABMapZoom, 0, 0 )  'Zero map
         gmap.AnimateCamera(cp)
End Sub

Sub SunAndShadow
'Position Sun and Shadow
    If Mode = 0 Then
        IVSun.Top = pnlDispLatLng.Height + 2dip
        IVSun.Left = BTSIV.Left - 2.5dip
        IVShadowLine.Height = BTSIV.Top - pnlDispLatLng.Height - IVSun.Height
    End If
    If Mode = 1 Then
        IVSun.Top = pnlABDispLatLng.Height + 2dip
        IVSun.Left = AEndIV.Left - 2.5dip
        IVShadowLine.Height = AEndIV.Top - pnlABDispLatLng.Height - IVSun.Height
    End If
    IVShadowLine.Top = IVSun.Top + IVSun.Height
    IVShadowLine.Left = IVSun.Left + 11.5dip       
End Sub

Sub SolarBearing1
    Private Mylocation1 As LatLng
    Private twopi As Double = 2*cPI
    Private rad  As Double = cPI/180
    Private dEarthMeanRadius As Double =   6371.01    ' In km
    Private dAstronomicalUnit As Double =   149597890    ' In km
    Private    dElapsedJulianDays, dDecimalHours, dEclipticLongitude, dEclipticObliquity, dSin_EclipticLongitude As Double
    Private dRightAscension, dDeclination, dDeclination, dJulianDate, Longitude1, Latitude1 As Double
    Private dx,dy, dMeanLongitude, dMeanAnomaly, dOmega, dGreenwichMeanSiderealTime, dLocalMeanSiderealTime As Double
    Private dLatitudeInRadians, dHourAngle, dCos_Latitude, dSin_Latitude, dCos_HourAngle, dParallax, dZenithAngle As Double
    Private    liAux1, liAux2, tnow As Long

'Find current location   
    Mylocation1 = gmap.MyLocation
    Latitude1 = Mylocation1.Latitude
    Longitude1 = Mylocation1.Longitude
   
    If Abs(Latitude1) > 24 Then                                            'Inaccurate at latitudes close to the equator
    ''Position Sun and Shadow
        SunAndShadow       
       
    ' 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
'        tnow = DateTime.Now
        tnow = DateTime.Now - DateTime.TimeZoneOffset * DateTime.TicksPerHour
        dDecimalHours = DateTime.GetHour(tnow) + DateTime.GetMinute(tnow)/60 _
                        + DateTime.GetSecond(tnow)/3600

    '    Calculate current Julian Day
        liAux1 = (DateTime.GetMonth(tnow)-14)/12
        liAux2 = (1461*(DateTime.GetYear(tnow) + 4800 + liAux1)) / 4 + (367 * (DateTime.GetMonth(tnow) _
            - 2 -12*liAux1)) / 12- (3 * ((DateTime.GetYear(tnow) + 4900 _
            + liAux1) / 100)) / 4 + DateTime.GetDayOfMonth(tnow) - 32075
       
        dJulianDate = liAux2 - 0.5 + dDecimalHours / 24.0
    '        // Calculate difference between current Julian Day AND JD 2451545.0
        dElapsedJulianDays = dJulianDate - 2451545.0
           
    '    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)   
        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)
       
        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
        dGreenwichMeanSiderealTime = 6.6974243242 +    0.0657098283 * dElapsedJulianDays + dDecimalHours
        dLocalMeanSiderealTime = (dGreenwichMeanSiderealTime * 15 + Longitude1) * rad
        dHourAngle = dLocalMeanSiderealTime - dRightAscension
        dLatitudeInRadians = Latitude1 * 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
        Azimuth1 = ATan2( dy, dx )
        If  Azimuth1< 0.0 Then Azimuth1 = Azimuth1 + twopi
        Azimuth1 = Azimuth1/rad

'        // Parallax Correction
        dParallax = (dEarthMeanRadius/dAstronomicalUnit) * Sin(dZenithAngle)
        dZenithAngle = (dZenithAngle + dParallax)/rad
        Elevation1 = 90 - dZenithAngle

'        Rotate Map
        DegreesSDRoation = Azimuth1
        If Mode = 0 Then
            cp.Initialize2(MapLat, MapLng, MapZoom, 0, 0 )                          'Zero map
            cp.Initialize2(MapLat, MapLng, MapZoom, DegreesSDRoation, 0 )              'Rotate map
        End If
        If Mode = 1 Then
            cp.Initialize2(ABMapLat, ABMapLng, ABMapZoom, 0, 0 )                      'Zero map
            cp.Initialize2(ABMapLat, ABMapLng, ABMapZoom, DegreesSDRoation, 0 )      'Rotate map
        End If
         gmap.AnimateCamera(cp)
    Else
        ToastMessageShow ( "Unavailable at this Latitude.", True)
        SunDialOff
    End If
End Sub

#End Region

Regards Roger