Skip to content

Commit 43379ea

Browse files
twpoltwpol
committed
style: Format sun/moon position code
1 parent 654dfd1 commit 43379ea

File tree

1 file changed

+114
-100
lines changed

1 file changed

+114
-100
lines changed

‎Source/RunActivity/Viewer3D/Common/SunMoonPos.cs‎

Lines changed: 114 additions & 100 deletions
Original file line numberDiff line numberDiff line change
@@ -1,32 +1,29 @@
1-
// COPYRIGHT 2010 by the Open Rails project.
2-
//
1+
// COPYRIGHT 2009 - 2023 by the Open Rails project.
2+
//
33
// This file is part of Open Rails.
4-
//
4+
//
55
// Open Rails is free software: you can redistribute it and/or modify
66
// it under the terms of the GNU General Public License as published by
77
// the Free Software Foundation, either version 3 of the License, or
88
// (at your option) any later version.
9-
//
9+
//
1010
// Open Rails is distributed in the hope that it will be useful,
1111
// but WITHOUT ANY WARRANTY; without even the implied warranty of
1212
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1313
// GNU General Public License for more details.
14-
//
14+
//
1515
// You should have received a copy of the GNU General Public License
1616
// along with Open Rails. If not, see <http://www.gnu.org/licenses/>.
1717

1818
/*
1919
* Contains equations to calculate the sun and moon position for any latitude
2020
* and longitude on any date
21-
* Solar postion equations are NOAA equations adapted from "Astronomical Algorithms," by Jean Meeus,
21+
* Solar position equations are NOAA equations adapted from "Astronomical Algorithms," by Jean Meeus,
2222
* Lunar equations are adapted by Keith Burnett from the US Naval Observatory Astronomical Almanac.
2323
*/
24-
// Principal Author:
25-
// Rick Grout
26-
//
2724

28-
using Microsoft.Xna.Framework;
2925
using System;
26+
using Microsoft.Xna.Framework;
3027

3128
namespace Orts.Viewer3D.Common
3229
{
@@ -36,67 +33,67 @@ static class SunMoonPos
3633
/// Calculates the solar direction vector.
3734
/// Used for locating the sun graphic and as the location of the main scenery light source.
3835
/// </summary>
39-
/// <param name="latitude">latitude</param>
40-
/// <param name="longitude">longitude</param>
41-
/// <param name="clockTime">wall clock time since start of activity, days</param>
42-
/// <param name="date">structure made up of day, month, year and ordinal date</param>
36+
/// <param name="latitude">Latitude, radians.</param>
37+
/// <param name="longitude">Longitude, radians.</param>
38+
/// <param name="clockTime">Wall clock time since start of activity, days.</param>
39+
/// <param name="date">Structure made up of day, month, year and ordinal date.</param>
40+
/// <returns>A normalize 3D vector indicating the sun's position from the viewer's location.</returns>
4341
public static Vector3 SolarAngle(double latitude, double longitude, float clockTime, SkyViewer.SkyDate date)
4442
{
4543
Vector3 sunDirection;
4644

47-
// For these calculations, west longitude is in positive degrees,
48-
float NOAAlongitude = -MathHelper.ToDegrees((float)longitude);
45+
// For these calculations, west longitude is in positive degrees
46+
var noaaLongitude = -MathHelper.ToDegrees((float)longitude);
47+
4948
// Fractional year, radians
50-
double fYear = (MathHelper.TwoPi / 365) * (date.OrdinalDate - 1 + (clockTime - 0.5));
49+
var fYear = (MathHelper.TwoPi / 365) * (date.OrdinalDate - 1 + (clockTime - 0.5));
50+
5151
// Equation of time, minutes
52-
double eqTime = 229.18 * (0.000075
53-
+ 0.001868 * Math.Cos(fYear)
54-
- 0.032077 * Math.Sin(fYear)
55-
- 0.014615 * Math.Cos(2 * fYear)
56-
- 0.040849 * Math.Sin(2 * fYear));
52+
var eqTime = 229.18 * (0.000075
53+
+ (0.001868 * Math.Cos(fYear))
54+
- (0.032077 * Math.Sin(fYear))
55+
- (0.014615 * Math.Cos(2 * fYear))
56+
- (0.040849 * Math.Sin(2 * fYear)));
57+
5758
// Solar declination, radians
58-
double solarDeclination = 0.006918
59-
- 0.399912 * Math.Cos(fYear)
60-
+ 0.070257 * Math.Sin(fYear)
61-
- 0.006758 * Math.Cos(2 * fYear)
62-
+ 0.000907 * Math.Sin(2 * fYear)
63-
- 0.002697 * Math.Cos(3 * fYear)
64-
+ 0.001480 * Math.Sin(3 * fYear);
59+
var solarDeclination = 0.006918
60+
- (0.399912 * Math.Cos(fYear))
61+
+ (0.070257 * Math.Sin(fYear))
62+
- (0.006758 * Math.Cos(2 * fYear))
63+
+ (0.000907 * Math.Sin(2 * fYear))
64+
- (0.002697 * Math.Cos(3 * fYear))
65+
+ (0.001480 * Math.Sin(3 * fYear));
66+
6567
// Time offset at present longitude, minutes
66-
double timeOffset = eqTime - 4 * NOAAlongitude + 60 * Math.Round(NOAAlongitude / 15);
68+
var timeOffset = eqTime - (4 * noaaLongitude) + (60 * Math.Round(noaaLongitude / 15));
69+
6770
// True solar time, minutes (since midnight)
68-
double trueSolar = clockTime * 24 * 60 + timeOffset;
71+
var trueSolar = (clockTime * 24 * 60) + timeOffset;
72+
6973
// Solar hour angle, radians
70-
double solarHourAngle = MathHelper.ToRadians((float)(trueSolar / 4) - 180);
74+
var solarHourAngle = MathHelper.ToRadians((float)(trueSolar / 4) - 180);
7175

7276
// Solar zenith cosine. This is the Y COORDINATE of the solar Vector.
73-
double solarZenithCosine = Math.Sin(latitude)
74-
* Math.Sin(solarDeclination)
75-
+ Math.Cos(latitude)
76-
* Math.Cos(solarDeclination)
77-
* Math.Cos(solarHourAngle);
77+
var solarZenithCosine = (Math.Sin(latitude) * Math.Sin(solarDeclination))
78+
+ (Math.Cos(latitude) * Math.Cos(solarDeclination) * Math.Cos(solarHourAngle));
7879

7980
// Solar elevation angle, radians. Currently not used.
8081
// double solarElevationAngle = MathHelper.PiOver2 - Math.Acos(solarZenithCosine);
8182

8283
// Solar azimuth cosine. This is the Z COORDINATE of the solar Vector.
83-
double solarAzimuthCosine = -(Math.Sin(latitude)
84-
* solarZenithCosine
85-
- Math.Sin(solarDeclination)) / (
86-
+Math.Cos(latitude)
87-
* Math.Sin(Math.Acos(solarZenithCosine)));
84+
var solarAzimuthCosine = -((Math.Sin(latitude) * solarZenithCosine) - Math.Sin(solarDeclination))
85+
/ (Math.Cos(latitude) * Math.Sin(Math.Acos(solarZenithCosine)));
8886

8987
// Running at 64 bit solarAzimuthCosine can be slightly below -1, generating NaN results
90-
if (solarAzimuthCosine > 1.0d) solarAzimuthCosine = 1.0d;
91-
if (solarAzimuthCosine < -1.0d) solarAzimuthCosine = -1.0d;
88+
solarAzimuthCosine = MathHelper.Clamp((float)solarAzimuthCosine, -1, 1);
9289

9390
// Solar azimuth angle, radians. Currently not used.
9491
// double solarAzimuthAngle = Math.Acos(solarAzimuthCosine);
9592
// if (clockTime > 0.5)
9693
// solarAzimuthAngle = MathHelper.TwoPi - solarAzimuthAngle;
9794

9895
// Solar azimuth sine. This is the X COORDINATE of the solar Vector.
99-
double solarAzimuthSine = Math.Sin(Math.Acos(solarAzimuthCosine)) * (clockTime > 0.5 ? 1 : -1);
96+
var solarAzimuthSine = Math.Sin(Math.Acos(solarAzimuthCosine)) * (clockTime > 0.5 ? 1 : -1);
10097

10198
sunDirection.X = -(float)solarAzimuthSine;
10299
sunDirection.Y = (float)solarZenithCosine;
@@ -106,83 +103,102 @@ public static Vector3 SolarAngle(double latitude, double longitude, float clockT
106103
}
107104

108105
/// <summary>
109-
/// Calculates the lunar direction vector.
106+
/// Calculates the lunar direction vector.
110107
/// </summary>
111-
/// <param name="latitude">latitude</param>
112-
/// <param name="longitude">longitude</param>
113-
/// <param name="clockTime">wall clock time since start of activity</param>
114-
/// <param name="date">structure made up of day, month, year and ordinal date</param>
108+
/// <param name="latitude">Latitude, radians.</param>
109+
/// <param name="longitude">Longitude, radians.</param>
110+
/// <param name="clockTime">Wall clock time since start of activity, days.</param>
111+
/// <param name="date">Structure made up of day, month, year and ordinal date.</param>
112+
/// <returns>A normalize 3D vector indicating the moon's position from the viewer's location.</returns>
115113
public static Vector3 LunarAngle(double latitude, double longitude, float clockTime, SkyViewer.SkyDate date)
116114
{
117115
Vector3 moonDirection;
118116

119117
// Julian day.
120-
double JEday = (float)367 * date.Year - 7 * (date.Year + (date.Month + 9) / 12) / 4 + 275 * date.Month / 9 + date.Day - 730531.5 + clockTime;
118+
var jEday = (367f * date.Year) - (7 * (date.Year + ((date.Month + 9) / 12)) / 4) + (275 * date.Month / 9) + date.Day - 730531.5 + clockTime;
119+
121120
// Fractional time within current 100-year epoch, days
122-
double Ftime = JEday / 36525;
121+
var fTime = jEday / 36525;
122+
123123
// Ecliptic longitude, radians
124-
double GeoLong = Normalize(Normalize(3.8104 + 8399.71 * Ftime, MathHelper.TwoPi)
125-
+ 0.10978 * Math.Sin(Normalize(2.3544 + 8328.69 * Ftime, MathHelper.TwoPi))
126-
- 0.02216 * Math.Sin(Normalize(4.5239 - 7214.12 * Ftime, MathHelper.TwoPi))
127-
+ 0.0115 * Math.Sin(Normalize(4.11374 + 15542.75 * Ftime, MathHelper.TwoPi))
128-
+ 0.003665 * Math.Sin(Normalize(4.7106 + 16657.38 * Ftime, MathHelper.TwoPi))
129-
- 0.003316 * Math.Sin(Normalize(6.2396 + 628.3 * Ftime, MathHelper.TwoPi))
130-
- 0.00192 * Math.Sin(Normalize(3.2568 + 16866.93 * Ftime, MathHelper.TwoPi)), MathHelper.TwoPi);
124+
var geoLong = Normalize(
125+
Normalize(3.8104 + (8399.71 * fTime), MathHelper.TwoPi)
126+
+ (0.10978 * Math.Sin(Normalize(2.3544 + (8328.69 * fTime), MathHelper.TwoPi)))
127+
- (0.02216 * Math.Sin(Normalize(4.5239 - (7214.12 * fTime), MathHelper.TwoPi)))
128+
+ (0.0115 * Math.Sin(Normalize(4.11374 + (15542.75 * fTime), MathHelper.TwoPi)))
129+
+ (0.003665 * Math.Sin(Normalize(4.7106 + (16657.38 * fTime), MathHelper.TwoPi)))
130+
- (0.003316 * Math.Sin(Normalize(6.2396 + (628.3 * fTime), MathHelper.TwoPi)))
131+
- (0.00192 * Math.Sin(Normalize(3.2568 + (16866.93 * fTime), MathHelper.TwoPi))), MathHelper.TwoPi);
132+
131133
// Ecliptic latitude, radians
132-
double GeoLat = 0.08954 * Math.Sin(Normalize(1.6284 + 8433.47 * Ftime, MathHelper.TwoPi))
133-
+ 0.004887 * Math.Sin(Normalize(3.9828 + 16762.16 * Ftime, MathHelper.TwoPi))
134-
- 0.004887 * Math.Sin(Normalize(5.5554 + 104.76 * Ftime, MathHelper.TwoPi))
135-
- 0.002967 * Math.Sin(Normalize(3.7978 - 7109.29 * Ftime, MathHelper.TwoPi));
134+
var geoLat = (0.08954 * Math.Sin(Normalize(1.6284 + (8433.47 * fTime), MathHelper.TwoPi)))
135+
+ (0.004887 * Math.Sin(Normalize(3.9828 + (16762.16 * fTime), MathHelper.TwoPi)))
136+
- (0.004887 * Math.Sin(Normalize(5.5554 + (104.76 * fTime), MathHelper.TwoPi)))
137+
- (0.002967 * Math.Sin(Normalize(3.7978 - (7109.29 * fTime), MathHelper.TwoPi)));
138+
136139
// Parallax, radians
137-
double Parallax = 0.0165946
138-
+ 0.0009041 * Math.Cos(Normalize(2.3544 + 832869 * Ftime, MathHelper.TwoPi))
139-
+ 0.000166 * Math.Cos(Normalize(4.5238 - 7214.06 * Ftime, MathHelper.TwoPi))
140-
+ 0.000136 * Math.Cos(Normalize(4.1137 + 15542.75 * Ftime, MathHelper.TwoPi))
141-
+ 0.000489 * Math.Cos(Normalize(4.7106 + 16657.38 * Ftime, MathHelper.TwoPi));
140+
var parallax = 0.0165946
141+
+ (0.0009041 * Math.Cos(Normalize(2.3544 + (832869 * fTime), MathHelper.TwoPi)))
142+
+ (0.000166 * Math.Cos(Normalize(4.5238 - (7214.06 * fTime), MathHelper.TwoPi)))
143+
+ (0.000136 * Math.Cos(Normalize(4.1137 + (15542.75 * fTime), MathHelper.TwoPi)))
144+
+ (0.000489 * Math.Cos(Normalize(4.7106 + (16657.38 * fTime), MathHelper.TwoPi)));
145+
142146
// Geocentric distance, dimensionless
143-
double GeoDist = 1 / Math.Sin(Parallax);
147+
var geoDist = 1 / Math.Sin(parallax);
148+
144149
// Geocentric vector coordinates, dimensionless
145-
double geoX = GeoDist * Math.Cos(GeoLong) * Math.Cos(GeoLat);
146-
double geoY = GeoDist * Math.Sin(GeoLong) * Math.Cos(GeoLat);
147-
double geoZ = GeoDist * Math.Sin(GeoLat);
150+
var geoX = geoDist * Math.Cos(geoLong) * Math.Cos(geoLat);
151+
var geoY = geoDist * Math.Sin(geoLong) * Math.Cos(geoLat);
152+
var geoZ = geoDist * Math.Sin(geoLat);
153+
148154
// Mean obliquity of ecliptic, radians
149-
double Ecliptic = (0.4091 - 0.00000000623 * JEday);
155+
var ecliptic = 0.4091 - (0.00000000623 * jEday);
156+
150157
// Equator of ordinalDate vector coordinates, dimensionless
151-
double eclX = geoX;
152-
double eclY = geoY * Math.Cos(Ecliptic) - geoZ * Math.Sin(Ecliptic);
153-
double eclZ = geoY * Math.Sin(Ecliptic) + geoZ * Math.Cos(Ecliptic);
158+
var eclX = geoX;
159+
var eclY = (geoY * Math.Cos(ecliptic)) - (geoZ * Math.Sin(ecliptic));
160+
var eclZ = (geoY * Math.Sin(ecliptic)) + (geoZ * Math.Cos(ecliptic));
161+
154162
// Right ascension and declination, radians
155-
double RightAsc = Math.Atan(eclY / eclX);
163+
var rightAsc = Math.Atan(eclY / eclX);
156164
if (eclX < 0)
157-
RightAsc += MathHelper.Pi;
165+
{
166+
rightAsc += MathHelper.Pi;
167+
}
168+
158169
if ((eclY < 0) && (eclX > 0))
159-
RightAsc += MathHelper.TwoPi;
160-
double Declination = Math.Atan(eclZ / (Math.Pow((Math.Pow(eclX, 2) + Math.Pow(eclY, 2)), 0.5)));
170+
{
171+
rightAsc += MathHelper.TwoPi;
172+
}
173+
174+
var declination = Math.Atan(eclZ / Math.Pow(Math.Pow(eclX, 2) + Math.Pow(eclY, 2), 0.5));
175+
161176
// Convert right ascension and declination to altitude and azimuth.
162177
// Equations by Stephen R. Schmitt.
163178
// Greenwich Mean Sidereal Time, degrees
164-
double GMST = 280.46061837 + 360.98564736629 * JEday;
165-
GMST = Normalize(GMST, 360);
179+
var gMST = 280.46061837 + (360.98564736629 * jEday);
180+
gMST = Normalize(gMST, 360);
181+
166182
// Local Mean Sidereal Time, degrees
167-
double LMST = GMST + MathHelper.ToDegrees((float)longitude) / 15;
183+
var lMST = gMST + (MathHelper.ToDegrees((float)longitude) / 15);
184+
168185
// Hour angle, degrees
169-
double HA = LMST - MathHelper.ToDegrees((float)RightAsc) / 15;
186+
var hA = lMST - (MathHelper.ToDegrees((float)rightAsc) / 15);
187+
170188
// Hour angle, radians
171-
double lunarHourAngle = MathHelper.ToRadians((float)HA);
189+
var lunarHourAngle = MathHelper.ToRadians((float)hA);
190+
172191
// Lunar Altitude, radians from its Sine
173-
double lunarAltSin = Math.Sin(latitude) * Math.Sin(Declination) + Math.Cos(latitude) * Math.Cos(lunarHourAngle);
174-
double lunarAltitude = Math.Asin(lunarAltSin);
192+
var lunarAltSin = (Math.Sin(latitude) * Math.Sin(declination)) + (Math.Cos(latitude) * Math.Cos(lunarHourAngle));
193+
var lunarAltitude = Math.Asin(lunarAltSin);
194+
175195
// Lunar Azimuth, radians from its Cosine
176-
double lunarAltCos = Math.Cos(lunarAltitude);
177-
double hourAngleSin = Math.Sin(lunarHourAngle);
178-
double lunarAzCos = (Math.Sin(Declination) - lunarAltSin * Math.Sin(latitude)) / (lunarAltCos * Math.Cos(latitude));
196+
var lunarAltCos = Math.Cos(lunarAltitude);
197+
var hourAngleSin = Math.Sin(lunarHourAngle);
198+
var lunarAzCos = (Math.Sin(declination) - (lunarAltSin * Math.Sin(latitude))) / (lunarAltCos * Math.Cos(latitude));
179199
lunarAzCos = lunarAzCos < -1 ? -1 : lunarAzCos;
180200
lunarAzCos = lunarAzCos > 1 ? 1 : lunarAzCos;
181-
double lunarAzimuth;
182-
if (hourAngleSin < 0)
183-
lunarAzimuth = Math.Acos(lunarAzCos);
184-
else
185-
lunarAzimuth = MathHelper.TwoPi - Math.Acos(lunarAzCos);
201+
var lunarAzimuth = hourAngleSin < 0 ? Math.Acos(lunarAzCos) : MathHelper.TwoPi - Math.Acos(lunarAzCos);
186202

187203
moonDirection.X = (float)Math.Sin(lunarAzimuth);
188204
moonDirection.Y = (float)lunarAltSin;
@@ -194,13 +210,11 @@ public static Vector3 LunarAngle(double latitude, double longitude, float clockT
194210
/// <summary>
195211
/// Removes all multiples of "divisor" from the input number.
196212
/// </summary>
197-
/// <param name="input">the raw number</param>
198-
/// <param name="divisor">the number, or its multiples, we want to remove</param>
213+
/// <param name="input">the raw number.</param>
214+
/// <param name="divisor">the number, or its multiples, we want to remove.</param>
199215
static double Normalize(double input, double divisor)
200216
{
201-
double output = input - divisor * Math.Floor(input / divisor);
202-
return output;
217+
return input - (divisor * Math.Floor(input / divisor));
203218
}
204-
205219
}
206220
}

0 commit comments

Comments
 (0)