We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Nightshade w/ refraction=-6.0 yields ValueError: GEOSGeom_createLinearRing_r returned a NULL pointer
### Description

I'm trying use Nightshade to draw twilight lines (and sub solar point).

If I do Nightshade w/ refraction=-6.0, I get ValueError: GEOSGeom_createLinearRing_r returned a NULL pointer

Nightshade w/ refraction=-6.0001 works OK

Also, refraction=89 causes extraneous lines.

#### Code to reproduce

```
import matplotlib.pyplot as plt
fig = plt.figure()

import cartopy
ax = fig.add_subplot(projection=cartopy.crs.PlateCarree())
ax.coastlines()

import datetime
date = datetime.datetime(1999, 12, 31, 12)

from cartopy.feature.nightshade import Nightshade

ax.add_feature(Nightshade(date, color='none'))
ax.add_feature(Nightshade(date, refraction=-6.0, color='none')) # Civil Twilight
#ax.add_feature(Nightshade(date, refraction=-6.0001, color='none')) # Civil Twilight
ax.add_feature(Nightshade(date, refraction=-12, color='none')) # Nautical Twilight

ax.add_feature(Nightshade(date, refraction=89, color='none')) # sub solar circle ...

plt.show()
```

#### Traceback 

```
python /tmp/nightshade-bug.py 
/usr/lib64/python3.10/site-packages/cartopy/feature/nightshade.py:80: RuntimeWarning: invalid value encountered in arccos
  omega0 = np.rad2deg(np.arccos(np.sin(np.deg2rad(refraction)) /
IllegalArgumentException: Points of LinearRing do not form a closed linestring
Traceback (most recent call last):
  File "/tmp/nightshade-bug.py", line 15, in <module>
    ax.add_feature(Nightshade(date, refraction=-6, color='none')) # Civil Twilight
  File "/usr/lib64/python3.10/site-packages/cartopy/feature/nightshade.py", line 93, in __init__
    geom = sgeom.Polygon(np.column_stack((x, y)))
  File "/usr/lib64/python3.10/site-packages/shapely/geometry/polygon.py", line 243, in __init__
    ret = geos_polygon_from_py(shell, holes)
  File "/usr/lib64/python3.10/site-packages/shapely/geometry/polygon.py", line 509, in geos_polygon_from_py
    ret = geos_linearring_from_py(shell)
  File "shapely/speedups/_speedups.pyx", line 408, in shapely.speedups._speedups.geos_linearring_from_py
ValueError: GEOSGeom_createLinearRing_r returned a NULL pointer
```

<details>
 <summary>Full environment definition</summary>

<!-- fill in the following information as appropriate -->

### Operating system

Fedora 35

### Cartopy version

python3-cartopy-0.20.1-1.fc35.x86_64

python3-shapely-1.7.1-12.fc35.x86_64
python3-matplotlib-3.5.1-3.fc35.x86_64

</details>

Thanks !

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `lib/cartopy/feature/nightshade.py`
Content:
```
1 # Copyright Cartopy Contributors
2 #
3 # This file is part of Cartopy and is released under the LGPL license.
4 # See COPYING and COPYING.LESSER in the root of the repository for full
5 # licensing details.
6 
7 import datetime
8 
9 import numpy as np
10 import shapely.geometry as sgeom
11 
12 from . import ShapelyFeature
13 from .. import crs as ccrs
14 
15 
16 class Nightshade(ShapelyFeature):
17     def __init__(self, date=None, delta=0.1, refraction=-0.83,
18                  color="k", alpha=0.5, **kwargs):
19         """
20         Shade the darkside of the Earth, accounting for refraction.
21 
22         Parameters
23         ----------
24         date : datetime
25             A UTC datetime object used to calculate the position of the sun.
26             Default: datetime.datetime.utcnow()
27         delta : float
28             Stepsize in degrees to determine the resolution of the
29             night polygon feature (``npts = 180 / delta``).
30         refraction : float
31             The adjustment in degrees due to refraction,
32             thickness of the solar disc, elevation etc...
33 
34         Note
35         ----
36             Matplotlib keyword arguments can be used when drawing the feature.
37             This allows standard Matplotlib control over aspects such as
38             'color', 'alpha', etc.
39 
40         """
41         if date is None:
42             date = datetime.datetime.utcnow()
43 
44         # make sure date is UTC, or naive with respect to time zones
45         if date.utcoffset():
46             raise ValueError(
47                 f'datetime instance must be UTC, not {date.tzname()}')
48 
49         # Returns the Greenwich hour angle,
50         # need longitude (opposite direction)
51         lat, lon = _solar_position(date)
52         pole_lon = lon
53         if lat > 0:
54             pole_lat = -90 + lat
55             central_lon = 180
56         else:
57             pole_lat = 90 + lat
58             central_lon = 0
59 
60         rotated_pole = ccrs.RotatedPole(pole_latitude=pole_lat,
61                                         pole_longitude=pole_lon,
62                                         central_rotated_longitude=central_lon)
63 
64         npts = int(180/delta)
65         x = np.empty(npts*2)
66         y = np.empty(npts*2)
67 
68         # Solve the equation for sunrise/sunset:
69         # https://en.wikipedia.org/wiki/Sunrise_equation#Generalized_equation
70         # NOTE: In the generalized equation on Wikipedia,
71         #       delta == 0. in the rotated pole coordinate system.
72         #       Therefore, the max/min latitude is +/- (90+refraction)
73 
74         # Fill latitudes up and then down
75         y[:npts] = np.linspace(-(90+refraction), 90+refraction, npts)
76         y[npts:] = y[:npts][::-1]
77 
78         # Solve the generalized equation for omega0, which is the
79         # angle of sunrise/sunset from solar noon
80         omega0 = np.rad2deg(np.arccos(np.sin(np.deg2rad(refraction)) /
81                                       np.cos(np.deg2rad(y))))
82 
83         # Fill the longitude values from the offset for midnight.
84         # This needs to be a closed loop to fill the polygon.
85         # Negative longitudes
86         x[:npts] = -(180 - omega0[:npts])
87         # Positive longitudes
88         x[npts:] = 180 - omega0[npts:]
89 
90         kwargs.setdefault('facecolor', color)
91         kwargs.setdefault('alpha', alpha)
92 
93         geom = sgeom.Polygon(np.column_stack((x, y)))
94         return super().__init__(
95             [geom], rotated_pole, **kwargs)
96 
97 
98 def _julian_day(date):
99     """
100     Calculate the Julian day from an input datetime.
101 
102     Parameters
103     ----------
104     date
105         A UTC datetime object.
106 
107     Note
108     ----
109     Algorithm implemented following equations from Chapter 3 (Algorithm 14):
110     Vallado, David 'Fundamentals of Astrodynamics and Applications', (2007)
111 
112     Julian day epoch is: noon on January 1, 4713 BC (proleptic Julian)
113                          noon on November 24, 4714 BC (proleptic Gregorian)
114 
115     """
116     year = date.year
117     month = date.month
118     day = date.day
119     hour = date.hour
120     minute = date.minute
121     second = date.second
122 
123     # January/February correspond to months 13/14 respectively
124     # for the constants to work out properly
125     if month < 3:
126         month += 12
127         year -= 1
128 
129     B = 2 - year // 100 + (year // 100) // 4
130     C = ((second/60 + minute)/60 + hour)/24
131 
132     JD = (int(365.25*(year + 4716)) + int(30.6001*(month+1)) +
133           day + B - 1524.5 + C)
134     return JD
135 
136 
137 def _solar_position(date):
138     """
139     Calculate the latitude and longitude point where the sun is
140     directly overhead for the given date.
141 
142     Parameters
143     ----------
144     date
145         A UTC datetime object.
146 
147     Returns
148     -------
149     (latitude, longitude) in degrees
150 
151     Note
152     ----
153     Algorithm implemented following equations from Chapter 5 (Algorithm 29):
154     Vallado, David 'Fundamentals of Astrodynamics and Applications', (2007)
155 
156     """
157     # NOTE: Constants are in degrees in the textbook,
158     #       so we need to convert the values from deg2rad when taking sin/cos
159 
160     # Centuries from J2000
161     T_UT1 = (_julian_day(date) - 2451545.0)/36525
162 
163     # solar longitude (deg)
164     lambda_M_sun = (280.460 + 36000.771*T_UT1) % 360
165 
166     # solar anomaly (deg)
167     M_sun = (357.5277233 + 35999.05034*T_UT1) % 360
168 
169     # ecliptic longitude
170     lambda_ecliptic = (lambda_M_sun + 1.914666471*np.sin(np.deg2rad(M_sun)) +
171                        0.019994643*np.sin(np.deg2rad(2*M_sun)))
172 
173     # obliquity of the ecliptic (epsilon in Vallado's notation)
174     epsilon = 23.439291 - 0.0130042*T_UT1
175 
176     # declination of the sun
177     delta_sun = np.rad2deg(np.arcsin(np.sin(np.deg2rad(epsilon)) *
178                                      np.sin(np.deg2rad(lambda_ecliptic))))
179 
180     # Greenwich mean sidereal time (seconds)
181     theta_GMST = (67310.54841 +
182                   (876600*3600 + 8640184.812866)*T_UT1 +
183                   0.093104*T_UT1**2 -
184                   6.2e-6*T_UT1**3)
185     # Convert to degrees
186     theta_GMST = (theta_GMST % 86400)/240
187 
188     # Right ascension calculations
189     numerator = (np.cos(np.deg2rad(epsilon)) *
190                  np.sin(np.deg2rad(lambda_ecliptic)) /
191                  np.cos(np.deg2rad(delta_sun)))
192     denominator = (np.cos(np.deg2rad(lambda_ecliptic)) /
193                    np.cos(np.deg2rad(delta_sun)))
194 
195     alpha_sun = np.rad2deg(np.arctan2(numerator, denominator))
196 
197     # longitude is opposite of Greenwich Hour Angle (GHA)
198     # GHA == theta_GMST - alpha_sun
199     lon = -(theta_GMST-alpha_sun)
200     if lon < -180:
201         lon += 360
202 
203     return (delta_sun, lon)
204
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/lib/cartopy/feature/nightshade.py b/lib/cartopy/feature/nightshade.py
--- a/lib/cartopy/feature/nightshade.py
+++ b/lib/cartopy/feature/nightshade.py
@@ -77,8 +77,12 @@
 
         # Solve the generalized equation for omega0, which is the
         # angle of sunrise/sunset from solar noon
-        omega0 = np.rad2deg(np.arccos(np.sin(np.deg2rad(refraction)) /
-                                      np.cos(np.deg2rad(y))))
+        # We need to clip the input to arccos to [-1, 1] due to floating
+        # point precision and arccos creating nans for values outside
+        # of the domain
+        arccos_tmp = np.clip(np.sin(np.deg2rad(refraction)) /
+                             np.cos(np.deg2rad(y)), -1, 1)
+        omega0 = np.rad2deg(np.arccos(arccos_tmp))
 
         # Fill the longitude values from the offset for midnight.
         # This needs to be a closed loop to fill the polygon.

{"golden_diff": "diff --git a/lib/cartopy/feature/nightshade.py b/lib/cartopy/feature/nightshade.py\n--- a/lib/cartopy/feature/nightshade.py\n+++ b/lib/cartopy/feature/nightshade.py\n@@ -77,8 +77,12 @@\n \n         # Solve the generalized equation for omega0, which is the\n         # angle of sunrise/sunset from solar noon\n-        omega0 = np.rad2deg(np.arccos(np.sin(np.deg2rad(refraction)) /\n-                                      np.cos(np.deg2rad(y))))\n+        # We need to clip the input to arccos to [-1, 1] due to floating\n+        # point precision and arccos creating nans for values outside\n+        # of the domain\n+        arccos_tmp = np.clip(np.sin(np.deg2rad(refraction)) /\n+                             np.cos(np.deg2rad(y)), -1, 1)\n+        omega0 = np.rad2deg(np.arccos(arccos_tmp))\n \n         # Fill the longitude values from the offset for midnight.\n         # This needs to be a closed loop to fill the polygon.\n", "issue": "Nightshade w/ refraction=-6.0 yields ValueError: GEOSGeom_createLinearRing_r returned a NULL pointer\n### Description\r\n\r\nI'm trying use Nightshade to draw twilight lines (and sub solar point).\r\n\r\nIf I do Nightshade w/ refraction=-6.0, I get ValueError: GEOSGeom_createLinearRing_r returned a NULL pointer\r\n\r\nNightshade w/ refraction=-6.0001 works OK\r\n\r\nAlso, refraction=89 causes extraneous lines.\r\n\r\n#### Code to reproduce\r\n\r\n```\r\nimport matplotlib.pyplot as plt\r\nfig = plt.figure()\r\n\r\nimport cartopy\r\nax = fig.add_subplot(projection=cartopy.crs.PlateCarree())\r\nax.coastlines()\r\n\r\nimport datetime\r\ndate = datetime.datetime(1999, 12, 31, 12)\r\n\r\nfrom cartopy.feature.nightshade import Nightshade\r\n\r\nax.add_feature(Nightshade(date, color='none'))\r\nax.add_feature(Nightshade(date, refraction=-6.0, color='none')) # Civil Twilight\r\n#ax.add_feature(Nightshade(date, refraction=-6.0001, color='none')) # Civil Twilight\r\nax.add_feature(Nightshade(date, refraction=-12, color='none')) # Nautical Twilight\r\n\r\nax.add_feature(Nightshade(date, refraction=89, color='none')) # sub solar circle ...\r\n\r\nplt.show()\r\n```\r\n\r\n#### Traceback \r\n\r\n```\r\npython /tmp/nightshade-bug.py \r\n/usr/lib64/python3.10/site-packages/cartopy/feature/nightshade.py:80: RuntimeWarning: invalid value encountered in arccos\r\n  omega0 = np.rad2deg(np.arccos(np.sin(np.deg2rad(refraction)) /\r\nIllegalArgumentException: Points of LinearRing do not form a closed linestring\r\nTraceback (most recent call last):\r\n  File \"/tmp/nightshade-bug.py\", line 15, in <module>\r\n    ax.add_feature(Nightshade(date, refraction=-6, color='none')) # Civil Twilight\r\n  File \"/usr/lib64/python3.10/site-packages/cartopy/feature/nightshade.py\", line 93, in __init__\r\n    geom = sgeom.Polygon(np.column_stack((x, y)))\r\n  File \"/usr/lib64/python3.10/site-packages/shapely/geometry/polygon.py\", line 243, in __init__\r\n    ret = geos_polygon_from_py(shell, holes)\r\n  File \"/usr/lib64/python3.10/site-packages/shapely/geometry/polygon.py\", line 509, in geos_polygon_from_py\r\n    ret = geos_linearring_from_py(shell)\r\n  File \"shapely/speedups/_speedups.pyx\", line 408, in shapely.speedups._speedups.geos_linearring_from_py\r\nValueError: GEOSGeom_createLinearRing_r returned a NULL pointer\r\n```\r\n\r\n<details>\r\n <summary>Full environment definition</summary>\r\n\r\n<!-- fill in the following information as appropriate -->\r\n\r\n### Operating system\r\n\r\nFedora 35\r\n\r\n### Cartopy version\r\n\r\npython3-cartopy-0.20.1-1.fc35.x86_64\r\n\r\npython3-shapely-1.7.1-12.fc35.x86_64\r\npython3-matplotlib-3.5.1-3.fc35.x86_64\r\n\r\n</details>\r\n\r\nThanks !\n", "before_files": [{"content": "# Copyright Cartopy Contributors\n#\n# This file is part of Cartopy and is released under the LGPL license.\n# See COPYING and COPYING.LESSER in the root of the repository for full\n# licensing details.\n\nimport datetime\n\nimport numpy as np\nimport shapely.geometry as sgeom\n\nfrom . import ShapelyFeature\nfrom .. import crs as ccrs\n\n\nclass Nightshade(ShapelyFeature):\n    def __init__(self, date=None, delta=0.1, refraction=-0.83,\n                 color=\"k\", alpha=0.5, **kwargs):\n        \"\"\"\n        Shade the darkside of the Earth, accounting for refraction.\n\n        Parameters\n        ----------\n        date : datetime\n            A UTC datetime object used to calculate the position of the sun.\n            Default: datetime.datetime.utcnow()\n        delta : float\n            Stepsize in degrees to determine the resolution of the\n            night polygon feature (``npts = 180 / delta``).\n        refraction : float\n            The adjustment in degrees due to refraction,\n            thickness of the solar disc, elevation etc...\n\n        Note\n        ----\n            Matplotlib keyword arguments can be used when drawing the feature.\n            This allows standard Matplotlib control over aspects such as\n            'color', 'alpha', etc.\n\n        \"\"\"\n        if date is None:\n            date = datetime.datetime.utcnow()\n\n        # make sure date is UTC, or naive with respect to time zones\n        if date.utcoffset():\n            raise ValueError(\n                f'datetime instance must be UTC, not {date.tzname()}')\n\n        # Returns the Greenwich hour angle,\n        # need longitude (opposite direction)\n        lat, lon = _solar_position(date)\n        pole_lon = lon\n        if lat > 0:\n            pole_lat = -90 + lat\n            central_lon = 180\n        else:\n            pole_lat = 90 + lat\n            central_lon = 0\n\n        rotated_pole = ccrs.RotatedPole(pole_latitude=pole_lat,\n                                        pole_longitude=pole_lon,\n                                        central_rotated_longitude=central_lon)\n\n        npts = int(180/delta)\n        x = np.empty(npts*2)\n        y = np.empty(npts*2)\n\n        # Solve the equation for sunrise/sunset:\n        # https://en.wikipedia.org/wiki/Sunrise_equation#Generalized_equation\n        # NOTE: In the generalized equation on Wikipedia,\n        #       delta == 0. in the rotated pole coordinate system.\n        #       Therefore, the max/min latitude is +/- (90+refraction)\n\n        # Fill latitudes up and then down\n        y[:npts] = np.linspace(-(90+refraction), 90+refraction, npts)\n        y[npts:] = y[:npts][::-1]\n\n        # Solve the generalized equation for omega0, which is the\n        # angle of sunrise/sunset from solar noon\n        omega0 = np.rad2deg(np.arccos(np.sin(np.deg2rad(refraction)) /\n                                      np.cos(np.deg2rad(y))))\n\n        # Fill the longitude values from the offset for midnight.\n        # This needs to be a closed loop to fill the polygon.\n        # Negative longitudes\n        x[:npts] = -(180 - omega0[:npts])\n        # Positive longitudes\n        x[npts:] = 180 - omega0[npts:]\n\n        kwargs.setdefault('facecolor', color)\n        kwargs.setdefault('alpha', alpha)\n\n        geom = sgeom.Polygon(np.column_stack((x, y)))\n        return super().__init__(\n            [geom], rotated_pole, **kwargs)\n\n\ndef _julian_day(date):\n    \"\"\"\n    Calculate the Julian day from an input datetime.\n\n    Parameters\n    ----------\n    date\n        A UTC datetime object.\n\n    Note\n    ----\n    Algorithm implemented following equations from Chapter 3 (Algorithm 14):\n    Vallado, David 'Fundamentals of Astrodynamics and Applications', (2007)\n\n    Julian day epoch is: noon on January 1, 4713 BC (proleptic Julian)\n                         noon on November 24, 4714 BC (proleptic Gregorian)\n\n    \"\"\"\n    year = date.year\n    month = date.month\n    day = date.day\n    hour = date.hour\n    minute = date.minute\n    second = date.second\n\n    # January/February correspond to months 13/14 respectively\n    # for the constants to work out properly\n    if month < 3:\n        month += 12\n        year -= 1\n\n    B = 2 - year // 100 + (year // 100) // 4\n    C = ((second/60 + minute)/60 + hour)/24\n\n    JD = (int(365.25*(year + 4716)) + int(30.6001*(month+1)) +\n          day + B - 1524.5 + C)\n    return JD\n\n\ndef _solar_position(date):\n    \"\"\"\n    Calculate the latitude and longitude point where the sun is\n    directly overhead for the given date.\n\n    Parameters\n    ----------\n    date\n        A UTC datetime object.\n\n    Returns\n    -------\n    (latitude, longitude) in degrees\n\n    Note\n    ----\n    Algorithm implemented following equations from Chapter 5 (Algorithm 29):\n    Vallado, David 'Fundamentals of Astrodynamics and Applications', (2007)\n\n    \"\"\"\n    # NOTE: Constants are in degrees in the textbook,\n    #       so we need to convert the values from deg2rad when taking sin/cos\n\n    # Centuries from J2000\n    T_UT1 = (_julian_day(date) - 2451545.0)/36525\n\n    # solar longitude (deg)\n    lambda_M_sun = (280.460 + 36000.771*T_UT1) % 360\n\n    # solar anomaly (deg)\n    M_sun = (357.5277233 + 35999.05034*T_UT1) % 360\n\n    # ecliptic longitude\n    lambda_ecliptic = (lambda_M_sun + 1.914666471*np.sin(np.deg2rad(M_sun)) +\n                       0.019994643*np.sin(np.deg2rad(2*M_sun)))\n\n    # obliquity of the ecliptic (epsilon in Vallado's notation)\n    epsilon = 23.439291 - 0.0130042*T_UT1\n\n    # declination of the sun\n    delta_sun = np.rad2deg(np.arcsin(np.sin(np.deg2rad(epsilon)) *\n                                     np.sin(np.deg2rad(lambda_ecliptic))))\n\n    # Greenwich mean sidereal time (seconds)\n    theta_GMST = (67310.54841 +\n                  (876600*3600 + 8640184.812866)*T_UT1 +\n                  0.093104*T_UT1**2 -\n                  6.2e-6*T_UT1**3)\n    # Convert to degrees\n    theta_GMST = (theta_GMST % 86400)/240\n\n    # Right ascension calculations\n    numerator = (np.cos(np.deg2rad(epsilon)) *\n                 np.sin(np.deg2rad(lambda_ecliptic)) /\n                 np.cos(np.deg2rad(delta_sun)))\n    denominator = (np.cos(np.deg2rad(lambda_ecliptic)) /\n                   np.cos(np.deg2rad(delta_sun)))\n\n    alpha_sun = np.rad2deg(np.arctan2(numerator, denominator))\n\n    # longitude is opposite of Greenwich Hour Angle (GHA)\n    # GHA == theta_GMST - alpha_sun\n    lon = -(theta_GMST-alpha_sun)\n    if lon < -180:\n        lon += 360\n\n    return (delta_sun, lon)\n", "path": "lib/cartopy/feature/nightshade.py"}], "after_files": [{"content": "# Copyright Cartopy Contributors\n#\n# This file is part of Cartopy and is released under the LGPL license.\n# See COPYING and COPYING.LESSER in the root of the repository for full\n# licensing details.\n\nimport datetime\n\nimport numpy as np\nimport shapely.geometry as sgeom\n\nfrom . import ShapelyFeature\nfrom .. import crs as ccrs\n\n\nclass Nightshade(ShapelyFeature):\n    def __init__(self, date=None, delta=0.1, refraction=-0.83,\n                 color=\"k\", alpha=0.5, **kwargs):\n        \"\"\"\n        Shade the darkside of the Earth, accounting for refraction.\n\n        Parameters\n        ----------\n        date : datetime\n            A UTC datetime object used to calculate the position of the sun.\n            Default: datetime.datetime.utcnow()\n        delta : float\n            Stepsize in degrees to determine the resolution of the\n            night polygon feature (``npts = 180 / delta``).\n        refraction : float\n            The adjustment in degrees due to refraction,\n            thickness of the solar disc, elevation etc...\n\n        Note\n        ----\n            Matplotlib keyword arguments can be used when drawing the feature.\n            This allows standard Matplotlib control over aspects such as\n            'color', 'alpha', etc.\n\n        \"\"\"\n        if date is None:\n            date = datetime.datetime.utcnow()\n\n        # make sure date is UTC, or naive with respect to time zones\n        if date.utcoffset():\n            raise ValueError(\n                f'datetime instance must be UTC, not {date.tzname()}')\n\n        # Returns the Greenwich hour angle,\n        # need longitude (opposite direction)\n        lat, lon = _solar_position(date)\n        pole_lon = lon\n        if lat > 0:\n            pole_lat = -90 + lat\n            central_lon = 180\n        else:\n            pole_lat = 90 + lat\n            central_lon = 0\n\n        rotated_pole = ccrs.RotatedPole(pole_latitude=pole_lat,\n                                        pole_longitude=pole_lon,\n                                        central_rotated_longitude=central_lon)\n\n        npts = int(180/delta)\n        x = np.empty(npts*2)\n        y = np.empty(npts*2)\n\n        # Solve the equation for sunrise/sunset:\n        # https://en.wikipedia.org/wiki/Sunrise_equation#Generalized_equation\n        # NOTE: In the generalized equation on Wikipedia,\n        #       delta == 0. in the rotated pole coordinate system.\n        #       Therefore, the max/min latitude is +/- (90+refraction)\n\n        # Fill latitudes up and then down\n        y[:npts] = np.linspace(-(90+refraction), 90+refraction, npts)\n        y[npts:] = y[:npts][::-1]\n\n        # Solve the generalized equation for omega0, which is the\n        # angle of sunrise/sunset from solar noon\n        # We need to clip the input to arccos to [-1, 1] due to floating\n        # point precision and arccos creating nans for values outside\n        # of the domain\n        arccos_tmp = np.clip(np.sin(np.deg2rad(refraction)) /\n                             np.cos(np.deg2rad(y)), -1, 1)\n        omega0 = np.rad2deg(np.arccos(arccos_tmp))\n\n        # Fill the longitude values from the offset for midnight.\n        # This needs to be a closed loop to fill the polygon.\n        # Negative longitudes\n        x[:npts] = -(180 - omega0[:npts])\n        # Positive longitudes\n        x[npts:] = 180 - omega0[npts:]\n\n        kwargs.setdefault('facecolor', color)\n        kwargs.setdefault('alpha', alpha)\n\n        geom = sgeom.Polygon(np.column_stack((x, y)))\n        return super().__init__(\n            [geom], rotated_pole, **kwargs)\n\n\ndef _julian_day(date):\n    \"\"\"\n    Calculate the Julian day from an input datetime.\n\n    Parameters\n    ----------\n    date\n        A UTC datetime object.\n\n    Note\n    ----\n    Algorithm implemented following equations from Chapter 3 (Algorithm 14):\n    Vallado, David 'Fundamentals of Astrodynamics and Applications', (2007)\n\n    Julian day epoch is: noon on January 1, 4713 BC (proleptic Julian)\n                         noon on November 24, 4714 BC (proleptic Gregorian)\n\n    \"\"\"\n    year = date.year\n    month = date.month\n    day = date.day\n    hour = date.hour\n    minute = date.minute\n    second = date.second\n\n    # January/February correspond to months 13/14 respectively\n    # for the constants to work out properly\n    if month < 3:\n        month += 12\n        year -= 1\n\n    B = 2 - year // 100 + (year // 100) // 4\n    C = ((second/60 + minute)/60 + hour)/24\n\n    JD = (int(365.25*(year + 4716)) + int(30.6001*(month+1)) +\n          day + B - 1524.5 + C)\n    return JD\n\n\ndef _solar_position(date):\n    \"\"\"\n    Calculate the latitude and longitude point where the sun is\n    directly overhead for the given date.\n\n    Parameters\n    ----------\n    date\n        A UTC datetime object.\n\n    Returns\n    -------\n    (latitude, longitude) in degrees\n\n    Note\n    ----\n    Algorithm implemented following equations from Chapter 5 (Algorithm 29):\n    Vallado, David 'Fundamentals of Astrodynamics and Applications', (2007)\n\n    \"\"\"\n    # NOTE: Constants are in degrees in the textbook,\n    #       so we need to convert the values from deg2rad when taking sin/cos\n\n    # Centuries from J2000\n    T_UT1 = (_julian_day(date) - 2451545.0)/36525\n\n    # solar longitude (deg)\n    lambda_M_sun = (280.460 + 36000.771*T_UT1) % 360\n\n    # solar anomaly (deg)\n    M_sun = (357.5277233 + 35999.05034*T_UT1) % 360\n\n    # ecliptic longitude\n    lambda_ecliptic = (lambda_M_sun + 1.914666471*np.sin(np.deg2rad(M_sun)) +\n                       0.019994643*np.sin(np.deg2rad(2*M_sun)))\n\n    # obliquity of the ecliptic (epsilon in Vallado's notation)\n    epsilon = 23.439291 - 0.0130042*T_UT1\n\n    # declination of the sun\n    delta_sun = np.rad2deg(np.arcsin(np.sin(np.deg2rad(epsilon)) *\n                                     np.sin(np.deg2rad(lambda_ecliptic))))\n\n    # Greenwich mean sidereal time (seconds)\n    theta_GMST = (67310.54841 +\n                  (876600*3600 + 8640184.812866)*T_UT1 +\n                  0.093104*T_UT1**2 -\n                  6.2e-6*T_UT1**3)\n    # Convert to degrees\n    theta_GMST = (theta_GMST % 86400)/240\n\n    # Right ascension calculations\n    numerator = (np.cos(np.deg2rad(epsilon)) *\n                 np.sin(np.deg2rad(lambda_ecliptic)) /\n                 np.cos(np.deg2rad(delta_sun)))\n    denominator = (np.cos(np.deg2rad(lambda_ecliptic)) /\n                   np.cos(np.deg2rad(delta_sun)))\n\n    alpha_sun = np.rad2deg(np.arctan2(numerator, denominator))\n\n    # longitude is opposite of Greenwich Hour Angle (GHA)\n    # GHA == theta_GMST - alpha_sun\n    lon = -(theta_GMST-alpha_sun)\n    if lon < -180:\n        lon += 360\n\n    return (delta_sun, lon)\n", "path": "lib/cartopy/feature/nightshade.py"}]}