celestial-navigation/starfixdata_stat_1_na.py at main · alinnman/celestial-navigation · GitHub

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
''' This is a sample for celestial navigation for a stationary observer
    © August Linnman, 2025, email: august@linnman.net
    MIT License (see LICENSE file)

    This sample uses an algorithm for better accuracy using repeated
    refinements of the DR position.
'''

from time import time
from starfix import Sight, SightCollection, get_representation,\
                    get_google_map_string, IntersectError, LatLonGeodetic,\
                    show_or_display_file


def get_starfixes (drp_pos : LatLonGeodetic,
                   alt_sigma : float = 0.0,
                   time_sigma : float = 0.0) -> SightCollection :
    ''' Returns a list of used star fixes (SightCollection) '''

    Sight.set_estimated_position (drp_pos)
    Sight.set_alt_diff           (alt_sigma)
    Sight.set_time_diff          (time_sigma)

    temperature = 10
    dt_dh = -0.01

    a = Sight (object_name          = "Sun",
               set_time             = "2024-05-05 15:55:18+00:00",
               gha_time_0           = "45:50.4",
               gha_time_1           = "60:50.4",
               decl_time_0          = "16:30.6",
               decl_time_1          = "16:31.3",
               measured_alt         = "55:8:1.1",
               temperature          = temperature,
               dt_dh                = dt_dh
               )

    b = Sight (object_name          = "Sun",
               set_time             = "2024-05-05 23:01:19+00:00",
               gha_time_0           = "165:50.8",
               gha_time_1           = "180:50.8",
               decl_time_0          = "16:36.2",
               decl_time_1          = "16:36.9",
               measured_alt         = "19:28:19",
               temperature          = temperature,
               dt_dh                = dt_dh
               )

    c = Sight (object_name          = "Vega",
               set_time             = "2024-05-06 04:04:13+00:00",
               gha_time_0           = "284:30.4",
               gha_time_1           = "299:32.9",
               decl_time_0          = "38:48.1",
               measured_alt         = "30:16:23.7",
               sha_diff             = "80:33.4",
               temperature          = temperature,
               dt_dh                = dt_dh
               )
    return SightCollection ([a, b, c])

def main ():
    ''' Main body of script.'''

    starttime = time ()
    # the_pos = LatLonGeodetic (-40, 90) # BAD DRP position
    the_pos = LatLonGeodetic (40, -90) # Rough DRP position
    assume_good_pos = True
    # The exact position is 41°51'00.1"N 87°39'00.2"W
    intersections = collection = taken_ms = None
    try:
        intersections, _, _, collection, _ =\
              SightCollection.get_intersections_conv (return_geodetic=True,
                                                      estimated_position=the_pos,
                                                      get_starfixes=get_starfixes,
                                                      assume_good_estimated_position=\
                                                      assume_good_pos)
        assert intersections is not None
        assert collection is not None
        endtime = time ()
        taken_ms = round((endtime-starttime)*1000,3)
        print (get_representation(intersections,1))
        assert isinstance (intersections, LatLonGeodetic)
        print ("Google Maps Coordinate = " + get_google_map_string(intersections,4))

        # Check azimuth
        assert isinstance (intersections, LatLonGeodetic)
        counter = 0
        for s in collection.get_sf_list ():
            counter += 1
            az = s.get_azimuth (intersections)
            print ("Azimuth " + str(counter) + " = " + str(round(az,2)))

        # Diagnostics for map rendering etc.
        print ("Some useful data follows")
        counter = 0
        for s in collection.get_sf_list ():
            assert isinstance (s, Sight)
            counter += 1
            print (str(counter) + " radius = " +\
                    str(round(s.get_circle(geodetic=True).get_radius (),1)))
            print (str(counter) + " GP     = " +\
                    get_google_map_string(LatLonGeodetic(ll=s.get_gp()),4))

    except IntersectError as ve:
        print ("Cannot perform a sight reduction. Bad sight data.\n" + str(ve))
        if ve.coll_object is not None:
            if isinstance (ve.coll_object, SightCollection):
                collection = ve.coll_object

    if collection is not None and not isinstance (intersections, tuple):
        the_map = collection.render_folium (intersections, draw_grid=False)
        file_name = "./map.html"
        the_map.save (file_name)
        show_or_display_file (file_name)

    if taken_ms is not None:
        print ("Time taken = " +str(taken_ms)+" ms")

if __name__ == '__main__':
    main()