-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathrefview_struct.jl
2158 lines (1838 loc) · 73.3 KB
/
refview_struct.jl
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
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
### A Pluto.jl notebook ###
# v0.19.25
#> custom_attrs = ["hide-enabled"]
using Markdown
using InteractiveUtils
# ╔═╡ 7d23f727-1907-4965-a940-cc873f6b2191
using PlutoDevMacros
# ╔═╡ 73e00cef-9734-439a-b89b-7c1d99aab74e
# ╠═╡ skip_as_script = true
#=╠═╡
begin
using BenchmarkTools
using PlutoTest
using PlutoExtras
using PlutoPlotly
using PlutoUI
end
╠═╡ =#
# ╔═╡ 7729ce27-df74-4393-ab70-c4e2864c85f5
@fromparent begin
import *
using >.Rotations
using >.StaticArrays
using >.LinearAlgebra
using >.DocStringExtensions
end
# ╔═╡ 3bda5426-c0de-493f-9514-30b6fe762463
# ╠═╡ skip_as_script = true
#=╠═╡
md"""
# Packages
"""
╠═╡ =#
# ╔═╡ dcc81988-903b-4707-a70c-09c38682c80f
# ╠═╡ skip_as_script = true
#=╠═╡
ExtendedTableOfContents()
╠═╡ =#
# ╔═╡ 3fd1046c-fabf-4264-9638-ba41301b1804
# ╠═╡ skip_as_script = true
#=╠═╡
md"""
## load other notebook
"""
╠═╡ =#
# ╔═╡ a57e3983-21de-4a2e-a227-8265fee6b56b
# ╠═╡ skip_as_script = true
#=╠═╡
md"""
# Exports
"""
╠═╡ =#
# ╔═╡ a5112e66-c2a2-4ed2-9951-5f97bc1745d5
md"""
# ReferenceView
"""
# ╔═╡ 9e65ad17-95bd-46fa-bc18-9d5e8c501d9a
md"""
## face_rotation
"""
# ╔═╡ f1d1295e-7fa1-44d0-bdac-8b7830da8a61
@enum Faces begin
PositiveX = 1
PositiveY = 2
PositiveZ = 3
NegativeX = -1
NegativeY = -2
NegativeZ = -3
end
# ╔═╡ 9e8f0786-1daa-4b9e-9172-fc0767582c7e
begin
function to_face(s::Symbol)
s === :PositiveX && return PositiveX
s === :PositiveY && return PositiveY
s === :PositiveZ && return PositiveZ
s === :NegativeX && return NegativeX
s === :NegativeY && return NegativeY
s === :NegativeZ && return NegativeZ
error("Unrecognized")
end
to_face(n::Int) = Faces(n)
to_face(face::Faces) = face
end
# ╔═╡ 93222642-f2a6-4de7-8c92-21c96ef009a4
function face_rotation(face::Faces)::RotMatrix{3, Float64}
face === PositiveZ && return one(RotMatrix{3, Float64}) # No Rotation
face === NegativeZ && return RotY(180°)
face === PositiveX && return RotY(+90°)
face === NegativeX && return RotY(-90°)
face === PositiveY && return RotX(-90°)
face === NegativeY && return RotX(+90°)
end
# ╔═╡ cb0b070b-f70a-458e-bf72-0a4d2e93ec41
face_rotation(face) = face_rotation(to_face(face))
# ╔═╡ 28c868ca-9954-47ff-8949-a41fb7fc6d41
#=╠═╡
@benchmark Ref(face_rotation(:NegativeX))[]
╠═╡ =#
# ╔═╡ f4c3876a-a81b-42f3-870a-43526e4c116e
md"""
## boresight versor
"""
# ╔═╡ cc51ab70-0167-477f-a62d-88567e94fed9
begin
function boresight_versor(face::Faces)
face === PositiveZ && return SA_F64[0,0,1]
face === NegativeZ && return SA_F64[0,0,-1]
face === PositiveX && return SA_F64[1,0,0]
face === NegativeX && return SA_F64[-1,0,0]
face === PositiveY && return SA_F64[0,1,0]
face === NegativeY && return SA_F64[0,-1,0]
end
boresight_versor(face) = boresight_versor(to_face(face))
# convert SVectors and NTuple
boresight_versor(s::Union{SVector{3}, NTuple{3, <:Number}}) = SVector{3, Float64}(s) |> normalize
# Handle normal vectors
function boresight_versor(v::AbstractVector{<:Number})
@assert length(v) == 3 "Only Vectors with exactly three elements are supported for specifying the boresight versor. Consider using Tuples or SVector from StaticArrays"
# Revert to SVector method
boresight_versor(SVector{3, Float64}(v))
end
boresight_versor(x,y,z) = boresight_versor(SVector{3, Float64}(x,y,z))
end
# ╔═╡ c54622f5-835d-4957-bdc6-88afed9c9d2a
# ╠═╡ skip_as_script = true
#=╠═╡
boresight_versor(-3)
╠═╡ =#
# ╔═╡ 90811965-8bf0-4a4c-86f4-9bc5c86be71c
# ╠═╡ skip_as_script = true
#=╠═╡
boresight_versor((1,2,3))
╠═╡ =#
# ╔═╡ 278a536f-82dc-4f4e-8000-cc6a064ab2ee
# ╠═╡ skip_as_script = true
#=╠═╡
boresight_versor(1,2,3)
╠═╡ =#
# ╔═╡ 5dc71a85-a20d-4448-98b4-5065a249df1d
md"""
## Struct Definition
"""
# ╔═╡ c6ee08ba-3546-48ea-9801-edc00dfd25f0
begin
"""
rv = ReferenceView{T}(lla_or_ecef::Union{LLA, Point3D},earthmodel::EarthModel; kwargs...)
A `ReferenceView` instance `rv` tracks the instantaneous position and orientation of an object in space or on earth. It is used to provide convenience methods to assess pointing, distance and visibility between obejcts. It can also be used to extract 3D coordinates (LLA or ECEF) based on pointing angles from the `ReferenceView` instance.
When initializing `rv`, two convenience aliases can be used:
- `UserView === ReferenceView{:User}` can be used to represent users
- `SatView === ReferenceView{:Satellite}` can be used to represent satellites
The main difference between the two is that the default local CRS for computing pointing and visibilities is West-North-Down (+Z towards nadir) for `SatView` while it's East-North-Up (+Z towards Zenith) for `UserView`.
# Fields
$TYPEDFIELDS
When multiple `ReferenceView` objects have to be modelled/tracked, the `earthmodel` field of all of them should point to the same `EarthModel` instance, so it should be generated once and then passed to the constructor of all the ReferenceView instances to ensure that all the objects are referring to the same earth model.\\
Doing so will enforce the same Ellipsoid is shared between all satellites even when it is changed from one of the SatView instances.
To understand how the fields `rot_z`, `rot_y` and `rot_x` are used to model the
rotation applied to the default local CRS to the intended one, see the
documentation of [`crs_rotation`](@ref)
The possible values for the reference `face` (defined with respect to the axes
of the local CRS) used for pointing/visibility computations can be found looking
at [`change_reference_face!`](@ref)
See also: [`change_position!`](@ref), [`change_attitude!`](@ref),
[`change_reference_face!`](@ref), [`get_range`](@ref), [`get_pointing`](@ref),
[`get_lla`](@ref), [`get_ecef`](@ref), [`geod_inverse`](@ref),
[`get_distance_on_earth`](@ref), [`get_nadir_beam_diameter`](@ref),
[`ExtraOutput`](@ref), [`get_visibility`](@ref),
[`get_mutual_visibility`](@ref).
"""
Base.@kwdef mutable struct ReferenceView{T}
"ECEF coordinates of the current satellite position"
ecef::SVector{3,Float64}
"LLA coordinates of the current satellite position"
lla::LLA
"Reference EarthModel used for the projections and for the geodesic computations"
earthmodel::EarthModel
"Rotation Matrix to go from the default local CRS to the ECEF CRS"
R::RotMatrix3{Float64}
"Rotation angle [rad] about the z axis to bring the default local CRS to the current one"
rot_z::Float64 = 0
"Rotation angle [rad] about the y axis to bring the default local CRS to the current one"
rot_y::Float64 = 0
"Rotation angle [rad] about the x axis to bring the default local CRS to the current one"
rot_x::Float64 = 0
"Reference race of the satellite used for the pointing computations"
face::Faces = PositiveZ
end
# Version allow to specify face as Symbol
function ReferenceView{T}(ecef, lla, earthmodel, R, rot_z, rot_y, rot_x, face::Union{Symbol, Int}) where T
ReferenceView{T}(ecef, lla, earthmodel, R, rot_z, rot_y, rot_x, to_face(face))
end
const SatView = ReferenceView{:Satellite}
const UserView = ReferenceView{:User}
function _refview_rotation(s::Symbol)
s === :Satellite && return :ECEFfromWND
s === :User && return :ECEFfromENU
error("The provided type of ReferenceView $s is not supported")
end
# Custom constructor
function ReferenceView{T}(lla_or_ecef::Union{LLA, Point3D},em::EarthModel; kwargs...) where T
lla, ecef = if lla_or_ecef isa LLA
lla = lla_or_ecef
ecef = ECEFfromLLA(em.ellipsoid)(lla)
lla, ecef
else
ecef = lla_or_ecef
lla = LLAfromECEF(em.ellipsoid)(ecef)
lla, ecef
end
R = _rotation_matrix(_refview_rotation(T), lla.lat, lla.lon)
ReferenceView{T}(;ecef,lla,earthmodel=em,R, kwargs...)
end
end
# ╔═╡ a34ca715-e7f2-4fa7-ba94-8ad3f9b1f4cd
function Base.getproperty(sv::ReferenceView, name::Symbol)
if name ∈ (:ellipsoid, :geod)
return getfield(getfield(sv,:earthmodel),name)
else
return getfield(sv, name)
end
end
# ╔═╡ 227c0f3e-6e45-48ac-92d4-d1c1c730e4e0
function Base.setproperty!(sv::V, name::Symbol, x) where V <: ReferenceView
error("You can't change fields of $V directly, use change_position! or change_attitude!")
end
# ╔═╡ 6eb9424e-3dd2-46d4-b4d2-81596bb81668
# ╠═╡ skip_as_script = true
#=╠═╡
begin
em = EarthModel()
sv = SatView(LLA(0,0,600km), em)
end
╠═╡ =#
# ╔═╡ 5372f3fe-699a-4f00-8e8e-36cbea224963
#=╠═╡
let
lla = LLA(0,0,600km)
ecef = ECEFfromLLA(em.ellipsoid)(lla)
@benchmark SatView($lla, $em)
end
╠═╡ =#
# ╔═╡ 5f1fd82d-f441-4a1b-9840-773a8635d3db
#=╠═╡
@benchmark Ref(getproperty($sv, :geod))[]
╠═╡ =#
# ╔═╡ fb418edb-4937-4668-acf3-e22e3a818a99
md"""
## CRS / Attitude rotation
"""
# ╔═╡ 12f31d2c-1eac-47bc-bda2-a32395ff8265
"""
R = crs_rotation(;rot_x=0.0, rot_y=0.0, rot_z=0.0)
Create the z-y'-x'' _intrinsic_ [Tait-Bryan](https://en.wikipedia.org/wiki/Euler_angles#Tait%E2%80%93Bryan_angles) rotation matrix used to transform/rotate a starting CRS to another one.
This is used within [`ReferenceView`](@ref) to account for satellite attitude or user pointing when defining the final satellite/user centric CRS.
Taking as example a `SatView`, the default CRS without any attitude rotation is always the `WND` one centered on the satellite, with the `X` axis pointing towards West, `Y` axis towards north and `Z` axis towards nadir.
Let's call this default CRS as `WND` and let's assume to have a satellite in a LEO orbit with an inclination of `70` degrees.
Let's assume that the satellite is at its ascending node, and that we want to rotate the CRS of the satellite so that the `Y` axis points in the direction of propagation. To create this new CRS that we will call `XYZ`, we have to rotate the `WND` one by `20` degrees clock-wise around the `Z` axis of the `WND` CRS.
Assuming to have a point `xyz` in coordinates relative to the `XYZ` CRS, we can express the same point as `wnd` relative to the `WND` CRS with the following equation:
```
wnd = crs_rotation(;rot_z=20°) * xyz
```
"""
crs_rotation(;rot_x = 0.0, rot_y = 0.0, rot_z = 0.0) = RotZYX(rot_z,rot_y,rot_x)
# ╔═╡ 091e4ec2-ea9e-411e-8f39-73aeb73c0214
# ╠═╡ skip_as_script = true
#=╠═╡
md"""
## Change Position
"""
╠═╡ =#
# ╔═╡ 3c2781b8-e8df-446f-95ed-597112edabec
begin
"""
change_position!(rv::ReferenceView, ecef, lla::LLA, R, rot_z, rot_x, rot_y)
change_position!(rv::ReferenceView, ecef::Point3D, lla::LLA; rot_z = 0.0, rot_y = 0.0, rot_x = 0.0)
change_position!(rv::ReferenceView, lla::LLA; rot_z = 0.0, rot_y = 0.0, rot_x = 0.0)
change_position!(rv::ReferenceView, ecef::Point3D; rot_z = 0.0, rot_y = 0.0, rot_x = 0.0)
Change the position of a [`ReferenceView`](@ref) object `rv`, also returning as output the modified
`rv`. The function mutates the `ecef`, `lla`, `R`, `rot_z`, `rot_y`, `rot_x` fields of the `rv`
object with the values provided as arguments (when using the 1st method above).\\
For the remaining methods, the missing arguments are computed automatically.
One would normally use one of the last 3 methods so that the rotation matrix (and eventually either ECEF or LLA) are correctly computed by the function.
The first method avoids any computations but does not validate that provided inputs are
correct/consistent and refer to the same position in space. For this reason it should
only be used if those values are correctly pre-computed elsewhere and one wants to avoid the duplicate computations.
Note: `rot_z`, `rot_y` and `rot_x` are used to obtain the actual reference CRS
to use for the pointing computations. These angles are used to transform the
default CRS (which is WND for SatView and ENU for UserView) into the desired
CRS. This is done performing a z-y-x _intrinsic_ rotation of the default CRS following the
[Tait-Bryan
Definition](https://en.wikipedia.org/wiki/Euler_angles#Tait%E2%80%93Bryan_angles).
See also: [`ReferenceView`](@ref), [`get_range`](@ref), [`get_pointing`](@ref), [`get_lla`](@ref),
[`get_ecef`](@ref), [`get_ecef`](@ref).
"""
function change_position!(rv::ReferenceView, ecef, lla::LLA, R, rot_z, rot_y, rot_x)
setfield!(rv,:ecef,ecef)
setfield!(rv,:lla,lla)
setfield!(rv,:R,R)
setfield!(rv,:rot_z,rot_z)
setfield!(rv,:rot_y,rot_y)
setfield!(rv,:rot_x,rot_x)
return rv
end
function change_position!(rv::ReferenceView, ecef::StaticVector{3}, lla::LLA; rot_z = 0.0, rot_y = 0.0, rot_x = 0.0)
R = _rotation_matrix(:ECEFfromUV, lla.lat, lla.lon)
R *= crs_rotation(;rot_z, rot_y, rot_x)
change_position!(rv, ecef, lla, R, rot_z, rot_y, rot_x)
end
function change_position!(rv::ReferenceView, lla::LLA; kwargs...)
ecef = ECEFfromLLA(rv.ellipsoid)(lla)
change_position!(rv, ecef, lla;kwargs...)
end
function change_position!(rv::ReferenceView, ecef::StaticVector{3}; kwargs...)
lla = LLAfromECEF(rv.ellipsoid)(ecef)
change_position!(rv, ecef, lla;kwargs...)
end
end
# ╔═╡ b17a4660-3938-4816-8182-6feae471b50d
#=╠═╡
let
rot_z = rot_x = rot_y = 0.0
rv = SatView(LLA(0,0,600km), em)
new_lla = LLA(0.1,0,600km)
new_ecef = ECEFfromLLA(em.ellipsoid)(new_lla)
_R = _rotation_matrix(:ECEFfromUV, new_lla.lat, new_lla.lon)
R = _R * crs_rotation(;rot_z, rot_y, rot_x)
@benchmark change_position!($rv, $new_ecef, $new_lla, $R, $rot_z, $rot_y, $rot_x)
end
╠═╡ =#
# ╔═╡ b4777d24-5ea7-4524-8d2a-ac5d070d58d5
#=╠═╡
let
lla = LLA(0,0,600km)
@benchmark _rotation_matrix(:ECEFfromUV, $lla.lat, $lla.lon)
end
╠═╡ =#
# ╔═╡ cddcb68d-5400-4cb9-90ef-b3e149900f8a
#=╠═╡
let
# The benchmark below is wrong as phantom allocations are appearing from a call to _rotation_matrix (see cell above). The same does not happen outside of Pluto where this function call is just around 30ns
rot_z = rot_x = rot_y = 0.0
rv = SatView(LLA(0,0,600km), em)
new_lla = LLA(0.1,0,600km)
new_ecef = ECEFfromLLA(em.ellipsoid)(new_lla)
@benchmark change_position!($rv, $new_ecef, $new_lla)
end
╠═╡ =#
# ╔═╡ 64919bbb-3b52-48e1-a7ea-c0916c47994d
#=╠═╡
let
rv = SatView(LLA(0,0,600km), em)
xyz = SA_F64[0,100,0] # 100m towards the direction of propagation
wnd = crs_rotation(;rot_z = 20°) * xyz
end
╠═╡ =#
# ╔═╡ 5c1d0beb-c00c-40a8-b255-c887be99e706
md"""
## Change Attitude
"""
# ╔═╡ f127481d-c60e-43d0-ab9f-4d4f35984015
"""
change_attitude!(rv::ReferenceView; rot_z = 0.0, rot_y = 0.0, rot_x = 0.0)
Modify the attitude (local CRS orientation) of `rv` by recomputing its internal rotation matrix as a product of the default location based one and the attitude rotation based on `rot_z`, `rot_y` and `rot_x`.
See [`crs_rotation`](@ref) for details on how the attitude rotation matrix is defined and computed.
See also [`ReferenceView`](ref), [`change_position!`](@ref), [`change_Reference_face!`](@ref)
"""
function change_attitude!(rv::ReferenceView; rot_z = 0.0, rot_y = 0.0, rot_x = 0.0)
(;ecef, lla) = rv
change_position!(rv, ecef, lla ;rot_z, rot_y, rot_x)
end
# ╔═╡ 3a745709-8f5b-4f22-848a-2f9754ab27d8
#=╠═╡
let
sat_lla = LLA(0°, 0°, 600km)
rv = SatView(sat_lla, em)
wnd_R = rv.R
sat_ecef = rv.ecef
wnd2ecef = ECEFfromWND(rv.ecef, wnd_R, em.ellipsoid)
# Compute the ECEF position of a point 100km below the satellite
ref_ecef = wnd2ecef(SA_F64[0,0,100e3])
## Now we change the attitude of the satellite
# 90 degrees around the X axis means that the ref point should be towards the positive Y axis
change_attitude!(rv; rot_x = deg2rad(90))
@test WNDfromECEF(rv.ecef, rv.R', em.ellipsoid)(ref_ecef) ≈ SA_F64[0,100e3,0]
# 90 degrees around the Y axis means that the ref point should be towards -X
change_attitude!(rv; rot_y = deg2rad(90))
@test WNDfromECEF(rv.ecef, rv.R', em.ellipsoid)(ref_ecef) ≈ SA_F64[-100e3,0,0]
# 90 degrees around the Z and 90 degrees around the Y axis means that the ref point should be towards -X
change_attitude!(rv; rot_z = deg2rad(90), rot_y = deg2rad(90))
@test WNDfromECEF(rv.ecef, rv.R', em.ellipsoid)(ref_ecef) ≈ SA_F64[-100e3,0,0]
# 45 degrees around the Z and 45 degrees around the Y axis means that the ref point should be at 45 degrees between the -X and +Z
change_attitude!(rv; rot_z = deg2rad(45), rot_y = deg2rad(45))
@test WNDfromECEF(rv.ecef, rv.R', em.ellipsoid)(ref_ecef) ≈ SA_F64[-√(100e3^2/2),0,√(100e3^2/2)]
end
╠═╡ =#
# ╔═╡ 851f1994-22f4-4347-b106-2fa3ea75ebf2
md"""
## Change Reference Face
"""
# ╔═╡ 74d39d15-8d6c-4fb8-8d04-18f32c393ad4
"""
change_reference_face!(rv::ReferenceView, face)
Modify the reference face that is used to compute the visibilities from object `rv`.
The face can be specified using an instance of `TelecomUtils.Faces` (not exported), a Symbol or an Int as follows:
- TelecomUtils.PositiveX **or** :PositiveX **or** 1
- TelecomUtils.PositiveY **or** :PositiveY **or** 2
- TelecomUtils.PositiveZ **or** :PositiveZ **or** 3
- TelecomUtils.NegativeX **or** :NegativeX **or** -1
- TelecomUtils.NegativeY **or** :NegativeY **or** -2
- TelecomUtils.NegativeZ **or** :NegativeZ **or** -3
See also [`ReferenceView`](ref), [`change_position!`](@ref), [`change_attitude!`](@ref)
"""
function change_reference_face!(rv::ReferenceView, face)
setfield!(rv, :face, to_face(face))
return rv
end
# ╔═╡ 8b68c7c2-5cbd-4edd-9739-0d2e8c7f5449
#=╠═╡
let
rv = SatView(LLA(0,0,500km), em)
change_reference_face!(rv, :NegativeZ)
end
╠═╡ =#
# ╔═╡ 84769564-8ba8-46f5-b494-b0689d9abd65
# ╠═╡ skip_as_script = true
#=╠═╡
md"""
## Get Range
"""
╠═╡ =#
# ╔═╡ 0cde0a71-7f27-4290-88cd-2cccf627926b
begin
"""
get_range(rv::ReferenceView,uv::Point2D[, ::ExtraOutput]; h = 0.0, face = rv.face, R = nothing)
get_range(rv::ReferenceView,target::Union{LLA, Point3D, ReferenceView}[, ::ExtraOutput]; face = rv.face, R = nothing)
Get the range [in m] between the reference view `rv` (Satellite or User) and a target point. The target point can be identified in two ways:
- Providing the uv pointing `uv` and the reference altitude `h` [m] above the ellipsoid (first method)
- directly passing in a point expressed either as `LLA`, as a `Point3D` representing its ECEF coordinates, or as another instance of `ReferenceView` (second method).
The pointing is assumed to be referred to the specified `face` of the
`ReferenceView`. When the location identified by the provided pointing is not
visible either because it's blocked by earth or because it's in a direction not
visible from the specified `face`, `NaN` is returned.
See [`change_reference_face!`](@ref) for supported identification of the `face` kwarg.
When called with an instance of `TelecomUtils.ExtraOutput` as last argument, the
function also returns the coordinated of the identified point in the local CRS
of `rv`.
The kwarg `R` represents the 3D Rotation Matrix that translates a vector from
ECEF coordinates to the coordinates of the desired local CRS around `rv`. By
default (if `R === nothing`) this rotation matrix is computed based on the
rotation matrix of the `rv` object and on the selected reference face.
See also: [`ReferenceView`](@ref), [`get_era`](@ref), [`get_mutual_pointing`](@ref), [`get_pointing`](@ref), [`get_lla`](@ref), [`get_ecef`](@ref), [`get_distance_on_earth`](@ref).
"""
function get_range(sv::ReferenceView,uv::Point2D, eo::ExtraOutput; h = 0.0, face = sv.face, R = nothing)
_R = isnothing(R) ? inv(sv.R * face_rotation(face)) : R
# Find the ecef coordinate of the target point on earth
ecef, r = ECEFfromUV(sv.ecef, _R', sv.ellipsoid)(uv, eo;h)
# Compute the point coordinates in the ReferenceView local CRS
xyz = r * SVector(uv[1], uv[2], sqrt(1 - (uv[1]^2 + uv[2]^2)))
# Return the distance betwteen the satellite and the point
return r, xyz
end
# LLA or ECEF version
function get_range(sv::ReferenceView, lla_or_ecef::Union{LLA, Point3D}, ::ExtraOutput; face = sv.face, R = nothing)
ecef = if lla_or_ecef isa LLA
ECEFfromLLA(sv.ellipsoid)(lla_or_ecef)
else
lla_or_ecef
end
_R = isnothing(R) ? inv(sv.R * face_rotation(face)) : R
Δecef = ecef - sv.ecef
dist = norm(Δecef)
# Find if the target point is below the satellite, we do this by checking the last coordinate of the WND coordinates of the point
xyz = _R * Δecef
xyz[3] < 0 && return NaN, xyz
# We have to check that the given lla is visible from the satellite, this happens if either there is no intersection with earth in the direction of pointing, or if the first intersection happens for a range greater than th computed one
t₁, t₂ = _intersection_solutions(Δecef/dist, sv.ecef, sv.ellipsoid.a, sv.ellipsoid.b)
ref_t = t₁ > 0 ? t₁ : t₂ # t₁ always < t₂, so if positive that is the shortest distance, if negative t₂ might be also negative or positive, so we choose that
r = (
ref_t > 0 && # Check that the intersection is in the direction of pointing
ref_t < dist + 1e-5 # We aim for a precision of 1e-5 m. if ref_t is lower than the distance, it means that the target is behind earth
) ? NaN : dist
return r, xyz
end
# Single Output Version
get_range(rv, p; kwargs...) = get_range(rv, p, ExtraOutput(); kwargs...)[1]
# 2 ReferenceView version
function get_range(rv1::ReferenceView, rv2::ReferenceView, args...; kwargs...)
@assert rv1.earthmodel === rv2.earthmodel "When computing range between ReferenceView objects, the `EarthModel` they use internally must be the same"
get_range(rv1, rv2.ecef, args...; kwargs...)
end
end
# ╔═╡ 33e4c937-4443-4d97-a3ed-81479ede1e11
#=╠═╡
get_range(sv, SatView(LLA(0,0,700km), em); face = :NegativeZ)
╠═╡ =#
# ╔═╡ bd62bdd6-4de4-449c-b5f1-fb1b4f695cda
#=╠═╡
let
lla_ref = LLA(0,0,10km)
@benchmark get_range($sv, $lla_ref)
end
╠═╡ =#
# ╔═╡ 21bc362a-960c-4c5f-9118-7e1451edb996
#=╠═╡
let
lla_ref = LLA(0,0,10km)
ecef_ref = ECEFfromLLA(sv.ellipsoid)(lla_ref)
@benchmark get_range($sv, $ecef_ref)
end
╠═╡ =#
# ╔═╡ fd3a1dcb-2e64-47f8-a932-452742090ac1
#=╠═╡
let
sv = SatView(LLA(0,0,600km), em)
get_range(sv, (.1,0), ExtraOutput(); h = 100e3, face = :PositiveZ)
end
╠═╡ =#
# ╔═╡ 3381783a-b6d7-43c7-947d-c1de4437e25b
#=╠═╡
let
uv = UserView(LLA(0,0,0km), em)
get_range(uv, (0,0); h = 100e3)
end
╠═╡ =#
# ╔═╡ 449b49de-2951-41fc-ba46-89eaa6c52e79
# ╠═╡ skip_as_script = true
#=╠═╡
get_range(sv,LLA(0°,5°,10km),ExtraOutput())
╠═╡ =#
# ╔═╡ f99b984a-b7cf-4318-8f74-aacb644bc146
#=╠═╡
let
sv = SatView(LLA(0,0,600km), em)
get_range(sv, (0,0); face = :PositiveZ)
end
╠═╡ =#
# ╔═╡ b389a1db-c444-481c-b168-3069c7367276
#=╠═╡
let
sv = SatView(LLA(0,0,600km), em)
@benchmark get_range($sv,LLA(0,0,10km))
end
╠═╡ =#
# ╔═╡ 39a1850b-f64a-4157-8f07-d7a78918fea1
# ╠═╡ skip_as_script = true
#=╠═╡
md"""
## Get Pointing
"""
╠═╡ =#
# ╔═╡ 83634223-87d0-4c31-801a-af8a7f9f678a
begin
function _get_pointing(rv::ReferenceView, lla_or_ecef::Union{LLA, Point3D}, eo::ExtraOutput, pointing_type::PointingType; face = rv.face, R = nothing)
ecef = if lla_or_ecef isa LLA
ECEFfromLLA(rv.ellipsoid)(lla_or_ecef)
else
lla_or_ecef
end
# Check if earth is blocking
blocked, normalized_pdiff, pdiff = block_data = earth_blocking(rv.ecef,ecef, rv.ellipsoid.a,rv.ellipsoid.b, eo)
# If there is an intersection, we just return false
blocked && return SA_F64[NaN, NaN], SA_F64[NaN, NaN, NaN], block_data
_R = isnothing(R) ? inv(rv.R * face_rotation(face)) : R
xyz = _R * pdiff
x,y,z = normalize(xyz)
# If the target is behind, we return NaN as first output
z < 0 && return SA_F64[NaN, NaN], xyz, block_data
if pointing_type isa ThetaPhi
return SA_F64[acos(z), atan(y,x)], xyz, block_data
else
return SA_F64[x, y], xyz, block_data
end
end
"""
get_pointing(rv::ReferenceView, target::Union{LLA, Point3D, ReferenceView}[, ::ExtraOutput]; pointing_type::Symbol=:uv, face = rv.face, R=nothing)
Provide the 2-D angular pointing at which the target point (specified as LLA, ECEF or as another `ReferenceView`) is seen from the ReferenceView object `rv`.
`pointing_type` is used to select whether the output is returned in UV or
ThetaPhi [rad] coordinates. The following symbols are supported for this kwarg:
- `:thetaphi`, `:ThetaPhi` and `:θφ` can be used to represent pointing in ThetaPhi [rad]
- `:UV` and `:uv` can be used to represent pointing in UV
When called with an instance of `TelecomUtils.ExtraOutput` as last argument, the
function also returns the coordinated of the identified point in the local CRS
of `rv` as second output, and a NamedTuple containing information on the earth
blockage status as third argument.
For details on how to modify the reference pointing direction using the kwargs `face` and `R` look at the documentation of [`get_range`](@ref)
See also: [`get_mutual_pointing`](@ref), [`ReferenceView`](@ref),
[`get_range`](@ref), [`get_pointing`](@ref), [`get_lla`](@ref),
[`get_ecef`](@ref), [`get_distance_on_earth`](@ref), [`get_visibility`](@ref),
[`get_mutual_visibility`](@ref).
"""
function get_pointing(sv::ReferenceView, lla_or_ecef::Union{LLA, Point3D}, eo::ExtraOutput; pointing_type::Symbol=:uv, kwargs...)
_get_pointing(sv, lla_or_ecef, eo, PointingType(pointing_type); kwargs...)
end
# Single output version
get_pointing(sv, p; kwargs...) = get_pointing(sv, p, ExtraOutput(); kwargs...)[1]
# 2 ReferenceView version
function get_pointing(rv1::ReferenceView, rv2::ReferenceView, args...; kwargs...)
@assert rv1.earthmodel === rv2.earthmodel "When computing visibility angles between ReferenceView objects, the `EarthModel` they use internally must be the same"
get_pointing(rv1, rv2.ecef, args...; kwargs...)
end
end
# ╔═╡ d4ff395e-55bd-4240-89e3-a4be2dcf9ebe
#=╠═╡
let
lla_ref = LLA(1°, 0.5°, 0km)
ecef_ref = ECEFfromLLA(em.ellipsoid)(lla_ref)
@benchmark get_pointing($sv, $lla_ref)
# @code_warntype get_pointing(sv, lla_ref)
end
╠═╡ =#
# ╔═╡ cb4d8cd7-ff3e-43e6-ad38-8db343087214
#=╠═╡
let
lla_ref = LLA(1°, 0.5°, 0km)
ecef_ref = ECEFfromLLA(em.ellipsoid)(lla_ref)
@benchmark get_pointing($sv, $ecef_ref)
end
╠═╡ =#
# ╔═╡ f3b8b9ba-c428-4868-96b9-1983eabb3774
#=╠═╡
let
lla_ref = LLA(1°, 0.5°, 0km)
rv2 = SatView(lla_ref, em)
@benchmark get_pointing($sv, $rv2)
end
╠═╡ =#
# ╔═╡ 951fc2e2-ada9-4aad-876a-4cbd46200f6c
#=╠═╡
@bind test_face let
map(Symbol, collect(instances(Faces))) |> Select
end
╠═╡ =#
# ╔═╡ a93e2354-34d1-4d1d-ac3c-f4c99099820a
#=╠═╡
let
sat_lla = LLA(0°,0°,600km)
sv = SatView(sat_lla, em)
directions = (;
PositiveZ = LLA(0°,0°,500km),
PositiveX = LLA(0°,-1°,600km),
PositiveY = LLA(1°,0°,600km),
NegativeZ = LLA(0°,0°,700km),
NegativeX = LLA(0°,1°,600km),
NegativeY = LLA(-1°,0°,600km),
)
face = test_face
ref_uv = get_pointing(sv, getproperty(directions,face); face = face)
end
╠═╡ =#
# ╔═╡ 46eb3a4d-c80b-41a0-9333-aaf6411a010c
#=╠═╡
let
sv = SatView(LLA(0,0,600km),em)
lla = map(rand(1000), rand(1000)) do lat,lon
LLA(lat*°,lon*°,0km)
end
ecef = map(lla) do lla
ECEFfromLLA(sv.ellipsoid)(lla)
end
faces = rand([:PositiveZ, :PositiveY, :NegativeX], 1000)
map((x,y) -> get_pointing(sv,x; pointing_type= :thetaphi, face = y), ecef, faces)
end
╠═╡ =#
# ╔═╡ 314c4c2a-6f3b-4a7d-aa83-9073eb67897b
#=╠═╡
let
rv1 = SatView(LLA(0°, 10°, 800km), em)
lat = range(0°, 90°; step = 0.1°)
θ = map(lat) do lat
rv2 = SatView(LLA(lat, 10°, 600km), em; face = -3)
# v = get_visibility(rv1, rv2; fov = 60°)
p = get_pointing(rv1, rv2; pointing_type = :thetaphi)[1] |> rad2deg
end
x = to_degrees(lat)
plot(scatter(;x, y = θ))
# v,p
end
╠═╡ =#
# ╔═╡ a61e66cf-b430-4d84-a753-faf42f4b6337
md"""
## Get Mutual Pointing
"""
# ╔═╡ 98f3f83d-fcde-48ba-8855-c30643776d81
begin
"""
p₁, p₂ = get_mutual_pointing(rv1::ReferenceView, rv2::ReferenceView[, ::ExtraOutput]; pointing_type::Symbol=:uv, faces = (rv1.face, rv2.face), Rs=(nothing, nothing), short_circuit = true)
Provide the 2-D angular pointing in both directions between `rv1` and `rv2`:
- `p₁` is the pointing of `rv2` with respect to `rv1`
- `p₂` is the pointing of `rv1` with respect to `rv2`
`pointing_type` is used to select whether the outputs are returned in UV or
ThetaPhi [rad] coordinates. The following symbols are supported for this kwarg:
- `:thetaphi`, `:ThetaPhi` and `:θφ` can be used to represent pointing in ThetaPhi [rad]
- `:UV` and `:uv` can be used to represent pointing in UV
When called with an instance of `TelecomUtils.ExtraOutput` as last argument, the
function also returns the coordinated of the identified point in the local CRS
of `rv1` (or `rv2`). In this case:
- `p₁` is a tuple containing the pointing as first argument, the local CRS coordinates of `rv2` with respect to `rv1` as second argument and a NamedTuple with earth blockage data as third argument.
- `p₂` is a tuple containing the pointing as first argument, the local CRS coordinates of `rv1` with respect to `rv2` as second argument and a NamedTuple with earth blockage data as third argument.
`faces` and `Rs` are tuples containing the values of `face` and `R` for the two
ReferenceView objects. `faces[1]` is used as reference face for `rv1` while
`faces[2]` is used for `rv2`. Similarly for `Rs`.
The `short_circuit` kwarg defaults to true and permits to return early from the
function to save time. When `short_circuit == true`, if the fwd pointing is
false because of earth blockage, the function will directly return without
computing the pointing on the return direction.
For details on how to modify the reference pointing direction using `face` and
`R` look at the documentation of [`get_range`](@ref)
See also: [`ReferenceView`](@ref), [`get_range`](@ref), [`get_pointing`](@ref),
[`get_lla`](@ref), [`get_ecef`](@ref), [`get_distance_on_earth`](@ref),
[`get_visibility`](@ref), [`get_mutual_visibility`](@ref).
"""
function get_mutual_pointing(rv1::ReferenceView, rv2::ReferenceView, eo::ExtraOutput; pointing_type = :uv, faces = (rv1.face, rv2.face), Rs = (nothing, nothing), short_circuit = true)
# Pointing of rv2 as seen from rv1
p₁, xyz₁, block_data = data1 = get_pointing(rv1, rv2, eo; pointing_type, face = faces[1], R = Rs[1])
block_data.blocked && short_circuit && return data1, (SA_F64[NaN, NaN], SA_F64[NaN, NaN, NaN], block_data)
block_data = (;block_data..., pdiff = -block_data.pdiff, normalized_pdiff = -block_data.normalized_pdiff)
R = Rs[2]
_R = isnothing(R) ? inv(rv2.R * face_rotation(faces[2])) : R
xyz = _R * block_data.pdiff
x,y,z = normalize(xyz)
# If the target is behind, we return NaN as first output
z < 0 && return data1, (SA_F64[NaN, NaN], xyz, block_data)
data2 = if pointing_type isa ThetaPhi
SA_F64[acos(z), atan(y,x)], xyz, block_data
else
SA_F64[x, y], xyz, block_data
end
return data1, data2
end
get_mutual_pointing(rv1, rv2; kwargs...) = first.(get_mutual_pointing(rv1, rv2, ExtraOutput(); kwargs...))
end
# ╔═╡ cda9503a-5262-4318-9f5d-aea7910ab42e
#=╠═╡
let
sv1 = sv
sv2 = SatView(LLA(0°,1°, 1000km), em)
@benchmark get_mutual_pointing($sv1, $sv2; faces = (-3,3))
end
╠═╡ =#
# ╔═╡ 868bb972-aa0b-484d-b906-1628412b1646
#=╠═╡
let
sv1 = sv
sv2 = SatView(LLA(0°,1°, 1000km), em)
@benchmark get_mutual_pointing($sv1, $sv2, ExtraOutput();pointing_type = :thetaphi, faces = (-3,3))
end
╠═╡ =#
# ╔═╡ 83098fe4-2089-44be-8db9-6e2db74d886d
#=╠═╡
let
sv1 = sv
sv2 = SatView(LLA(0°,1°, 1000km), em)
get_mutual_pointing(sv1, sv2, ExtraOutput();pointing_type = :thetaphi, faces = (-3,3))
end
╠═╡ =#
# ╔═╡ 1de548c0-3830-4161-b80c-be72dd894e85
md"""
## Get Visibility
"""
# ╔═╡ 664be7e6-28f5-4ab3-b8f0-b7618bcd5fe5
# ╠═╡ skip_as_script = true
#=╠═╡
begin
@addmethod function earth_blocking(ecef1, ecef2, a, b, ::ExtraOutput)
@inline
# Check if the given ecef coordinate is visible from the satellite position or is obstructed from earth
pdiff = (ecef2 - ecef1)
# Find the magnitude of the difference to compare with the intersection solutions
t = norm(pdiff)
normalized_pdiff = pdiff ./ t
t₁,t₂ = _intersection_solutions(normalized_pdiff,ecef1,a,b)
# If there is an intersection, we return true and the normalized_pdiff
blocked = t₁ > 0 && t > t₁+1e-3
return (;blocked, normalized_pdiff, pdiff)
end
# Single Output Version
@addmethod earth_blocking(ecef1, ecef2, a, b) = earth_blocking(ecef1, ecef2, a, b, ExtraOutput())[1]
end
╠═╡ =#
# ╔═╡ 387d2c76-1a08-441c-97fc-7b0a90a95c9a
begin
"""
get_visibility(rv::ReferenceView, target::Union{LLA, Point3D, ReferenceView}[, ::ExtraOutput]; boresight = rv.face, fov = 90°)
Returns `true` if `target` is visible from `rv` assuming an antenna boresight
direction specified by `boresight` and a maximum Field of View from the
boresight specified by `fov`.
The `boresight` kwarg can be specified either as a `face` (compatible with the
input types specified in [`change_reference_face!`](@ref)) or as a 3D pointing
vector using either 3 elements Tuple, Vector or SVector. If a 3D vector is
provided, the function normalizes to unitary norm automatically.
The `fov` kwarg must be specified as an angle in radians if provided as a simple
Number. To provide a fov in degress, directly use ° from Unitful re-exported by
TelecomUtils.
When called with an instance of `TelecomUtils.ExtraOutput` as last argument, the
function also returns the pointing angle θ between `rv` and `target` as second
argument, and the direction in ECEF coordinates from `rv` to `target`.
See also: [`get_pointing`](@ref), [`get_mutual_pointing`](@ref),
[`ReferenceView`](@ref), [`get_range`](@ref), [`get_pointing`](@ref),
[`get_lla`](@ref), [`get_ecef`](@ref), [`get_distance_on_earth`](@ref),
[`get_mutual_visibility`](@ref).
"""
function get_visibility(rv::ReferenceView, lla_or_ecef::Union{LLA, Point3D}, eo::ExtraOutput; boresight = rv.face, fov = 90°)
ecef = if lla_or_ecef isa LLA
ECEFfromLLA(rv.ellipsoid)(lla_or_ecef)
else
lla_or_ecef
end
blocked, normalized_pdiff, _ = block_data = earth_blocking(rv.ecef,ecef, rv.ellipsoid.a,rv.ellipsoid.b, eo)
# # If there is an intersection, we just return false
blocked && return false, NaN, normalized_pdiff
xyz = rv.R' * normalized_pdiff
θ = acos(dot(xyz, boresight_versor(boresight)))
return θ <= fov, θ, normalized_pdiff
end
# Double RefView method
function get_visibility(rv1::ReferenceView, rv2::ReferenceView, args...; kwargs...)
@assert rv1.earthmodel === rv2.earthmodel "When computing visibility between ReferenceView objects, the `EarthModel` they use internally must be the same"
get_visibility(rv1, rv2.ecef, args...; kwargs...)
end
# Single Output
get_visibility(rv1::ReferenceView, target; kwargs...) = get_visibility(rv1, target, ExtraOutput(); kwargs...)[1]
end
# ╔═╡ abb3a7ee-ff60-421c-8a85-23c3ef1ce6af
#=╠═╡
let
rv1 = SatView(LLA(0°, 10°, 800km), em)
rv2 = SatView(LLA(35°, 10°, 600km), em; face = -3)
v = get_visibility(rv1, rv2; fov = 60°)
p = get_pointing(rv1, rv2, ExtraOutput(); pointing_type = :uv)
v,p
end
╠═╡ =#
# ╔═╡ b75df549-ed46-4def-9fe0-e63ac496b3f4
#=╠═╡
let
rv1 = SatView(LLA(0°, 10°, 800km), em)
rv2 = LLA(0°, 190°, 600km)
v = get_visibility(rv1, rv2; fov = 60°)
p = get_pointing(rv1, rv2, ExtraOutput(); pointing_type = :uv)
v,p
end
╠═╡ =#
# ╔═╡ 864cc72a-6252-45a6-8ed4-81d960c0b080
#=╠═╡
let
rv1 = SatView(LLA(1°, 10°, 800km), em)
rv2 = SatView(LLA(0°, 10°, 600km), em; face = -3)
v = get_visibility(rv1, rv2; fov = 60°)
p = get_pointing(rv1, rv2; pointing_type = :thetaphi) |> x -> rad2deg.(x)
v,p
end
╠═╡ =#
# ╔═╡ ecd8d6ac-27f6-4a50-81ff-9c8f0b7fed47
#=╠═╡
let
rv1 = SatView(LLA(50°, 10°, 800km), em)
rv2 = SatView(LLA(0°, 10°, 600km), em; face = -3)
@benchmark get_visibility($rv1, $rv2; fov = 60°)
end
╠═╡ =#
# ╔═╡ fe0680dd-95fb-40a0-8063-97c009f552dd
#=╠═╡