-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathcentralProgram.py
executable file
·171 lines (126 loc) · 7.69 KB
/
centralProgram.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
from mpl_toolkits import mplot3d
import numpy as np
import datetime
import matplotlib.pyplot as plt
from sideFunctions import time2Julian, earth2Sun, sun2Moon
from sideFunctions import leapYearAdjuster
from eclipsePlotter import EclipsePlotter
import math
import time
from angleCalculations import AngleCalculations as angleCalc
# this gets the position vectors of the stellar bodies that we are using
from jplephem.spk import SPK
kernel = SPK.open('de440.bsp')
yearEclipse = []
possibleSolarEclipseDates = np.zeros((13, 1000))
# trying to plot a line from center of sun to center of moon,
# as well as two lines which approximate the diameter of the earth
VecStart_x = [0 for rows in range(20)]
VecEnd_x = [0 for rows in range(20)]
VecStart_y = [0 for rows in range(20)]
VecEnd_y = [0 for rows in range(20)]
VecStart_z = [0 for rows in range(20)]
VecEnd_z = [0 for rows in range(20)]
VecConnection_x = [0 for rows in range(20)]
VecConnection_y = [0 for rows in range(20)]
VecConnection_z = [0 for rows in range(20)]
m = 0
# + These for loops check each date to see whether it meets the criteria below,
# basically that the moons shadow hits somewhere within the cross section of the earth.
# + I have been selecting a specific year in order to compare to resources online to see
# the accuracy of my model.
for TestYear in range(2023, 2024, 1):
for TestMonth in range(10, 11, 1):
counter = 0
print(TestMonth)
# this function handles if its a leap year, adjusting the days in that particular month
dayUpperRange = leapYearAdjuster(TestYear, TestMonth)
for TestDay in range(1, dayUpperRange, 1):
for TestHour in range(1, 24, 1):
for TestMinute in range(1, 60, 1):
# This gets the julian date for that specific time of year
julianTime = time2Julian(TestYear, TestMonth, TestDay, TestHour, TestMinute)
# position of earths center w.r.t. the sun
position1 = earth2Sun(julianTime)
xEarthCent = position1[0]
yEarthCent = position1[1]
zEarthCent = position1[2]
# position of moon w.r.t. the sun
position2 = sun2Moon(julianTime)
xMoonCent = position2[0]
yMoonCent = position2[1]
zMoonCent = position2[2]
# radius of the respective body (km)
r_earth = 6378
r_moon = 1737
# straight line from center of sun to center of moon
VecStart_x[0] = 0
VecStart_y[0] = 0
VecStart_z[0] = 0
VecEnd_x[0] = xMoonCent
VecEnd_y[0] = yMoonCent
VecEnd_z[0] = zMoonCent
# +line that is perpendicular to the line connecting the sun and the earth,
# parallel with the xy plane (horizontal)
# originally centered at the sun but I've translated it by adding x1,y1,z1 to
# it in order to get it centered on the earth
VecStart_x[1] = xEarthCent+r_earth
VecStart_y[1] = yEarthCent
VecStart_z[1] = zEarthCent
VecEnd_x[1] = xEarthCent - r_earth
VecEnd_y[1] = yEarthCent
VecEnd_z[1] = zEarthCent
# +line that is perpendicular to the line connecting the sun and the earth,
# parallel with the xz-plane (vertical)
# +originally centered at the sun but I've translated it by adding x1,y1,z1 to
# it in order to get it centered on the earth
VecStart_x[2] = xEarthCent
VecStart_y[2] = yEarthCent
VecStart_z[2] = zEarthCent + r_earth
VecEnd_x[2] = xEarthCent
VecEnd_y[2] = yEarthCent
VecEnd_z[2] = zEarthCent - r_earth
# + line that is perpendicular to the line connecting the sun and the moon,
# and I have translated it to the moon
# + this is horizontal line(parallel with the xz-plane)
VecStart_x[3] = xMoonCent + r_moon
VecStart_y[3] = yMoonCent
VecStart_z[3] = zMoonCent
VecEnd_x[3] = xMoonCent - r_moon
VecEnd_y[3] = yMoonCent
VecEnd_z[3] = zMoonCent
# + this is vertical line(parallel with the yz-plane) at the moon
vecStart4, vecEnd4 = [0, 0, 0], [0, 0, 0]
VecStart_x[4] = xMoonCent
VecStart_y[4] = yMoonCent
VecStart_z[4] = zMoonCent + r_moon
VecEnd_x[4] = xMoonCent
VecEnd_y[4] = yMoonCent
VecEnd_z[4] = zMoonCent - r_moon
# these are the angle theta from spherical coord.
# theta plus and minus are angles to get from one side of the earth to the other,horizontally
thetaPlusEarth = angleCalc(VecEnd_x[1], VecEnd_y[1], VecEnd_z[1]).theta()
thetaMinusEarth = angleCalc(VecStart_x[1], VecStart_y[1], VecStart_z[1]).theta()
# lines that are being drawn from center of sun to horizontal edges of the moon
thetaSun2Moon1 = angleCalc(VecEnd_x[3], VecEnd_y[3], VecEnd_z[3]).theta()
thetaSun2Moon2 = angleCalc(VecStart_x[3], VecStart_y[3], VecStart_z[3]).theta()
# phi plus and minus are the upper and lower points on the edges of the earth
phiPlusEarth = angleCalc(VecEnd_x[2], VecEnd_y[2], VecEnd_z[2]).phi()
phiMinusEarth = angleCalc(VecStart_x[2], VecStart_y[2], VecStart_z[2]).phi()
# lines that are being drawn from center of sun to vertical edges of the moon
phiSun2Moon1 = angleCalc(VecEnd_x[4], VecEnd_y[4], VecEnd_z[4]).phi()
phiSun2Moon2 = angleCalc(VecStart_x[4], VecStart_y[4], VecStart_z[4]).phi()
# +this checks whether the angle of the line from the sun's center to the moons horizontal
# and vertical edges, where I have split the moon's cross section into two perpendicular lines,
# is within the range of the angle from left/right
# and top/bottom parts of the earth, same thing being done to the earth as was done to the moon.
if ((phiMinusEarth <= phiSun2Moon1 <= phiPlusEarth) or (phiMinusEarth <= phiSun2Moon2 <= phiPlusEarth)) and \
(((thetaPlusEarth <= thetaSun2Moon1 <= thetaMinusEarth) or (thetaPlusEarth <= thetaSun2Moon2 <= thetaMinusEarth)) or \
((thetaMinusEarth <= thetaSun2Moon1 <= thetaPlusEarth) or (thetaMinusEarth <= thetaSun2Moon2 <= thetaPlusEarth))):
print(TestYear, TestMonth, TestDay, TestHour, TestMinute)
# this filters out the solar from the lunar eclipses
if (((xMoonCent**2+yMoonCent**2+zMoonCent)**(1/2)) < ((xEarthCent**2+yEarthCent**2+zEarthCent)**(1/2))):
possibleSolarEclipseDates[TestMonth][counter] = (julianTime)
counter += 1
lat, lon = EclipsePlotter(possibleSolarEclipseDates).latANDlonFinder()
blurb = EclipsePlotter(possibleSolarEclipseDates).plotter(lat, lon)