gistfile1.txt

import pvlib
import pandas as pd
import numpy as np
from scipy.optimize import minimize
import pytz
from pytz import timezone
from sunposition import observed_sunpos
from skyfield.api import Topos, load
import numpy as np
import folium


def find_latitudes_for_solar_angle(utc_time, solar_angle, longitude=0, tolerance=0.1):
    """
    Find possible latitudes for a given solar angle at a specific UTC time.

    Parameters:
        utc_time: tuple - UTC time (year, month, day, hour, minute, second)
        solar_angle: float - Desired solar elevation angle in degrees
        longitude: float - Longitude to check (default is 0 for Greenwich)
        tolerance: float - Tolerance for matching the solar angle

    Returns:
        list of latitudes where the solar angle matches
    """
    # Load timescale and ephemeris
    ts = load.timescale()
    planets = load("de421.bsp")
    earth = planets["earth"]
    sun = planets["sun"]

    # Define UTC time
    time = ts.utc(*utc_time)

    # Iterate through latitudes
    min_south_lat, max_south_lat = -90, 0
    min_north_lat, max_north_lat = 90, 0

    for lat in np.linspace(-90, 90, 1000):  # Adjust step size for precision
        # Define the observer's location
        observer = earth + Topos(latitude_degrees=lat, longitude_degrees=longitude)

        # Observe the Sun from the observer's location
        astrometric = observer.at(time).observe(sun).apparent()

        # Get the altitude (elevation) of the Sun
        altitude, _, _ = astrometric.altaz(temperature_C=24.0, pressure_mbar=1005)

        # Check if the altitude matches the desired solar angle within the tolerance
        if abs(altitude.degrees - solar_angle) <= tolerance:
            lat = float(lat)
            if lat < 0:
                if lat > min_south_lat:
                    min_south_lat = lat
                if lat < max_south_lat:
                    max_south_lat = lat
            if lat > 0:
                if lat < min_north_lat:
                    min_north_lat = lat
                if lat > max_north_lat:
                    max_north_lat = lat

    return {
        "min_south_lat": min_south_lat,
        "max_south_lat": max_south_lat,
        "min_north_lat": min_north_lat,
        "max_north_lat": max_north_lat,
    }


# REAL REFERENCE  DATA
# https://www.suncalc.org/#/-7.1457,-34.8044,19/2016.09.02/09:38/1/1

# Example usage
utc_time = (2016, 9, 2, 12, 38, 0)  # YYYY, MM, DD, HH, MM, SS in UTC
solar_angle = 60.87  # Solar elevation angle in degrees
lat_long_dict = {}
precision_degrees = 1
for long in range(-180, 180, precision_degrees):
    possible_latitudes = find_latitudes_for_solar_angle(utc_time, solar_angle, long)
    if (
        possible_latitudes["min_south_lat"] != -90
        or possible_latitudes["max_south_lat"] != 0
        or possible_latitudes["min_north_lat"] != 90
        or possible_latitudes["max_north_lat"] != 0
    ):
        print(
            f"Possible latitudes for long {long} at solar angle {solar_angle}° at {utc_time}: {possible_latitudes}"
        )
        lat_long_dict[long] = possible_latitudes


# Function to create polygons for the latitude range
def create_range_polygons(lat_long_dict):
    southern_polygon = []
    northern_polygon = []

    for long, data in sorted(lat_long_dict.items()):
        # Add points for southern range (min and max latitudes at this longitude)
        if data["min_south_lat"] != -90 and data["max_south_lat"] != 0:
            southern_polygon.append([data["min_south_lat"], long])
            southern_polygon.insert(
                0, [data["max_south_lat"], long]
            )  # Insert to close ring

        # Add points for northern range (min and max latitudes at this longitude)
        if data["min_north_lat"] != 90 and data["max_north_lat"] != 0:
            northern_polygon.append([data["min_north_lat"], long])
            northern_polygon.insert(
                0, [data["max_north_lat"], long]
            )  # Insert to close ring

    return southern_polygon, northern_polygon


# Function to add polygons to the map
def add_range_to_map(map_object, polygon, color, popup_text):
    if len(polygon) > 2:  # Ensure there are enough points to create a polygon
        folium.Polygon(
            locations=polygon,
            color=color,
            fill=True,
            fill_opacity=0.5,
            popup=folium.Popup(popup_text),
        ).add_to(map_object)


# Main function to plot latitude ranges as rings
def plot_latitude_range_rings(lat_long_dict):
    # Create a base map
    map_center = [0, 0]
    folium_map = folium.Map(location=map_center, zoom_start=2)

    # Generate polygons for southern and northern ranges
    southern_polygon, northern_polygon = create_range_polygons(lat_long_dict)

    # Add the southern and northern polygons to the map
    add_range_to_map(
        folium_map,
        southern_polygon,
        color="red",
        popup_text="Southern Hemisphere Range",
    )
    add_range_to_map(
        folium_map,
        northern_polygon,
        color="blue",
        popup_text="Northern Hemisphere Range",
    )

    # Save the map to an HTML file
    folium_map.save("latitude_range_ring_map.html")
    print(
        "Map saved as 'latitude_range_ring_map.html'. Open this file in a browser to view it."
    )


# Plot the latitude ranges as rings
plot_latitude_range_rings(lat_long_dict)