My agentic slop goes here. Not intended for anyone else!
anil.recoil.org
1(** Coordinate pair module *)
2
3type t = {
4 lat : Lat.t;
5 lon : Lon.t;
6}
7
8let create ~lat ~lon = {
9 lat = Lat.create lat;
10 lon = Lon.create lon;
11}
12
13let of_lat_lon lat lon = { lat; lon }
14
15let lat t = t.lat
16let lon t = t.lon
17let lat_float t = Lat.to_float t.lat
18let lon_float t = Lon.to_float t.lon
19
20(** Constants *)
21let earth_radius_km = 6371.0
22let deg_to_rad = Float.pi /. 180.0
23let rad_to_deg = 180.0 /. Float.pi
24
25(** Distance calculations *)
26
27let distance_haversine p1 p2 =
28 let lat1 = lat_float p1 *. deg_to_rad in
29 let lat2 = lat_float p2 *. deg_to_rad in
30 let lon1 = lon_float p1 *. deg_to_rad in
31 let lon2 = lon_float p2 *. deg_to_rad in
32
33 let dlat = lat2 -. lat1 in
34 let dlon = lon2 -. lon1 in
35
36 let a =
37 (sin (dlat /. 2.0) ** 2.0) +.
38 cos lat1 *. cos lat2 *. (sin (dlon /. 2.0) ** 2.0) in
39
40 let c = 2.0 *. atan2 (sqrt a) (sqrt (1.0 -. a)) in
41 earth_radius_km *. c
42
43let distance_vincenty p1 p2 =
44 (* WGS84 ellipsoid constants *)
45 let a = 6378.137 in (* semi-major axis in km *)
46 let b = 6356.752 in (* semi-minor axis in km *)
47 let f = 1.0 /. 298.257223563 in (* flattening *)
48
49 let lat1 = lat_float p1 *. deg_to_rad in
50 let lat2 = lat_float p2 *. deg_to_rad in
51 let lon1 = lon_float p1 *. deg_to_rad in
52 let lon2 = lon_float p2 *. deg_to_rad in
53
54 let l = lon2 -. lon1 in
55 let u1 = atan ((1.0 -. f) *. tan lat1) in
56 let u2 = atan ((1.0 -. f) *. tan lat2) in
57 let sin_u1 = sin u1 in
58 let cos_u1 = cos u1 in
59 let sin_u2 = sin u2 in
60 let cos_u2 = cos u2 in
61
62 let rec iterate lambda _prev_lambda iter =
63 if iter > 100 then
64 (* Failed to converge, fall back to haversine *)
65 distance_haversine p1 p2
66 else
67 let sin_lambda = sin lambda in
68 let cos_lambda = cos lambda in
69 let sin_sigma =
70 sqrt ((cos_u2 *. sin_lambda) ** 2.0 +.
71 (cos_u1 *. sin_u2 -. sin_u1 *. cos_u2 *. cos_lambda) ** 2.0) in
72
73 if sin_sigma = 0.0 then 0.0 (* coincident points *)
74 else
75 let cos_sigma = sin_u1 *. sin_u2 +. cos_u1 *. cos_u2 *. cos_lambda in
76 let sigma = atan2 sin_sigma cos_sigma in
77 let sin_alpha = cos_u1 *. cos_u2 *. sin_lambda /. sin_sigma in
78 let cos_sq_alpha = 1.0 -. sin_alpha ** 2.0 in
79 let cos_2sigma_m =
80 if cos_sq_alpha = 0.0 then 0.0 (* equatorial line *)
81 else cos_sigma -. 2.0 *. sin_u1 *. sin_u2 /. cos_sq_alpha in
82
83 let c = f /. 16.0 *. cos_sq_alpha *. (4.0 +. f *. (4.0 -. 3.0 *. cos_sq_alpha)) in
84 let lambda' = l +. (1.0 -. c) *. f *. sin_alpha *.
85 (sigma +. c *. sin_sigma *.
86 (cos_2sigma_m +. c *. cos_sigma *.
87 (-1.0 +. 2.0 *. cos_2sigma_m ** 2.0))) in
88
89 if abs_float (lambda' -. lambda) < 1e-12 then
90 (* Converged *)
91 let u_sq = cos_sq_alpha *. (a ** 2.0 -. b ** 2.0) /. (b ** 2.0) in
92 let cap_a = 1.0 +. u_sq /. 16384.0 *.
93 (4096.0 +. u_sq *. (-768.0 +. u_sq *. (320.0 -. 175.0 *. u_sq))) in
94 let cap_b = u_sq /. 1024.0 *. (256.0 +. u_sq *. (-128.0 +. u_sq *. (74.0 -. 47.0 *. u_sq))) in
95 let delta_sigma = cap_b *. sin_sigma *.
96 (cos_2sigma_m +. cap_b /. 4.0 *.
97 (cos_sigma *. (-1.0 +. 2.0 *. cos_2sigma_m ** 2.0) -.
98 cap_b /. 6.0 *. cos_2sigma_m *.
99 (-3.0 +. 4.0 *. sin_sigma ** 2.0) *.
100 (-3.0 +. 4.0 *. cos_2sigma_m ** 2.0))) in
101
102 b *. cap_a *. (sigma -. delta_sigma)
103 else
104 iterate lambda' lambda (iter + 1)
105 in
106
107 iterate l 0.0 0
108
109(** Bearing calculations *)
110
111let bearing p1 p2 =
112 let lat1 = lat_float p1 *. deg_to_rad in
113 let lat2 = lat_float p2 *. deg_to_rad in
114 let lon1 = lon_float p1 *. deg_to_rad in
115 let lon2 = lon_float p2 *. deg_to_rad in
116
117 let dlon = lon2 -. lon1 in
118 let y = sin dlon *. cos lat2 in
119 let x = cos lat1 *. sin lat2 -. sin lat1 *. cos lat2 *. cos dlon in
120
121 let bearing = atan2 y x *. rad_to_deg in
122 mod_float (bearing +. 360.0) 360.0
123
124let bearing_final p1 p2 =
125 (* Final bearing is the reverse of the initial bearing from p2 to p1 *)
126 let reverse_bearing = bearing p2 p1 in
127 mod_float (reverse_bearing +. 180.0) 360.0
128
129(** Geographic operations *)
130
131let offset t ~heading ~distance =
132 let lat1 = lat_float t *. deg_to_rad in
133 let lon1 = lon_float t *. deg_to_rad in
134 let brng = heading *. deg_to_rad in
135 let d = distance /. earth_radius_km in
136
137 let lat2 = asin (sin lat1 *. cos d +. cos lat1 *. sin d *. cos brng) in
138 let lon2 = lon1 +. atan2 (sin brng *. sin d *. cos lat1)
139 (cos d -. sin lat1 *. sin lat2) in
140
141 create ~lat:(lat2 *. rad_to_deg) ~lon:(lon2 *. rad_to_deg)
142
143let midpoint p1 p2 =
144 let lat1 = lat_float p1 *. deg_to_rad in
145 let lat2 = lat_float p2 *. deg_to_rad in
146 let lon1 = lon_float p1 *. deg_to_rad in
147 let lon2 = lon_float p2 *. deg_to_rad in
148
149 let dlon = lon2 -. lon1 in
150 let bx = cos lat2 *. cos dlon in
151 let by = cos lat2 *. sin dlon in
152
153 let lat3 = atan2 (sin lat1 +. sin lat2)
154 (sqrt ((cos lat1 +. bx) ** 2.0 +. by ** 2.0)) in
155 let lon3 = lon1 +. atan2 by (cos lat1 +. bx) in
156
157 create ~lat:(lat3 *. rad_to_deg) ~lon:(lon3 *. rad_to_deg)
158
159let interpolate p1 p2 fraction =
160 if fraction <= 0.0 then p1
161 else if fraction >= 1.0 then p2
162 else
163 let lat1 = lat_float p1 in
164 let lat2 = lat_float p2 in
165 let lon1 = lon_float p1 in
166 let lon2 = lon_float p2 in
167
168 (* Simple linear interpolation in degrees *)
169 let lat = lat1 +. (lat2 -. lat1) *. fraction in
170 let lon = lon1 +. (lon2 -. lon1) *. fraction in
171
172 create ~lat ~lon
173
174(** Comparison *)
175
176let equal p1 p2 =
177 Lat.equal p1.lat p2.lat && Lon.equal p1.lon p2.lon
178
179let ( = ) = equal
180
181let almost_equal ?(epsilon = 0.000001) p1 p2 =
182 distance_haversine p1 p2 < epsilon
183
184(** Formatting *)
185
186let to_string t =
187 Printf.sprintf "%s, %s" (Lat.to_string t.lat) (Lon.to_string t.lon)
188
189let to_dms_string t =
190 Printf.sprintf "%s, %s" (Lat.to_dms_string t.lat) (Lon.to_dms_string t.lon)
191
192let format ?(lat_precision = 6) ?(lon_precision = 6) t =
193 Printf.sprintf "%s, %s"
194 (Lat.format ~precision:lat_precision t.lat)
195 (Lon.format ~precision:lon_precision t.lon)