celestial-navigation/starfixdata_stat_2_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
''' 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,
                   time_sigma : float = 0.0,
                   alt_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)

    a = Sight (object_name          = "Capella",
               set_time             = "2024-09-17 23:36:13+00:00",
               #gha_time_0           = "342:21.9",
               #gha_time_1           = "357:24.4",
               #decl_time_0          = "46 :1.2",
               #decl_time_1          = "46 :1.3",
               #sha_diff             = "280:22.3",
               measured_alt         = "33 :9    :34"
              )

    b = Sight (object_name          = "Moon",
               set_time             = "2024-09-17 23:41:13+00:00",
               #gha_time_0           = "347:55.7" ,
               #gha_time_1           = "2  :24.6",
               #decl_time_0          = "-3 :43.5",
               #decl_time_1          = "-3 :25.3",
               #horizontal_parallax  = 61.2,
               measured_alt         = "48 :22  :5.2"
              )

    c = Sight (object_name          = "Vega",
               set_time             = "2024-09-17 23:46:13+00:00",
               #gha_time_0           = "342:21.9",
               #gha_time_1           = "357:24.4",
               #decl_time_0          = "38 :48.6",
               #sha_diff             = "80 :33.3",
               measured_alt         = "25 :39:4"
              )
    return SightCollection ([a, b, c])

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

    starttime = time ()
    the_pos = LatLonGeodetic (35, 10) # Rough DRP position
    # The exact position is 36° 45' 11.01", 10° 13' 8.00"

    intersections = collection = taken_ms = None
    try:
        intersections, _, _, collection, _=\
              SightCollection.get_intersections_conv (return_geodetic=True,
                                                      estimated_position=the_pos,
                                                      get_starfixes=get_starfixes)
        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)
        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()