-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpn.go
10273 lines (8525 loc) · 364 KB
/
pn.go
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
// Copyright 2022 HE Boliang
// All rights reserved.
package gofa
// SOFA Precession / Nutation / Polar Motion
/*
Bi00 Frame bias components, IAU 2000 Frame bias components of IAU 2000
precession-nutation models; part of the Mathews-Herring-Buffett (MHB2000)
nutation series, with additions.
Returned:
dpsibi,depsbi float64 longitude and obliquity corrections
dra float64 the ICRS RA of the J2000.0 mean equinox
Notes:
1) The frame bias corrections in longitude and obliquity (radians)
are required in order to correct for the offset between the GCRS
pole and the mean J2000.0 pole. They define, with respect to the
GCRS frame, a J2000.0 mean pole that is consistent with the rest
of the IAU 2000A precession-nutation model.
2) In addition to the displacement of the pole, the complete
description of the frame bias requires also an offset in right
ascension. This is not part of the IAU 2000A model, and is from
Chapront et al. (2002). It is returned in radians.
3) This is a supplemented implementation of one aspect of the IAU
2000A nutation model, formally adopted by the IAU General
Assembly in 2000, namely MHB2000 (Mathews et al. 2002).
References:
Chapront, J., Chapront-Touze, M. & Francou, G., Astron.
Astrophys., 387, 700, 2002.
Mathews, P.M., Herring, T.A., Buffet, B.A., "Modeling of nutation
and precession: New nutation series for nonrigid Earth and
insights into the Earth's interior", J.Geophys.Res., 107, B4,
2002. The MHB2000 code itself was obtained on 2002 September 9
from ftp://maia.usno.navy.mil/conv2000/chapter5/IAU2000A.
*/
func Bi00(dpsibi, depsbi, dra *float64) {
/* The frame bias corrections in longitude and obliquity */
const DPBIAS = -0.041775 * DAS2R
const DEBIAS = -0.0068192 * DAS2R
/* The ICRS RA of the J2000.0 equinox (Chapront et al., 2002) */
const DRA0 = -0.0146 * DAS2R
/* Return the results (which are fixed). */
*dpsibi = DPBIAS
*depsbi = DEBIAS
*dra = DRA0
}
/*
Bp00 Frame bias and precession matrices, IAU 2000 Frame bias and precession, IAU 2000.
Given:
date1,date2 float64 TT as a 2-part Julian Date (Note 1)
Returned:
rb [3][3]float64 frame bias matrix (Note 2)
rp [3][3]float64 precession matrix (Note 3)
rbp [3][3]float64 bias-precession matrix (Note 4)
Notes:
1) The TT date date1+date2 is a Julian Date, apportioned in any
convenient way between the two arguments. For example,
JD(TT)=2450123.7 could be expressed in any of these ways,
among others:
date1 date2
2450123.7 0.0 (JD method)
2451545.0 -1421.3 (J2000 method)
2400000.5 50123.2 (MJD method)
2450123.5 0.2 (date & time method)
The JD method is the most natural and convenient to use in
cases where the loss of several decimal digits of resolution
is acceptable. The J2000 method is best matched to the way
the argument is handled internally and will deliver the
optimum resolution. The MJD method and the date & time methods
are both good compromises between resolution and convenience.
2) The matrix rb transforms vectors from GCRS to mean J2000.0 by
applying frame bias.
3) The matrix rp transforms vectors from J2000.0 mean equator and
equinox to mean equator and equinox of date by applying
precession.
4) The matrix rbp transforms vectors from GCRS to mean equator and
equinox of date by applying frame bias then precession. It is
the product rp x rb.
5) It is permissible to re-use the same array in the returned
arguments. The arrays are filled in the order given.
Called:
Bi00 frame bias components, IAU 2000
Pr00 IAU 2000 precession adjustments
Ir initialize r-matrix to identity
Rx rotate around X-axis
Ry rotate around Y-axis
Rz rotate around Z-axis
Cr copy r-matrix
Rxr product of two r-matrices
Reference:
"Expressions for the Celestial Intermediate Pole and Celestial
Ephemeris Origin consistent with the IAU 2000A precession-
nutation model", Astron.Astrophys. 400, 1145-1154 (2003)
n.b. The celestial ephemeris origin (CEO) was renamed "celestial
intermediate origin" (CIO) by IAU 2006 Resolution 2.
*/
func Bp00(date1, date2 float64, rb, rp, rbp *[3][3]float64) {
/* J2000.0 obliquity (Lieske et al. 1977) */
const EPS0 = 84381.448 * DAS2R
var t, dpsibi, depsbi, dra0, psia77, oma77, chia,
dpsipr, depspr, psia, oma float64
var rbw [3][3]float64
/* Interval between fundamental epoch J2000.0 and current date (JC). */
t = ((date1 - DJ00) + date2) / DJC
/* Frame bias. */
Bi00(&dpsibi, &depsbi, &dra0)
/* Precession angles (Lieske et al. 1977) */
psia77 = (5038.7784 + (-1.07259+(-0.001147)*t)*t) * t * DAS2R
oma77 = EPS0 + ((0.05127+(-0.007726)*t)*t)*t*DAS2R
chia = (10.5526 + (-2.38064+(-0.001125)*t)*t) * t * DAS2R
/* Apply IAU 2000 precession corrections. */
Pr00(date1, date2, &dpsipr, &depspr)
psia = psia77 + dpsipr
oma = oma77 + depspr
/* Frame bias matrix: GCRS to J2000.0. */
Ir(&rbw)
Rz(dra0, &rbw)
Ry(dpsibi*sin(EPS0), &rbw)
Rx(-depsbi, &rbw)
Cr(rbw, rb)
/* Precession matrix: J2000.0 to mean of date. */
Ir(rp)
Rx(EPS0, rp)
Rz(-psia, rp)
Rx(-oma, rp)
Rz(chia, rp)
/* Bias-precession matrix: GCRS to mean of date. */
Rxr(*rp, rbw, rbp)
}
/*
Bp06 Frame bias and precession matrices, IAU 2006 Frame bias and precession, IAU 2006.
Given:
date1,date2 float64 TT as a 2-part Julian Date (Note 1)
Returned:
rb [3][3]float64 frame bias matrix (Note 2)
rp [3][3]float64 precession matrix (Note 3)
rbp [3][3]float64 bias-precession matrix (Note 4)
Notes:
1) The TT date date1+date2 is a Julian Date, apportioned in any
convenient way between the two arguments. For example,
JD(TT)=2450123.7 could be expressed in any of these ways,
among others:
date1 date2
2450123.7 0.0 (JD method)
2451545.0 -1421.3 (J2000 method)
2400000.5 50123.2 (MJD method)
2450123.5 0.2 (date & time method)
The JD method is the most natural and convenient to use in
cases where the loss of several decimal digits of resolution
is acceptable. The J2000 method is best matched to the way
the argument is handled internally and will deliver the
optimum resolution. The MJD method and the date & time methods
are both good compromises between resolution and convenience.
2) The matrix rb transforms vectors from GCRS to mean J2000.0 by
applying frame bias.
3) The matrix rp transforms vectors from mean J2000.0 to mean of
date by applying precession.
4) The matrix rbp transforms vectors from GCRS to mean of date by
applying frame bias then precession. It is the product rp x rb.
5) It is permissible to re-use the same array in the returned
arguments. The arrays are filled in the order given.
Called:
Pfw06 bias-precession F-W angles, IAU 2006
Fw2m F-W angles to r-matrix
Pmat06 PB matrix, IAU 2006
Tr transpose r-matrix
Rxr product of two r-matrices
Cr copy r-matrix
References:
Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981
*/
func Bp06(date1, date2 float64, rb, rp, rbp *[3][3]float64) {
var gamb, phib, psib, epsa float64
var rbpw, rbt [3][3]float64
/* B matrix. */
Pfw06(DJM0, DJM00, &gamb, &phib, &psib, &epsa)
Fw2m(gamb, phib, psib, epsa, rb)
/* PxB matrix (temporary). */
Pmat06(date1, date2, &rbpw)
/* P matrix. */
Tr(*rb, &rbt)
Rxr(rbpw, rbt, rp)
/* PxB matrix. */
Cr(rbpw, rbp)
}
/*
Bpn2xy CIP XY given Bias-precession-nutation matrix Extract from the
bias-precession-nutation matrix the X,Y coordinates of the Celestial Intermediate Pole.
Given:
rbpn [3][3]float64 celestial-to-true matrix (Note 1)
Returned:
x,y float64 Celestial Intermediate Pole (Note 2)
Notes:
1) The matrix rbpn transforms vectors from GCRS to true equator (and
CIO or equinox) of date, and therefore the Celestial Intermediate
Pole unit vector is the bottom row of the matrix.
2) The arguments x,y are components of the Celestial Intermediate
Pole unit vector in the Geocentric Celestial Reference System.
Reference:
"Expressions for the Celestial Intermediate Pole and Celestial
Ephemeris Origin consistent with the IAU 2000A precession-
nutation model", Astron.Astrophys. 400, 1145-1154
(2003)
n.b. The celestial ephemeris origin (CEO) was renamed "celestial
intermediate origin" (CIO) by IAU 2006 Resolution 2.
*/
func Bpn2xy(rbpn [3][3]float64, x, y *float64) {
/* Extract the X,Y coordinates. */
*x = rbpn[2][0]
*y = rbpn[2][1]
}
/*
C2i00a Celestial-to-intermediate matrix, IAU 2000A Form the
celestial-to-intermediate matrix for a given date using the
IAU 2000A precession-nutation model.
Given:
date1,date2 float64 TT as a 2-part Julian Date (Note 1)
Returned:
rc2i [3][3]float64 celestial-to-intermediate matrix (Note 2)
Notes:
1) The TT date date1+date2 is a Julian Date, apportioned in any
convenient way between the two arguments. For example,
JD(TT)=2450123.7 could be expressed in any of these ways,
among others:
date1 date2
2450123.7 0.0 (JD method)
2451545.0 -1421.3 (J2000 method)
2400000.5 50123.2 (MJD method)
2450123.5 0.2 (date & time method)
The JD method is the most natural and convenient to use in
cases where the loss of several decimal digits of resolution
is acceptable. The J2000 method is best matched to the way
the argument is handled internally and will deliver the
optimum resolution. The MJD method and the date & time methods
are both good compromises between resolution and convenience.
2) The matrix rc2i is the first stage in the transformation from
celestial to terrestrial coordinates:
[TRS] = RPOM * R_3(ERA) * rc2i * [CRS]
= rc2t * [CRS]
where [CRS] is a vector in the Geocentric Celestial Reference
System and [TRS] is a vector in the International Terrestrial
Reference System (see IERS Conventions 2003), ERA is the Earth
Rotation Angle and RPOM is the polar motion matrix.
3) A faster, but slightly less accurate, result (about 1 mas) can be
obtained by using instead the iauC2i00b function.
Called:
Pnm00a classical NPB matrix, IAU 2000A
C2ibpn celestial-to-intermediate matrix, given NPB matrix
References:
"Expressions for the Celestial Intermediate Pole and Celestial
Ephemeris Origin consistent with the IAU 2000A precession-
nutation model", Astron.Astrophys. 400, 1145-1154
(2003)
n.b. The celestial ephemeris origin (CEO) was renamed "celestial
intermediate origin" (CIO) by IAU 2006 Resolution 2.
McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
IERS Technical Note No. 32, BKG (2004)
*/
func C2i00a(date1, date2 float64, rc2i *[3][3]float64) {
var rbpn [3][3]float64
/* Obtain the celestial-to-true matrix (IAU 2000A). */
Pnm00a(date1, date2, &rbpn)
/* Form the celestial-to-intermediate matrix. */
C2ibpn(date1, date2, rbpn, rc2i)
}
/*
C2i00b Celestial-to-intermediate matrix, IAU 2000B
Form the celestial-to-intermediate matrix for a given date using the
IAU 2000B precession-nutation model.
Given:
date1,date2 float64 TT as a 2-part Julian Date (Note 1)
Returned:
rc2i [3][3]float64 celestial-to-intermediate matrix (Note 2)
Notes:
1) The TT date date1+date2 is a Julian Date, apportioned in any
convenient way between the two arguments. For example,
JD(TT)=2450123.7 could be expressed in any of these ways,
among others:
date1 date2
2450123.7 0.0 (JD method)
2451545.0 -1421.3 (J2000 method)
2400000.5 50123.2 (MJD method)
2450123.5 0.2 (date & time method)
The JD method is the most natural and convenient to use in
cases where the loss of several decimal digits of resolution
is acceptable. The J2000 method is best matched to the way
the argument is handled internally and will deliver the
optimum resolution. The MJD method and the date & time methods
are both good compromises between resolution and convenience.
2) The matrix rc2i is the first stage in the transformation from
celestial to terrestrial coordinates:
[TRS] = RPOM * R_3(ERA) * rc2i * [CRS]
= rc2t * [CRS]
where [CRS] is a vector in the Geocentric Celestial Reference
System and [TRS] is a vector in the International Terrestrial
Reference System (see IERS Conventions 2003), ERA is the Earth
Rotation Angle and RPOM is the polar motion matrix.
3) The present function is faster, but slightly less accurate (about
1 mas), than the iauC2i00a function.
Called:
Pnm00b classical NPB matrix, IAU 2000B
C2ibpn celestial-to-intermediate matrix, given NPB matrix
References:
"Expressions for the Celestial Intermediate Pole and Celestial
Ephemeris Origin consistent with the IAU 2000A precession-
nutation model", Astron.Astrophys. 400, 1145-1154
(2003)
n.b. The celestial ephemeris origin (CEO) was renamed "celestial
intermediate origin" (CIO) by IAU 2006 Resolution 2.
McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
IERS Technical Note No. 32, BKG (2004)
*/
func C2i00b(date1, date2 float64, rc2i *[3][3]float64) {
var rbpn [3][3]float64
/* Obtain the celestial-to-true matrix (IAU 2000B). */
Pnm00b(date1, date2, &rbpn)
/* Form the celestial-to-intermediate matrix. */
C2ibpn(date1, date2, rbpn, rc2i)
}
/*
C2i06a Celestial-to-intermediate matrix, IAU 2006/2000A
Form the celestial-to-intermediate matrix for a given date using the
IAU 2006 precession and IAU 2000A nutation models.
Given:
date1,date2 float64 TT as a 2-part Julian Date (Note 1)
Returned:
rc2i [3][3]float64 celestial-to-intermediate matrix (Note 2)
Notes:
1) The TT date date1+date2 is a Julian Date, apportioned in any
convenient way between the two arguments. For example,
JD(TT)=2450123.7 could be expressed in any of these ways,
among others:
date1 date2
2450123.7 0.0 (JD method)
2451545.0 -1421.3 (J2000 method)
2400000.5 50123.2 (MJD method)
2450123.5 0.2 (date & time method)
The JD method is the most natural and convenient to use in
cases where the loss of several decimal digits of resolution
is acceptable. The J2000 method is best matched to the way
the argument is handled internally and will deliver the
optimum resolution. The MJD method and the date & time methods
are both good compromises between resolution and convenience.
2) The matrix rc2i is the first stage in the transformation from
celestial to terrestrial coordinates:
[TRS] = RPOM * R_3(ERA) * rc2i * [CRS]
= RC2T * [CRS]
where [CRS] is a vector in the Geocentric Celestial Reference
System and [TRS] is a vector in the International Terrestrial
Reference System (see IERS Conventions 2003), ERA is the Earth
Rotation Angle and RPOM is the polar motion matrix.
Called:
Pnm06a classical NPB matrix, IAU 2006/2000A
Bpn2xy extract CIP X,Y coordinates from NPB matrix
S06 the CIO locator s, given X,Y, IAU 2006
C2ixys celestial-to-intermediate matrix, given X,Y and s
References:
McCarthy, D. D., Petit, G. (eds.), 2004, IERS Conventions (2003),
IERS Technical Note No. 32, BKG
*/
func C2i06a(date1, date2 float64, rc2i *[3][3]float64) {
var rbpn [3][3]float64
var x, y, s float64
/* Obtain the celestial-to-true matrix (IAU 2006/2000A). */
Pnm06a(date1, date2, &rbpn)
/* Extract the X,Y coordinates. */
Bpn2xy(rbpn, &x, &y)
/* Obtain the CIO locator. */
s = S06(date1, date2, x, y)
/* Form the celestial-to-intermediate matrix. */
C2ixys(x, y, s, rc2i)
}
/*
C2ibpn Celestial-to-intermediate matrix given B-P-N
Form the celestial-to-intermediate matrix for a given date given
the bias-precession-nutation matrix. IAU 2000.
Given:
date1,date2 float64 TT as a 2-part Julian Date (Note 1)
rbpn [3][3]float64 celestial-to-true matrix (Note 2)
Returned:
rc2i [3][3]float64 celestial-to-intermediate matrix (Note 3)
Notes:
1) The TT date date1+date2 is a Julian Date, apportioned in any
convenient way between the two arguments. For example,
JD(TT)=2450123.7 could be expressed in any of these ways,
among others:
date1 date2
2450123.7 0.0 (JD method)
2451545.0 -1421.3 (J2000 method)
2400000.5 50123.2 (MJD method)
2450123.5 0.2 (date & time method)
The JD method is the most natural and convenient to use in
cases where the loss of several decimal digits of resolution
is acceptable. The J2000 method is best matched to the way
the argument is handled internally and will deliver the
optimum resolution. The MJD method and the date & time methods
are both good compromises between resolution and convenience.
2) The matrix rbpn transforms vectors from GCRS to true equator (and
CIO or equinox) of date. Only the CIP (bottom row) is used.
3) The matrix rc2i is the first stage in the transformation from
celestial to terrestrial coordinates:
[TRS] = RPOM * R_3(ERA) * rc2i * [CRS]
= RC2T * [CRS]
where [CRS] is a vector in the Geocentric Celestial Reference
System and [TRS] is a vector in the International Terrestrial
Reference System (see IERS Conventions 2003), ERA is the Earth
Rotation Angle and RPOM is the polar motion matrix.
4) Although its name does not include "00", This function is in fact
specific to the IAU 2000 models.
Called:
Bpn2xy extract CIP X,Y coordinates from NPB matrix
C2ixy celestial-to-intermediate matrix, given X,Y
References:
"Expressions for the Celestial Intermediate Pole and Celestial
Ephemeris Origin consistent with the IAU 2000A precession-
nutation model", Astron.Astrophys. 400, 1145-1154 (2003)
n.b. The celestial ephemeris origin (CEO) was renamed "celestial
intermediate origin" (CIO) by IAU 2006 Resolution 2.
McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
IERS Technical Note No. 32, BKG (2004)
*/
func C2ibpn(date1, date2 float64, rbpn [3][3]float64, rc2i *[3][3]float64) {
var x, y float64
/* Extract the X,Y coordinates. */
Bpn2xy(rbpn, &x, &y)
/* Form the celestial-to-intermediate matrix (n.b. IAU 2000 specific). */
C2ixy(date1, date2, x, y, rc2i)
}
/*
C2ixy Celestial-to-intermediate matrix given CIP
Form the celestial to intermediate-frame-of-date matrix for a given
date when the CIP X,Y coordinates are known. IAU 2000.
Given:
date1,date2 float64 TT as a 2-part Julian Date (Note 1)
x,y float64 Celestial Intermediate Pole (Note 2)
Returned:
rc2i [3][3]float64 celestial-to-intermediate matrix (Note 3)
Notes:
1) The TT date date1+date2 is a Julian Date, apportioned in any
convenient way between the two arguments. For example,
JD(TT)=2450123.7 could be expressed in any of these ways,
among others:
date1 date2
2450123.7 0.0 (JD method)
2451545.0 -1421.3 (J2000 method)
2400000.5 50123.2 (MJD method)
2450123.5 0.2 (date & time method)
The JD method is the most natural and convenient to use in
cases where the loss of several decimal digits of resolution
is acceptable. The J2000 method is best matched to the way
the argument is handled internally and will deliver the
optimum resolution. The MJD method and the date & time methods
are both good compromises between resolution and convenience.
2) The Celestial Intermediate Pole coordinates are the x,y components
of the unit vector in the Geocentric Celestial Reference System.
3) The matrix rc2i is the first stage in the transformation from
celestial to terrestrial coordinates:
[TRS] = RPOM * R_3(ERA) * rc2i * [CRS]
= RC2T * [CRS]
where [CRS] is a vector in the Geocentric Celestial Reference
System and [TRS] is a vector in the International Terrestrial
Reference System (see IERS Conventions 2003), ERA is the Earth
Rotation Angle and RPOM is the polar motion matrix.
4) Although its name does not include "00", This function is in fact
specific to the IAU 2000 models.
Called:
C2ixys celestial-to-intermediate matrix, given X,Y and s
S00 the CIO locator s, given X,Y, IAU 2000A
Reference:
McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
IERS Technical Note No. 32, BKG (2004)
*/
func C2ixy(date1, date2 float64, x, y float64, rc2i *[3][3]float64) {
/* Compute s and then the matrix. */
C2ixys(x, y, S00(date1, date2, x, y), rc2i)
}
/*
C2ixys Celestial-to-intermediate matrix given CIP and s
Form the celestial to intermediate-frame-of-date matrix given the CIP
X,Y and the CIO locator s.
Given:
x,y float64 Celestial Intermediate Pole (Note 1)
s float64 the CIO locator s (Note 2)
Returned:
rc2i [3][3]float64 celestial-to-intermediate matrix (Note 3)
Notes:
1) The Celestial Intermediate Pole coordinates are the x,y
components of the unit vector in the Geocentric Celestial
Reference System.
2) The CIO locator s (in radians) positions the Celestial
Intermediate Origin on the equator of the CIP.
3) The matrix rc2i is the first stage in the transformation from
celestial to terrestrial coordinates:
[TRS] = RPOM * R_3(ERA) * rc2i * [CRS]
= RC2T * [CRS]
where [CRS] is a vector in the Geocentric Celestial Reference
System and [TRS] is a vector in the International Terrestrial
Reference System (see IERS Conventions 2003), ERA is the Earth
Rotation Angle and RPOM is the polar motion matrix.
Called:
Ir initialize r-matrix to identity
Rz rotate around Z-axis
Ry rotate around Y-axis
Reference:
McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
IERS Technical Note No. 32, BKG (2004)
*/
func C2ixys(x, y, s float64, rc2i *[3][3]float64) {
var r2, e, d float64
/* Obtain the spherical angles E and d. */
r2 = x*x + y*y
// e = (r2 > 0.0) ? atan2(y, x) : 0.0;
if r2 > 0.0 {
e = atan2(y, x)
} else {
e = 0.0
}
d = atan(sqrt(r2 / (1.0 - r2)))
/* Form the matrix. */
Ir(rc2i)
Rz(e, rc2i)
Ry(d, rc2i)
Rz(-(e + s), rc2i)
}
/*
C2t00a Celestial-to-terrestrial matrix, IAU 2000A
Form the celestial to terrestrial matrix given the date, the UT1 and
the polar motion, using the IAU 2000A precession-nutation model.
Given:
tta,ttb float64 TT as a 2-part Julian Date (Note 1)
uta,utb float64 UT1 as a 2-part Julian Date (Note 1)
xp,yp float64 CIP coordinates (radians, Note 2)
Returned:
rc2t [3][3]float64 celestial-to-terrestrial matrix (Note 3)
Notes:
1) The TT and UT1 dates tta+ttb and uta+utb are Julian Dates,
apportioned in any convenient way between the arguments uta and
utb. For example, JD(UT1)=2450123.7 could be expressed in any of
these ways, among others:
uta utb
2450123.7 0.0 (JD method)
2451545.0 -1421.3 (J2000 method)
2400000.5 50123.2 (MJD method)
2450123.5 0.2 (date & time method)
The JD method is the most natural and convenient to use in
cases where the loss of several decimal digits of resolution is
acceptable. The J2000 and MJD methods are good compromises
between resolution and convenience. In the case of uta,utb, the
date & time method is best matched to the Earth rotation angle
algorithm used: maximum precision is delivered when the uta
argument is for 0hrs UT1 on the day in question and the utb
argument lies in the range 0 to 1, or vice versa.
2) The arguments xp and yp are the coordinates (in radians) of the
Celestial Intermediate Pole with respect to the International
Terrestrial Reference System (see IERS Conventions 2003),
measured along the meridians 0 and 90 deg west respectively.
3) The matrix rc2t transforms from celestial to terrestrial
coordinates:
[TRS] = RPOM * R_3(ERA) * RC2I * [CRS]
= rc2t * [CRS]
where [CRS] is a vector in the Geocentric Celestial Reference
System and [TRS] is a vector in the International Terrestrial
Reference System (see IERS Conventions 2003), RC2I is the
celestial-to-intermediate matrix, ERA is the Earth rotation
angle and RPOM is the polar motion matrix.
4) A faster, but slightly less accurate, result (about 1 mas) can
be obtained by using instead the iauC2t00b function.
Called:
C2i00a celestial-to-intermediate matrix, IAU 2000A
Era00 Earth rotation angle, IAU 2000
Sp00 the TIO locator s', IERS 2000
Pom00 polar motion matrix
C2tcio form CIO-based celestial-to-terrestrial matrix
Reference:
McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
IERS Technical Note No. 32, BKG (2004)
*/
func C2t00a(tta, ttb float64, uta, utb float64, xp, yp float64, rc2t *[3][3]float64) {
var rc2i, rpom [3][3]float64
var era, sp float64
/* Form the celestial-to-intermediate matrix for this TT (IAU 2000A). */
C2i00a(tta, ttb, &rc2i)
/* Predict the Earth rotation angle for this UT1. */
era = Era00(uta, utb)
/* Estimate s'. */
sp = Sp00(tta, ttb)
/* Form the polar motion matrix. */
Pom00(xp, yp, sp, &rpom)
/* Combine to form the celestial-to-terrestrial matrix. */
C2tcio(rc2i, era, rpom, rc2t)
}
/*
C2t00b Celestial-to-terrestrial matrix, IAU 2000B
Form the celestial to terrestrial matrix given the date, the UT1 and
the polar motion, using the IAU 2000B precession-nutation model.
Given:
tta,ttb float64 TT as a 2-part Julian Date (Note 1)
uta,utb float64 UT1 as a 2-part Julian Date (Note 1)
xp,yp float64 coordinates of the pole (radians, Note 2)
Returned:
rc2t [3][3]float64 celestial-to-terrestrial matrix (Note 3)
Notes:
1) The TT and UT1 dates tta+ttb and uta+utb are Julian Dates,
apportioned in any convenient way between the arguments uta and
utb. For example, JD(UT1)=2450123.7 could be expressed in any of
these ways, among others:
uta utb
2450123.7 0.0 (JD method)
2451545.0 -1421.3 (J2000 method)
2400000.5 50123.2 (MJD method)
2450123.5 0.2 (date & time method)
The JD method is the most natural and convenient to use in
cases where the loss of several decimal digits of resolution is
acceptable. The J2000 and MJD methods are good compromises
between resolution and convenience. In the case of uta,utb, the
date & time method is best matched to the Earth rotation angle
algorithm used: maximum precision is delivered when the uta
argument is for 0hrs UT1 on the day in question and the utb
argument lies in the range 0 to 1, or vice versa.
2) The arguments xp and yp are the coordinates (in radians) of the
Celestial Intermediate Pole with respect to the International
Terrestrial Reference System (see IERS Conventions 2003),
measured along the meridians 0 and 90 deg west respectively.
3) The matrix rc2t transforms from celestial to terrestrial
coordinates:
[TRS] = RPOM * R_3(ERA) * RC2I * [CRS]
= rc2t * [CRS]
where [CRS] is a vector in the Geocentric Celestial Reference
System and [TRS] is a vector in the International Terrestrial
Reference System (see IERS Conventions 2003), RC2I is the
celestial-to-intermediate matrix, ERA is the Earth rotation
angle and RPOM is the polar motion matrix.
4) The present function is faster, but slightly less accurate (about
1 mas), than the iauC2t00a function.
Called:
C2i00b celestial-to-intermediate matrix, IAU 2000B
Era00 Earth rotation angle, IAU 2000
Pom00 polar motion matrix
C2tcio form CIO-based celestial-to-terrestrial matrix
Reference:
McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
IERS Technical Note No. 32, BKG (2004)
*/
func C2t00b(tta, ttb float64, uta, utb float64, xp, yp float64, rc2t *[3][3]float64) {
var rc2i, rpom [3][3]float64
var era float64
/* Form the celestial-to-intermediate matrix for this TT (IAU 2000B). */
C2i00b(tta, ttb, &rc2i)
/* Predict the Earth rotation angle for this UT1. */
era = Era00(uta, utb)
/* Form the polar motion matrix (neglecting s'). */
Pom00(xp, yp, 0.0, &rpom)
/* Combine to form the celestial-to-terrestrial matrix. */
C2tcio(rc2i, era, rpom, rc2t)
}
/*
C2t06a Celestial-to-terrestrial matrix, IAU 2006/2000A
Form the celestial to terrestrial matrix given the date, the UT1 and
the polar motion, using the IAU 2006/2000A precession-nutation
model.
Given:
tta,ttb float64 TT as a 2-part Julian Date (Note 1)
uta,utb float64 UT1 as a 2-part Julian Date (Note 1)
xp,yp float64 coordinates of the pole (radians, Note 2)
Returned:
rc2t [3][3]float64 celestial-to-terrestrial matrix (Note 3)
Notes:
1) The TT and UT1 dates tta+ttb and uta+utb are Julian Dates,
apportioned in any convenient way between the two arguments. For
example, JD(UT1)=2450123.7 could be expressed in any of
these ways, among others:
uta utb
2450123.7 0.0 (JD method)
2451545.0 -1421.3 (J2000 method)
2400000.5 50123.2 (MJD method)
2450123.5 0.2 (date & time method)
The JD method is the most natural and convenient to use in
cases where the loss of several decimal digits of resolution is
acceptable. The J2000 and MJD methods are good compromises
between resolution and convenience. In the case of uta,utb, the
date & time method is best matched to the Earth rotation angle
algorithm used: maximum precision is delivered when the uta
argument is for 0hrs UT1 on the day in question and the utb
argument lies in the range 0 to 1, or vice versa.
2) The arguments xp and yp are the coordinates (in radians) of the
Celestial Intermediate Pole with respect to the International
Terrestrial Reference System (see IERS Conventions 2003),
measured along the meridians 0 and 90 deg west respectively.
3) The matrix rc2t transforms from celestial to terrestrial
coordinates:
[TRS] = RPOM * R_3(ERA) * RC2I * [CRS]
= rc2t * [CRS]
where [CRS] is a vector in the Geocentric Celestial Reference
System and [TRS] is a vector in the International Terrestrial
Reference System (see IERS Conventions 2003), RC2I is the
celestial-to-intermediate matrix, ERA is the Earth rotation
angle and RPOM is the polar motion matrix.
Called:
C2i06a celestial-to-intermediate matrix, IAU 2006/2000A
Era00 Earth rotation angle, IAU 2000
Sp00 the TIO locator s', IERS 2000
Pom00 polar motion matrix
C2tcio form CIO-based celestial-to-terrestrial matrix
Reference:
McCarthy, D. D., Petit, G. (eds.), 2004, IERS Conventions (2003),
IERS Technical Note No. 32, BKG
*/
func C2t06a(tta, ttb float64, uta, utb float64, xp, yp float64, rc2t *[3][3]float64) {
var rc2i, rpom [3][3]float64
var era, sp float64
/* Form the celestial-to-intermediate matrix for this TT. */
C2i06a(tta, ttb, &rc2i)
/* Predict the Earth rotation angle for this UT1. */
era = Era00(uta, utb)
/* Estimate s'. */
sp = Sp00(tta, ttb)
/* Form the polar motion matrix. */
Pom00(xp, yp, sp, &rpom)
/* Combine to form the celestial-to-terrestrial matrix. */
C2tcio(rc2i, era, rpom, rc2t)
}
/*
C2tcio Form CIO-based Celestial-to-terrestrial matrix
Assemble the celestial to terrestrial matrix from CIO-based
components (the celestial-to-intermediate matrix, the Earth Rotation
Angle and the polar motion matrix).
Given:
rc2i [3][3]float64 celestial-to-intermediate matrix
era float64 Earth rotation angle (radians)
rpom [3][3]float64 polar-motion matrix
Returned:
rc2t [3][3]float64 celestial-to-terrestrial matrix
Notes:
1) This function constructs the rotation matrix that transforms
vectors in the celestial system into vectors in the terrestrial
system. It does so starting from precomputed components, namely
the matrix which rotates from celestial coordinates to the
intermediate frame, the Earth rotation angle and the polar motion
matrix. One use of the present function is when generating a
series of celestial-to-terrestrial matrices where only the Earth
Rotation Angle changes, avoiding the considerable overhead of
recomputing the precession-nutation more often than necessary to
achieve given accuracy objectives.
2) The relationship between the arguments is as follows:
[TRS] = RPOM * R_3(ERA) * rc2i * [CRS]
= rc2t * [CRS]
where [CRS] is a vector in the Geocentric Celestial Reference
System and [TRS] is a vector in the International Terrestrial
Reference System (see IERS Conventions 2003).
Called:
Cr copy r-matrix