import os
import glob
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.time import Time
from earthshadow import get_shadow_center, get_shadow_radius, dist_from_shadow_center
from astropy.coordinates import SkyCoord, SkyOffsetFrame
import astropy.units as u
from tqdm import tqdm
from scipy import stats
def time_str(date_obs, ut):
date_split = date_obs.split("/")
return f"19{date_split[2]}-{date_split[1]}-{date_split[0]} {ut}"
# See https://github.com/guynir42/earthshadow
def in_shadow( c, time):
center = get_shadow_center(time, obs='Palomar', orbit='GEO')
radius = get_shadow_radius(orbit='GEO')
dist = dist_from_shadow_center( c.ra.degree, c.dec.degree, time=time, obs='Palomar', orbit='GEO')
return dist < radius - 2 * u.deg
def randomize_time_within_exposure( t, exposure ):
return t + np.random.uniform( 0.0, exposure ) * u.minute
def randomize_position_on_plate_naive( center ):
# hardcoded approximate value for Palomar Schmidt 48" plates used for POSSI
plate_diameter_deg=6.6
# radius
r_max = plate_diameter_deg / 2.0
# random radius and angle
r = r_max * np.sqrt(np.random.random())
theta = 2 * np.pi * np.random.random()
#
xi = r * np.cos(theta) * u.deg
eta = r * np.sin(theta) * u.deg
off = SkyOffsetFrame(origin=center)
pos = SkyCoord(xi, eta, frame=off ).icrs
return pos
def extract_data(fits_files) :
extracted_data = []
for fits_file in fits_files:
with fits.open(fits_file) as hdul:
header = hdul[0].header
extracted_data.append({
'c_obj' : SkyCoord( # coordinate of the object (vanishing source)
ra=header["OBJCTRA"],
dec=header["OBJCTDEC"],
unit=(u.hourangle, u.deg),
frame='fk5',
equinox='J2000.0').transform_to('icrs'),
'c_plate' : SkyCoord( # coordinate of the center of the plate
ra=f"{header['PLTRAH']}h{header['PLTRAM']}m{header['PLTRAS']}s",
dec=f"{header['PLTDECSN']}{header['PLTDECD']}d{header['PLTDECM']}m{header['PLTDECS']}s",
unit=(u.hourangle, u.deg),
frame='fk5',
equinox='J2000.0').transform_to('icrs'),
'time' : Time( #observation time
time_str(
header["DATE-OBS"],
header["UT"]),
scale="ut1",
format="iso"),
'exposure' : float( header["EXPOSURE"] ) # length of exposure in minutes
} )
return extracted_data
def run_trial( data, randomize_coords=False) :
results = []
for obs in data:
t = randomize_time_within_exposure( obs[ 'time' ], obs[ 'exposure' ] )
c = randomize_position_on_plate_naive( obs[ 'c_plate' ] ) if randomize_coords else obs[ 'c_obj' ]
results.append( in_shadow( c, t )[ 0 ] )
return results
if __name__ == "__main__":
###################################################################
# N random trials
N = 1000
fits_files = glob.glob('fits_files/*.fits')
extracted_data = extract_data(fits_files)
# Random time within exposure window with real coords
real_percents = []
for i in tqdm(range(N), desc="Running real coordinate trials"):
real_results = run_trial(extracted_data, randomize_coords=False)
real_percents.append(np.mean(real_results) * 100)
real_percent_mean = np.mean( real_percents )
# Random time within exposure window with random coords on plate
rand_percents = []
for i in tqdm(range(N), desc="Running random coordinate trials"):
random_results = run_trial(extracted_data, randomize_coords=True)
rand_percents.append(np.mean(random_results) * 100)
rand_percent_mean = np.mean(rand_percents)
####################################################################
#Save the results
rand_percents = np.array( rand_percents ).astype( np.float64 )
rand_percents.tofile( "rand_percents.npy" )
real_percents = np.array( real_percents ).astype( np.float64 )
real_percents.tofile( "real_percents.npy" )
####################################################################
# significance
# Welch's t-test
t_stat, p_value_welch = stats.ttest_ind( real_percents, rand_percents, equal_var=False )
# Mann-Whitney U test
u_stat, p_value_mw = stats.mannwhitneyu(real_percents, rand_percents, alternative='two-sided')
# Cohen's d effect size
pooled_std = np.sqrt((np.var(real_percents, ddof=1) + np.var(rand_percents, ddof=1)) / 2)
cohens_d = (real_percent_mean - rand_percent_mean) / pooled_std
# Convert p-values to sigma
def p_to_sigma(p_value):
return stats.norm.ppf(1 - p_value/2)
sigma_welch = p_to_sigma(p_value_welch)
sigma_mw = p_to_sigma(p_value_mw)
print(f"Welch's t-test: t={t_stat:.4f}, p={p_value_welch:.6f} ({sigma_welch:.2f}σ)")
print(f"Mann-Whitney U: U={u_stat:.0f}, p={p_value_mw:.6f} ({sigma_mw:.2f}σ)")
print(f"Cohen's d: {cohens_d:.4f}")
####################################################################
# Plot the results
N_BINS = int( np.sqrt( 1.5*N ) )
bins = np.linspace(
min( np.min( rand_percents ), np.min( real_percents ) ),
max( np.max( rand_percents ), np.max( real_percents ) ),
N_BINS )
fig, ax = plt.subplots(figsize=(12, 8))
ax.hist(rand_percents, bins, color="seagreen", alpha=0.5, label='Random Coord Mean Shadow Percent' )
ax.hist(real_percents, bins, color="royalblue", alpha=0.5, label='Real Coord Mean Shadow Percent' )
ax.axvline(x=rand_percent_mean, color='green', linestyle='-', linewidth=2.5, label=f'Random Coord Mean: {rand_percent_mean:.3f}%')
ax.axvline(x=real_percent_mean, color='blue', linestyle='-', linewidth=2.5, label=f'Real Coord Mean: {real_percent_mean:.3f}%')
# Add statistics text box
stats_text = f"Welch's t-test: t={t_stat:.3f}, p={p_value_welch:.4f} ({sigma_welch:.2f}σ)\n"
stats_text += f"Mann-Whitney U: U={u_stat:.0f}, p={p_value_mw:.4f} ({sigma_mw:.2f}σ)\n"
stats_text += f"Cohen's d: {cohens_d:.3f}"
ax.text(0.95, 0.95, stats_text, transform=ax.transAxes, fontsize=11,
verticalalignment='top', horizontalalignment='left',
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
ax.set_title(f'Distribution Comparison for {N} Trials', fontsize=16)
ax.set_xlabel('Overall Percentage In Earth\'s Shadow (%)', fontsize=12)
ax.set_ylabel('Number of Trials', fontsize=12)
ax.set_ylim(0, ax.get_ylim()[1] * 1.2)
ax.legend(fontsize=10)
plt.savefig( "random_plate_pos_v_real_plate_pos.png" )
plt.show()
. I don't want to use my work laptop for this, but I want to get involved with your code 
