# Android QuestionCalculating Solar Azimuth [Solved]

#### Roger Daley

##### Well-Known Member
Longtime User
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:
B4X:
``````        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

#### Roger Daley

##### Well-Known Member
Longtime User
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?

B4X:
``````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``````

B4X:
``````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

##### Expert
Longtime User
The sub returns an array

B4X:
``AzimuthAstro = Astro2.SunPosition(...)``

should be
B4X:
``Private AzimuthAstro() As Double``

i guess

#### derez

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

#### Roger Daley

##### Well-Known Member
Longtime User
Hi Guys.

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

B4X:
``````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

##### Expert
Longtime User
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

##### Well-Known Member
Longtime User
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

##### Expert
Longtime User
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

##### Well-Known Member
Longtime User
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

#### Attachments

• ABTWIP.zip
126 KB · Views: 289

#### Roger Daley

##### Well-Known Member
Longtime User
To close out the question this is the code I finally used. Many thanks to Derez.

B4X:
``````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
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

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

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
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

' 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
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

'        // Parallax Correction
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

Replies
0
Views
249
Replies
5
Views
2K
Replies
51
Views
9K
Replies
14
Views
10K
Replies
36
Views
10K