#!/usr/bin/env python3
"""
Earth Shadow Field Intersection Analyzer
This script analyzes CSV files containing observation metadata to determine
if the field of view of each photographic plate intersects with Earth's shadow
cone in space, where satellites could be eclipsed.
Author: Analysis for eclipse region intersection with observation fields
"""
import os
import csv
import pandas as pd
from datetime import datetime, timezone
import numpy as np
from pathlib import Path
import json
class EarthShadowFieldAnalyzer:
def __init__(self, csv_directory="poss_1/red/"):
self.csv_directory = Path(csv_directory)
self.earth_radius = 6371 # km
self.au = 149597870.7 # km (astronomical unit)
self.results = []
def parse_observation_time(self, date_obs_str):
"""Parse DATE-OBS string into datetime object"""
try:
# Handle format: 1955-04-14T06:32:00
return datetime.fromisoformat(date_obs_str.replace('Z', '+00:00'))
except ValueError:
try:
# Try alternative formats
return datetime.strptime(date_obs_str, "%Y-%m-%dT%H:%M:%S")
except ValueError:
print(f"Warning: Could not parse date string: {date_obs_str}")
return None
def calculate_solar_position(self, dt):
"""Calculate solar position in equatorial coordinates"""
if dt is None:
return None, None
# Julian day calculation
a = (14 - dt.month) // 12
y = dt.year - a
m = dt.month + 12 * a - 3
jd = dt.day + (153 * m + 2) // 5 + 365 * y + y // 4 - y // 100 + y // 400 + 1721119
# Add time of day
jd += (dt.hour + dt.minute/60 + dt.second/3600) / 24
# Calculate solar coordinates (simplified)
n = jd - 2451545.0 # Days since J2000
L = (280.460 + 0.9856474 * n) % 360 # Solar longitude
g = np.radians((357.528 + 0.9856003 * n) % 360) # Solar anomaly
# Solar longitude with equation of center
lambda_sun = np.radians(L + 1.915 * np.sin(g) + 0.020 * np.sin(2*g))
# Obliquity of ecliptic
epsilon = np.radians(23.439 - 0.0000004 * n)
# Convert to equatorial coordinates
solar_ra = np.degrees(np.arctan2(np.cos(epsilon) * np.sin(lambda_sun), np.cos(lambda_sun)))
solar_dec = np.degrees(np.arcsin(np.sin(epsilon) * np.sin(lambda_sun)))
# Ensure RA is in 0-360 range
solar_ra = solar_ra % 360
return solar_ra, solar_dec
def calculate_antisolar_point(self, solar_ra, solar_dec):
"""Calculate the antisolar point (opposite direction from Sun)"""
antisolar_ra = (solar_ra + 180) % 360
antisolar_dec = -solar_dec
return antisolar_ra, antisolar_dec
def angular_separation(self, ra1, dec1, ra2, dec2):
"""Calculate angular separation between two points in degrees"""
ra1, dec1, ra2, dec2 = map(np.radians, [ra1, dec1, ra2, dec2])
cos_sep = (np.sin(dec1) * np.sin(dec2) +
np.cos(dec1) * np.cos(dec2) * np.cos(ra1 - ra2))
# Handle floating point errors
cos_sep = np.clip(cos_sep, -1, 1)
return np.degrees(np.arccos(cos_sep))
def calculate_shadow_intersection(self, obs_datetime, field_ra, field_dec,
field_size_deg=6, altitude_km=40):
"""
Calculate if the observation field intersects Earth's shadow cone
Parameters:
- obs_datetime: observation time
- field_ra, field_dec: field center coordinates (degrees)
- field_size_deg: approximate field of view diameter (degrees)
- altitude_km: satellite altitude for shadow calculation
Returns:
- dict with shadow intersection analysis
"""
if obs_datetime is None:
return {"error": "Invalid datetime"}
# Calculate solar position
solar_ra, solar_dec = self.calculate_solar_position(obs_datetime)
if solar_ra is None:
return {"error": "Could not calculate solar position"}
# Calculate antisolar point (center of Earth's shadow direction)
antisolar_ra, antisolar_dec = self.calculate_antisolar_point(solar_ra, solar_dec)
# Angular separation between field center and antisolar point
field_to_antisolar = self.angular_separation(field_ra, field_dec,
antisolar_ra, antisolar_dec)
# Earth's shadow geometry at satellite altitude
orbit_radius = self.earth_radius + altitude_km
# Shadow cone angle calculation
# Earth's umbra cone half-angle at satellite altitude
shadow_half_angle = np.degrees(np.arctan(self.earth_radius / orbit_radius))
# Penumbra considerations (larger shadow region)
sun_angular_radius = 0.267 # degrees (angular size of Sun from Earth)
penumbra_half_angle = shadow_half_angle + sun_angular_radius
# Field of view considerations
field_radius = field_size_deg / 2
# Intersection analysis
# Distance from field center to shadow center
shadow_center_distance = field_to_antisolar
# Focus only on umbra intersection
umbra_intersects = shadow_center_distance < (shadow_half_angle + field_radius)
field_entirely_in_umbra = shadow_center_distance + field_radius < shadow_half_angle
# Calculate precise umbra area coverage using circle intersection geometry
if field_entirely_in_umbra:
# Entire field is within umbra
umbra_area_fraction = 1.0
elif umbra_intersects:
# Partial intersection - calculate overlapping area of two circles
# Circle 1: field of view (radius = field_radius)
# Circle 2: umbra projection (radius = shadow_half_angle)
# Distance between centers: shadow_center_distance
r1 = field_radius # field radius
r2 = shadow_half_angle # umbra radius
d = shadow_center_distance # distance between centers
if d >= r1 + r2:
# No intersection
umbra_area_fraction = 0.0
elif d <= abs(r2 - r1):
# One circle entirely within the other
if r1 <= r2:
# Field entirely within umbra
umbra_area_fraction = 1.0
else:
# Umbra entirely within field
umbra_area_fraction = (r2 / r1) ** 2
else:
# Partial intersection - use circle intersection formula
# Area of intersection of two circles
alpha1 = 2 * np.arccos((d**2 + r1**2 - r2**2) / (2 * d * r1))
alpha2 = 2 * np.arccos((d**2 + r2**2 - r1**2) / (2 * d * r2))
area1 = 0.5 * r1**2 * (alpha1 - np.sin(alpha1)) # segment of circle 1
area2 = 0.5 * r2**2 * (alpha2 - np.sin(alpha2)) # segment of circle 2
intersection_area = area1 + area2
field_area = np.pi * r1**2
umbra_area_fraction = intersection_area / field_area
# Ensure fraction is within bounds
umbra_area_fraction = np.clip(umbra_area_fraction, 0.0, 1.0)
else:
# No intersection
umbra_area_fraction = 0.0
return {
"observation_datetime": obs_datetime.isoformat(),
"field_center_ra": field_ra,
"field_center_dec": field_dec,
"solar_ra": round(solar_ra, 3),
"solar_dec": round(solar_dec, 3),
"antisolar_ra": round(antisolar_ra, 3),
"antisolar_dec": round(antisolar_dec, 3),
"field_to_antisolar_deg": round(field_to_antisolar, 3),
"umbra_half_angle_deg": round(shadow_half_angle, 3),
"field_radius_deg": round(field_radius, 3),
"umbra_intersects": umbra_intersects,
"field_entirely_in_umbra": field_entirely_in_umbra,
"umbra_area_fraction": round(umbra_area_fraction, 4),
"umbra_area_percentage": round(umbra_area_fraction * 100, 2),
"satellites_could_be_eclipsed": umbra_intersects,
"assumed_altitude_km": altitude_km
}
def process_csv_file(self, csv_path):
"""Process a single CSV file and extract observation parameters"""
try:
# Read CSV file
df = pd.read_csv(csv_path)
# Extract key parameters
params = {}
for _, row in df.iterrows():
params[row['Key']] = row['Value']
# Get required parameters
date_obs = params.get('DATE-OBS')
plate_ra = float(params.get('PLATERA', 0))
plate_dec = float(params.get('PLATEDEC', 0))
region = params.get('REGION', 'Unknown')
plate_id = params.get('PLATEID', 'Unknown')
# Calculate field size from plate dimensions and scale
plate_scale = float(params.get('PLTSCALE', 67.2)) # arcsec/mm
plate_size_x = float(params.get('PLTSIZEX', 355)) # mm
plate_size_y = float(params.get('PLTSIZEY', 355)) # mm
# Field size in degrees (use diagonal)
field_size_arcsec = plate_scale * np.sqrt(plate_size_x**2 + plate_size_y**2)
field_size_deg = field_size_arcsec / 3600
# Parse observation time
obs_datetime = self.parse_observation_time(date_obs)
# Calculate shadow intersection
shadow_analysis = self.calculate_shadow_intersection(
obs_datetime, plate_ra, plate_dec, field_size_deg
)
# Combine all information
result = {
"csv_file": str(csv_path),
"region": region,
"plate_id": plate_id,
"date_obs": date_obs,
"calculated_field_size_deg": round(field_size_deg, 3),
**shadow_analysis
}
return result
except Exception as e:
print(f"Error processing {csv_path}: {e}")
return {
"csv_file": str(csv_path),
"error": str(e)
}
def process_all_csv_files(self):
"""Process all CSV files in the directory"""
if not self.csv_directory.exists():
print(f"Directory {self.csv_directory} does not exist!")
return
csv_files = list(self.csv_directory.glob("*.csv"))
print(f"Found {len(csv_files)} CSV files to process...")
for csv_file in csv_files:
print(f"Processing: {csv_file.name}")
result = self.process_csv_file(csv_file)
self.results.append(result)
# Calculate and display running average (only for valid results)
valid_results_so_far = [r for r in self.results if 'error' not in r]
if valid_results_so_far:
current_coverages = [r.get('umbra_area_fraction', 0) for r in valid_results_so_far]
running_avg = np.mean(current_coverages)
print(f" → Running average umbra coverage: {running_avg:.4f} ({running_avg*100:.2f}%)")
def save_results(self, output_file="shadow_intersection_analysis.json"):
"""Save results to JSON file"""
with open(output_file, 'w') as f:
json.dump(self.results, f, indent=2, default=str)
print(f"Results saved to {output_file}")
def generate_summary_report(self):
"""Generate a summary report of the shadow intersection analysis"""
if not self.results:
print("No results to summarize!")
return
total_observations = len(self.results)
valid_analyses = [r for r in self.results if 'error' not in r]
if not valid_analyses:
print("No valid analyses found!")
return
umbra_intersections = sum(1 for r in valid_analyses if r.get('umbra_intersects', False))
satellites_eclipsed_possible = sum(1 for r in valid_analyses if r.get('satellites_could_be_eclipsed', False))
field_entirely_in_umbra = sum(1 for r in valid_analyses if r.get('field_entirely_in_umbra', False))
significant_coverage = sum(1 for r in valid_analyses if r.get('umbra_area_fraction', 0) > 0.1)
# Calculate average coverage for fields that have any intersection
intersecting_fields = [r for r in valid_analyses if r.get('umbra_intersects', False)]
avg_coverage_intersecting = np.mean([r.get('umbra_area_fraction', 0) for r in intersecting_fields]) if intersecting_fields else 0
# Calculate overall average coverage across ALL plates (including zero coverage)
all_coverages = [r.get('umbra_area_fraction', 0) for r in valid_analyses]
overall_avg_coverage = np.mean(all_coverages) if all_coverages else 0
print("\n" + "="*70)
print("EARTH UMBRA FIELD INTERSECTION ANALYSIS")
print("="*70)
print(f"Total CSV files processed: {total_observations}")
print(f"Valid analyses: {len(valid_analyses)}")
print(f"Fields intersecting umbra: {umbra_intersections}")
print(f"Fields where satellites could be eclipsed: {satellites_eclipsed_possible}")
print(f"Fields entirely within umbra: {field_entirely_in_umbra}")
print(f"Fields with >10% umbra coverage: {significant_coverage}")
print(f"\nPercentages:")
print(f" {(satellites_eclipsed_possible/len(valid_analyses)*100):.1f}% of fields could contain eclipsed satellites")
print(f" {(umbra_intersections/len(valid_analyses)*100):.1f}% intersect Earth's umbra")
print(f" {(field_entirely_in_umbra/len(valid_analyses)*100):.1f}% are entirely within umbra")
print(f"\nUMBRA COVERAGE STATISTICS:")
print(f" Overall average coverage (all plates): {overall_avg_coverage:.1%} ({overall_avg_coverage*100:.2f}%)")
if intersecting_fields:
print(f" Average coverage (intersecting plates only): {avg_coverage_intersecting:.1%}")
print(f" Maximum coverage: {max(r.get('umbra_area_fraction', 0) for r in intersecting_fields):.1%}")
print(f" Total plates with umbra intersection: {len(intersecting_fields)}")
# Show best examples
print(f"\nFields with highest umbra coverage:")
shadow_fields = [r for r in valid_analyses if r.get('satellites_could_be_eclipsed', False)]
shadow_fields.sort(key=lambda x: x.get('umbra_area_fraction', 0), reverse=True)
for example in shadow_fields[:5]:
print(f" {example['plate_id']} ({example['region']}) - {example['date_obs']}")
print(f" Umbra coverage: {example['umbra_area_percentage']:.2f}% of plate area")
print(f" Distance to antisolar point: {example['field_to_antisolar_deg']:.1f}°")
def create_summary_csv(self, output_file="shadow_intersection_summary.csv"):
"""Create a CSV summary of all analyses"""
if not self.results:
print("No results to export!")
return
# Filter out error results and create DataFrame
valid_results = [r for r in self.results if 'error' not in r]
if valid_results:
df = pd.DataFrame(valid_results)
df.to_csv(output_file, index=False)
print(f"Summary CSV saved to {output_file}")
def main():
# Initialize analyzer
analyzer = EarthShadowFieldAnalyzer("poss_1/red/")
# Process all CSV files
analyzer.process_all_csv_files()
# Save detailed results
analyzer.save_results("shadow_intersection_analysis.json")
# Create summary CSV
analyzer.create_summary_csv("shadow_intersection_summary.csv")
# Generate summary report
analyzer.generate_summary_report()
print("\nAnalysis complete!")
print("Files created:")
print(" - shadow_intersection_analysis.json (detailed results)")
print(" - shadow_intersection_summary.csv (summary table)")
print("\nLook for fields with 'satellites_could_be_eclipsed' = True")
print("Focus on 'umbra_area_percentage' for coverage of plate area")
if __name__ == "__main__":
main()
). In fact, can you write out a complete list of your analysis plan? Doesn't have to be super detailed.
