Add geometric solar position method to get_solarposition

docs/sphinx/source/reference/solarposition.rst

@@ -17,6 +17,7 @@ unless you know that you need a different function.
    solarposition.get_solarposition
    solarposition.spa_python
    solarposition.ephemeris
+   solarposition.geometric
    solarposition.pyephem
    solarposition.spa_c
 

docs/sphinx/source/whatsnew/v0.16.2.rst

@@ -27,6 +27,10 @@ Enhancements
 * Add ``dataset`` parameter to :py:func:`~pvlib.iotools.get_era5`
   allowing for choosing between ERA5 and ERA5-Land datasets.
   (:pull:`2780`)
+* Add :py:func:`pvlib.solarposition.geometric`, a fast but approximate solar
+  position calculation based on spherical geometry, and make it available as
+  ``method='geometric'`` in :py:func:`pvlib.solarposition.get_solarposition`.
+  (:issue:`2914`, :pull:`2915`)
 
 Documentation
 ~~~~~~~~~~~~~
@@ -69,3 +73,4 @@ Contributors
 * jason-writer (:ghuser:`jason-writer`)
 * :ghuser:`bfs2021`
 * :ghuser:`xiajiahang`
+* Divyanshi Singh (:ghuser:`divyanshisingh987456321`)

pvlib/solarposition.py

@@ -69,6 +69,9 @@ def get_solarposition(time, latitude, longitude,
 
         'nrel_c' uses the NLR SPA C code [3]: :py:func:`spa_c`
 
+        'geometric' uses a fast but approximate spherical-geometry
+        calculation: :py:func:`geometric`
+
     temperature : float, default 12
         Degrees C.
 
@@ -119,6 +122,8 @@ def get_solarposition(time, latitude, longitude,
     elif method == 'ephemeris':
         ephem_df = ephemeris(time, latitude, longitude, pressure, temperature,
                              **kwargs)
+    elif method == 'geometric':
+        ephem_df = geometric(time, latitude, longitude, **kwargs)
     else:
         raise ValueError('Invalid solar position method')
 
@@ -878,6 +883,83 @@ def ephemeris(time, latitude, longitude, pressure=101325.0, temperature=12.0):
     return result
 
 
+def geometric(time, latitude, longitude):
+    """
+    Fast, approximate solar position from spherical geometry.
+
+    The solar declination and the equation of time are evaluated with the
+    Spencer (1971) Fourier series, and the zenith and azimuth angles follow
+    from spherical trigonometry. This is much faster than, but less
+    accurate than, :py:func:`spa_python`; use it only when speed matters more
+    than accuracy.
+
+    .. warning::
+        Atmospheric refraction is neglected, so ``apparent_zenith`` and
+        ``apparent_elevation`` are equal to ``zenith`` and ``elevation``.
+
+    Parameters
+    ----------
+    time : pandas.DatetimeIndex
+        Must be localized or UTC will be assumed.
+    latitude : float
+        Latitude in decimal degrees. Positive north of equator, negative
+        to south.
+    longitude : float
+        Longitude in decimal degrees. Positive east of prime meridian,
+        negative to west.
+
+    Returns
+    -------
+    DataFrame with the following columns:
+
+        * apparent_elevation : equal to ``elevation`` (no refraction).
+        * elevation : elevation of the sun in decimal degrees, 0 = on
+          horizon. The complement of the zenith angle.
+        * azimuth : Azimuth of the sun in decimal degrees East of North.
+        * apparent_zenith : equal to ``zenith`` (no refraction).
+        * zenith : Solar zenith angle in decimal degrees.
+
+    See Also
+    --------
+    spa_python, ephemeris, declination_spencer71,
+    equation_of_time_spencer71, solar_zenith_analytical,
+    solar_azimuth_analytical
+    """
+    time_utc = tools._pandas_to_utc(time)  # naive times assumed to be UTC
+
+    dayofyear = time_utc.dayofyear
+    declination = declination_spencer71(dayofyear)
+    eot = equation_of_time_spencer71(dayofyear)
+    utc_hours = (time_utc.hour + time_utc.minute / 60. +
+                 time_utc.second / 3600. + time_utc.microsecond / 3.6e9)
+    # hour angle in degrees, as in hour_angle() but computed from UTC time
+    # (no timezone offset to apply)
+    hourangle = np.radians(15. * (np.asarray(utc_hours) - 12.) + longitude +
+                           eot / 4.)
+    # wrap to [-pi, pi): the sign of the hour angle selects the azimuth branch
+    hourangle = np.mod(hourangle + np.pi, 2 * np.pi) - np.pi
+    latitude_rad = np.radians(latitude)
+
+    zenith = solar_zenith_analytical(latitude_rad, hourangle, declination)
+    azimuth = solar_azimuth_analytical(
+        latitude_rad, hourangle, declination, zenith)
+
+    zenith = np.degrees(zenith)
+    azimuth = np.mod(np.degrees(azimuth), 360)
+    elevation = 90 - zenith
+
+    return pd.DataFrame(
+        {
+            "apparent_elevation": elevation,
+            "elevation": elevation,
+            "azimuth": azimuth,
+            "apparent_zenith": zenith,
+            "zenith": zenith,
+        },
+        index=time,
+    )
+
+
 def calc_time(lower_bound, upper_bound, latitude, longitude, attribute, value,
               altitude=0.0, pressure=101325.0, temperature=12.0,
               horizon='+0:00', xtol=1.0e-12):

tests/test_solarposition.py

@@ -463,6 +463,58 @@ def test_ephemeris_physical_no_tz(expected_solpos, golden_mst):
     assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])
 
 
+@pytest.mark.parametrize('latitude, longitude, tz', [
+    (39.742476, -105.1786, 'America/Denver'),  # west of Greenwich
+    (28.46, 77.03, 'Asia/Kolkata'),  # east of Greenwich
+    (-33.87, 151.21, 'Australia/Sydney'),  # southern hemisphere
+    (69.65, 18.96, 'Europe/Oslo'),  # high latitude
+])
+def test_geometric_vs_spa(latitude, longitude, tz):
+    # the geometric method is approximate: compare with SPA across a whole
+    # year, including evenings at western longitudes where the hour angle
+    # leaves [-180, 180] degrees before wrapping
+    times = pd.date_range('2026-01-01', '2026-12-31 23:00', freq='1h', tz=tz)
+    expected = solarposition.spa_python(times, latitude, longitude)
+    result = solarposition.geometric(times, latitude, longitude)
+    assert_allclose(result['zenith'], expected['zenith'], atol=1)
+    assert_allclose(result['elevation'], expected['elevation'], atol=1)
+    daytime = expected['elevation'] > 0
+    azimuth_error = (result['azimuth'] - expected['azimuth'] + 180) % 360 - 180
+    assert (azimuth_error[daytime].abs() < 2).all()
+
+
+def test_geometric_output(golden):
+    result = solarposition.geometric(times_localized, golden.latitude,
+                                     golden.longitude)
+    assert result.index.equals(times_localized)
+    assert list(result.columns) == ['apparent_elevation', 'elevation',
+                                    'azimuth', 'apparent_zenith', 'zenith']
+    assert not result.isna().any().any()
+    assert result['azimuth'].between(0, 360).all()
+    assert_allclose(result['zenith'] + result['elevation'], 90)
+    # refraction is neglected
+    assert_series_equal(result['apparent_zenith'], result['zenith'],
+                        check_names=False)
+
+
+def test_geometric_no_tz():
+    # times without a timezone are assumed to be UTC
+    times_naive = pd.date_range('2026-06-21 00:00', periods=24, freq='1h')
+    expected = solarposition.geometric(times_naive.tz_localize('UTC'),
+                                       32.2, -111)
+    result = solarposition.geometric(times_naive, 32.2, -111)
+    assert_allclose(result.to_numpy(), expected.to_numpy())
+
+
+def test_get_solarposition_method_geometric(golden):
+    expected = solarposition.geometric(times_localized, golden.latitude,
+                                       golden.longitude)
+    result = solarposition.get_solarposition(
+        times_localized, golden.latitude, golden.longitude,
+        method='geometric')
+    assert_frame_equal(result, expected)
+
+
 def test_get_solarposition_error(golden):
     times = pd.date_range(datetime.datetime(2003, 10, 17, 13, 30, 30),
                           periods=1, freq='D', tz=golden.tz)