Napkin math climate calculation

climate.py

import numpy as np
from astropy.time import Time
from astropy.coordinates import get_body_barycentric, solar_system_ephemeris
from astropy import units as u

def calculate_earth_sun_distance(date):
    """
    Calculate the distance between Earth and the Sun on a given date.
    
    Parameters:
    date (str): Date in 'YYYY-MM-DD' format
    
    Returns:
    distance (float): Distance in astronomical units (AU)
    """
    # Set the solar system ephemeris to use
    solar_system_ephemeris.set('builtin')
    
    # Convert the input date to Astropy Time object
    time = Time(date)
    
    # Get the barycentric positions of Earth and Sun
    earth_pos = get_body_barycentric('earth', time)
    sun_pos = get_body_barycentric('sun', time)
    
    # Calculate the distance
    distance = np.linalg.norm(earth_pos.xyz - sun_pos.xyz)
    
    return distance.to(u.au).value

# Example usage
date = '2024-07-17'  # You can change this to any date
distance = calculate_earth_sun_distance(date)
print(f"On {date}, the Earth is approximately {distance:.6f} AU from the Sun.")

# Interactive input
user_date = input("Enter a date (YYYY-MM-DD) to calculate Earth-Sun distance: ")
user_distance = calculate_earth_sun_distance(user_date)
print(f"On {user_date}, the Earth is approximately {user_distance:.6f} AU from the Sun.")

# Plotting the Earth-Sun distance over a year
import matplotlib.pyplot as plt
from datetime import datetime, timedelta

start_date = datetime(2024, 1, 1)
dates = [start_date + timedelta(days=i) for i in range(366)]  # 2024 is a leap year
distances = [calculate_earth_sun_distance(date.strftime('%Y-%m-%d')) for date in dates]

plt.figure(figsize=(12, 6))
plt.plot(dates, distances)
plt.title("Earth-Sun Distance Over 2024")
plt.xlabel("Date")
plt.ylabel("Distance (AU)")
plt.grid(True)
plt.show()


#Claude estimate - need to calibrate it to sensor data
def estimate_aerosol_optical_depth(months_since_eruption, eruption_intensity):
    # Constants for a simple decay model
    max_optical_depth = eruption_intensity * 0.1  # Assuming max optical depth is proportional to eruption intensity
    decay_rate = 0.1  # Decay rate per month

    optical_depth = max_optical_depth * math.exp(-decay_rate * months_since_eruption)
    return optical_depth

# Example usage
sunspots = 50  # Number of sunspots
distance = 1  # Distance in Astronomical Units (AU)
eruption_intensity = 6  # On a scale of 1-8 (e.g., VEI scale)
months_since_eruption = 3

aerosol_optical_depth = estimate_aerosol_optical_depth(months_since_eruption, eruption_intensity)
heat = solar_heat_to_earth(sunspots, distance, aerosol_optical_depth)
print(f"Solar heat reaching Earth: {heat:.2f} W/m^2")
print(f"Estimated aerosol optical depth: {aerosol_optical_depth:.4f}")