-
Notifications
You must be signed in to change notification settings - Fork 0
/
cmpToDayMonthMethod.py
122 lines (99 loc) · 4.41 KB
/
cmpToDayMonthMethod.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
#!/usr/bin/env python3
#
# -*- coding: utf-8 -*-
#
# Copyright (C) 2019 CS Systèmes d'Information (CS SI)
#
# This file is part of the "Calceph Earth-Sun distance" project
#
# https://github.com/CS-SI/calceph_earth-sun_distance
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
import numpy as np
import math
import argparse
import matplotlib.pyplot as plt
from calcephpy import CalcephBin, Constants
def to_julian_date(year, month, day, hour, minute, second):
a = (14 - month) / 12
y = year + 4800 - a
m = month + 12 * a - 3
julian_date = day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045
# float division
julian_date += (hour-12) / 24. + minute / 1440. + second / 86400.
return julian_date
def get_solar_distance_from_day_month(day, month):
if month <= 2:
j = (month - 1) * 31 + day
elif month > 8:
j = (month - 1) * 31 - (month - 2) / 2 - 2 + day
else:
j = (month - 1) * 31 - (month - 1) / 2 - 2 + day
om = (j - 4) * .9856 * math.pi/180
d__1 = 1. - math.cos(om) * .01673
return d__1
if __name__ == "__main__":
parser = argparse.ArgumentParser(
description='''cmpToMonthDayMethod script help section. ''')
parser.add_argument('ephemeris_files', type=str, nargs='+', help='Paths to the ephemeris files to use. They must '
'be compatible with the calceph library. (see '
'https://www.imcce.fr/content/medias/recherche/'
'equipes/asd/calceph/html/python/calceph.multiple.'
'html#menu-calceph-open)')
parser.add_argument('year', type=int, help='Year at the moment of the acquisition.')
args = parser.parse_args()
# open ephemeris file
try:
peph = CalcephBin.open(args.ephemeris_files)
except RuntimeError:
raise Exception('Cannot open ephemeris files')
solar_dist_eph = []
solar_dist_day_month = []
diff = []
rel_diff = []
for cur_month in range(12):
nb_days = 31
if cur_month == 1:
nb_days = 28
elif cur_month % 2 == 1:
nb_days = 30
for cur_day in range(nb_days):
# compute julian date
cur_julian_date = to_julian_date(args.year, cur_month+1, cur_day+1, 0, 0, 0)
julian_date_ent = int(cur_julian_date)
julian_date_dec = cur_julian_date - julian_date_ent
# compute earth/sun distance using the calceph library
out = peph.compute_unit(julian_date_ent, julian_date_dec, 11, 3, Constants.UNIT_AU + Constants.UNIT_SEC)
norm_eph = math.sqrt(pow(out[0], 2.) + pow(out[1], 2.) + pow(out[2], 2.))
# compute earth/sun distance using the elliptical model
norm_day_month = get_solar_distance_from_day_month(cur_day+1, cur_month+1)
solar_dist_eph.append(norm_eph)
solar_dist_day_month.append(norm_day_month)
diff.append(norm_eph-norm_day_month)
rel_diff.append(abs(norm_eph-norm_day_month)/norm_day_month*100.)
# plot the results
plt.plot(solar_dist_eph, 'r', label='using ephemeris')
plt.plot(solar_dist_day_month, 'b', label='using day/month')
plt.ylabel('solar distance (in AU)')
plt.legend(loc='upper right')
plt.show()
# print the errors statistics
print('difference max : {}'.format(np.max(np.abs(diff))))
print('difference mean : {}'.format(np.mean(diff)))
print('difference standard deviation : {}'.format(np.std(diff)))
print('relative difference max (in %): {}'.format(np.max(rel_diff)))
print('relative difference mean (in %): {}'.format(np.mean(rel_diff)))
print('relative difference standard deviation (in %): {}'.format(np.std(rel_diff)))
# close ephemeris file
peph.close()