1 SUBROUTINE sla_CLDJ (IY, IM, ID, DJM, J) 2*+ 3* - - - - - 4* C L D J 5* - - - - - 6* 7* Gregorian Calendar to Modified Julian Date 8* 9* Given: 10* IY,IM,ID int year, month, day in Gregorian calendar 11* 12* Returned: 13* DJM dp modified Julian Date (JD-2400000.5) for 0 hrs 14* J int status: 15* 0 = OK 16* 1 = bad year (MJD not computed) 17* 2 = bad month (MJD not computed) 18* 3 = bad day (MJD computed) 19* 20* The year must be -4699 (i.e. 4700BC) or later. 21* 22* The algorithm is adapted from Hatcher 1984 (QJRAS 25, 53-55). 23* 24* Last revision: 27 July 2004 25* 26* Copyright P.T.Wallace. All rights reserved. 27* 28* License: 29* This program is free software; you can redistribute it and/or modify 30* it under the terms of the GNU General Public License as published by 31* the Free Software Foundation; either version 2 of the License, or 32* (at your option) any later version. 33* 34* This program is distributed in the hope that it will be useful, 35* but WITHOUT ANY WARRANTY; without even the implied warranty of 36* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 37* GNU General Public License for more details. 38* 39* You should have received a copy of the GNU General Public License 40* along with this program (see SLA_CONDITIONS); if not, write to the 41* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 42* Boston, MA 02110-1301 USA 43* 44*- 45 46 IMPLICIT NONE 47 48 INTEGER IY,IM,ID 49 DOUBLE PRECISION DJM 50 INTEGER J 51 52* Month lengths in days 53 INTEGER MTAB(12) 54 DATA MTAB / 31,28,31,30,31,30,31,31,30,31,30,31 / 55 56 57 58* Preset status. 59 J = 0 60 61* Validate year. 62 IF ( IY .LT. -4699 ) THEN 63 J = 1 64 ELSE 65 66* Validate month. 67 IF ( IM.GE.1 .AND. IM.LE.12 ) THEN 68 69* Allow for leap year. 70 IF ( MOD(IY,4) .EQ. 0 ) THEN 71 MTAB(2) = 29 72 ELSE 73 MTAB(2) = 28 74 END IF 75 IF ( MOD(IY,100).EQ.0 .AND. MOD(IY,400).NE.0 ) 76 : MTAB(2) = 28 77 78* Validate day. 79 IF ( ID.LT.1 .OR. ID.GT.MTAB(IM) ) J=3 80 81* Modified Julian Date. 82 DJM = DBLE ( ( 1461 * ( IY - (12-IM)/10 + 4712 ) ) / 4 83 : + ( 306 * MOD ( IM+9, 12 ) + 5 ) / 10 84 : - ( 3 * ( ( IY - (12-IM)/10 + 4900 ) / 100 ) ) / 4 85 : + ID - 2399904 ) 86 87* Bad month. 88 ELSE 89 J=2 90 END IF 91 92 END IF 93 94 END 95 DOUBLE PRECISION FUNCTION sla_DAT (UTC) 96*+ 97* - - - - 98* D A T 99* - - - - 100* 101* Increment to be applied to Coordinated Universal Time UTC to give 102* International Atomic Time TAI (double precision) 103* 104* Given: 105* UTC d UTC date as a modified JD (JD-2400000.5) 106* 107* Result: TAI-UTC in seconds 108* 109* Notes: 110* 111* 1 The UTC is specified to be a date rather than a time to indicate 112* that care needs to be taken not to specify an instant which lies 113* within a leap second. Though in most cases UTC can include the 114* fractional part, correct behaviour on the day of a leap second 115* can only be guaranteed up to the end of the second 23:59:59. 116* 117* 2 For epochs from 1961 January 1 onwards, the expressions from the 118* file ftp://maia.usno.navy.mil/ser7/tai-utc.dat are used. 119* 120* 3 The 5ms time step at 1961 January 1 is taken from 2.58.1 (p87) of 121* the 1992 Explanatory Supplement. 122* 123* 4 UTC began at 1960 January 1.0 (JD 2436934.5) and it is improper 124* to call the routine with an earlier epoch. However, if this 125* is attempted, the TAI-UTC expression for the year 1960 is used. 126* 127* 128* :-----------------------------------------: 129* : : 130* : IMPORTANT : 131* : : 132* : This routine must be updated on each : 133* : occasion that a leap second is : 134* : announced : 135* : : 136* : Latest leap second: 2015 July 1 : 137* : : 138* :-----------------------------------------: 139* 140* Last revision: 5 July 2008 141* 142* Copyright P.T.Wallace. All rights reserved. 143*- 144 145 IMPLICIT NONE 146 147 DOUBLE PRECISION UTC 148 149 DOUBLE PRECISION DT 150 151 152 153 IF (.FALSE.) THEN 154 155* - - - - - - - - - - - - - - - - - - - - - - * 156* Add new code here on each occasion that a * 157* leap second is announced, and update the * 158* preamble comments appropriately. * 159* - - - - - - - - - - - - - - - - - - - - - - * 160 161* 2015 July 1 162 ELSE IF (UTC.GE.57204D0) THEN 163 DT=36D0 164 165* 2012 July 1 166 ELSE IF (UTC.GE.56109D0) THEN 167 DT=35D0 168 169* 2009 January 1 170 ELSE IF (UTC.GE.54832D0) THEN 171 DT=34D0 172 173* 2006 January 1 174 ELSE IF (UTC.GE.53736D0) THEN 175 DT=33D0 176 177* 1999 January 1 178 ELSE IF (UTC.GE.51179D0) THEN 179 DT=32D0 180 181* 1997 July 1 182 ELSE IF (UTC.GE.50630D0) THEN 183 DT=31D0 184 185* 1996 January 1 186 ELSE IF (UTC.GE.50083D0) THEN 187 DT=30D0 188 189* 1994 July 1 190 ELSE IF (UTC.GE.49534D0) THEN 191 DT=29D0 192 193* 1993 July 1 194 ELSE IF (UTC.GE.49169D0) THEN 195 DT=28D0 196 197* 1992 July 1 198 ELSE IF (UTC.GE.48804D0) THEN 199 DT=27D0 200 201* 1991 January 1 202 ELSE IF (UTC.GE.48257D0) THEN 203 DT=26D0 204 205* 1990 January 1 206 ELSE IF (UTC.GE.47892D0) THEN 207 DT=25D0 208 209* 1988 January 1 210 ELSE IF (UTC.GE.47161D0) THEN 211 DT=24D0 212 213* 1985 July 1 214 ELSE IF (UTC.GE.46247D0) THEN 215 DT=23D0 216 217* 1983 July 1 218 ELSE IF (UTC.GE.45516D0) THEN 219 DT=22D0 220 221* 1982 July 1 222 ELSE IF (UTC.GE.45151D0) THEN 223 DT=21D0 224 225* 1981 July 1 226 ELSE IF (UTC.GE.44786D0) THEN 227 DT=20D0 228 229* 1980 January 1 230 ELSE IF (UTC.GE.44239D0) THEN 231 DT=19D0 232 233* 1979 January 1 234 ELSE IF (UTC.GE.43874D0) THEN 235 DT=18D0 236 237* 1978 January 1 238 ELSE IF (UTC.GE.43509D0) THEN 239 DT=17D0 240 241* 1977 January 1 242 ELSE IF (UTC.GE.43144D0) THEN 243 DT=16D0 244 245* 1976 January 1 246 ELSE IF (UTC.GE.42778D0) THEN 247 DT=15D0 248 249* 1975 January 1 250 ELSE IF (UTC.GE.42413D0) THEN 251 DT=14D0 252 253* 1974 January 1 254 ELSE IF (UTC.GE.42048D0) THEN 255 DT=13D0 256 257* 1973 January 1 258 ELSE IF (UTC.GE.41683D0) THEN 259 DT=12D0 260 261* 1972 July 1 262 ELSE IF (UTC.GE.41499D0) THEN 263 DT=11D0 264 265* 1972 January 1 266 ELSE IF (UTC.GE.41317D0) THEN 267 DT=10D0 268 269* 1968 February 1 270 ELSE IF (UTC.GE.39887D0) THEN 271 DT=4.2131700D0+(UTC-39126D0)*0.002592D0 272 273* 1966 January 1 274 ELSE IF (UTC.GE.39126D0) THEN 275 DT=4.3131700D0+(UTC-39126D0)*0.002592D0 276 277* 1965 September 1 278 ELSE IF (UTC.GE.39004D0) THEN 279 DT=3.8401300D0+(UTC-38761D0)*0.001296D0 280 281* 1965 July 1 282 ELSE IF (UTC.GE.38942D0) THEN 283 DT=3.7401300D0+(UTC-38761D0)*0.001296D0 284 285* 1965 March 1 286 ELSE IF (UTC.GE.38820D0) THEN 287 DT=3.6401300D0+(UTC-38761D0)*0.001296D0 288 289* 1965 January 1 290 ELSE IF (UTC.GE.38761D0) THEN 291 DT=3.5401300D0+(UTC-38761D0)*0.001296D0 292 293* 1964 September 1 294 ELSE IF (UTC.GE.38639D0) THEN 295 DT=3.4401300D0+(UTC-38761D0)*0.001296D0 296 297* 1964 April 1 298 ELSE IF (UTC.GE.38486D0) THEN 299 DT=3.3401300D0+(UTC-38761D0)*0.001296D0 300 301* 1964 January 1 302 ELSE IF (UTC.GE.38395D0) THEN 303 DT=3.2401300D0+(UTC-38761D0)*0.001296D0 304 305* 1963 November 1 306 ELSE IF (UTC.GE.38334D0) THEN 307 DT=1.9458580D0+(UTC-37665D0)*0.0011232D0 308 309* 1962 January 1 310 ELSE IF (UTC.GE.37665D0) THEN 311 DT=1.8458580D0+(UTC-37665D0)*0.0011232D0 312 313* 1961 August 1 314 ELSE IF (UTC.GE.37512D0) THEN 315 DT=1.3728180D0+(UTC-37300D0)*0.001296D0 316 317* 1961 January 1 318 ELSE IF (UTC.GE.37300D0) THEN 319 DT=1.4228180D0+(UTC-37300D0)*0.001296D0 320 321* Before that 322 ELSE 323 DT=1.4178180D0+(UTC-37300D0)*0.001296D0 324 325 END IF 326 327 sla_DAT=DT 328 329 END 330 SUBROUTINE sla_DC62S (V, A, B, R, AD, BD, RD) 331*+ 332* - - - - - - 333* D C 6 2 S 334* - - - - - - 335* 336* Conversion of position & velocity in Cartesian coordinates 337* to spherical coordinates (double precision) 338* 339* Given: 340* V d(6) Cartesian position & velocity vector 341* 342* Returned: 343* A d longitude (radians) 344* B d latitude (radians) 345* R d radial coordinate 346* AD d longitude derivative (radians per unit time) 347* BD d latitude derivative (radians per unit time) 348* RD d radial derivative 349* 350* P.T.Wallace Starlink 28 April 1996 351* 352* Copyright (C) 1996 Rutherford Appleton Laboratory 353* 354* License: 355* This program is free software; you can redistribute it and/or modify 356* it under the terms of the GNU General Public License as published by 357* the Free Software Foundation; either version 2 of the License, or 358* (at your option) any later version. 359* 360* This program is distributed in the hope that it will be useful, 361* but WITHOUT ANY WARRANTY; without even the implied warranty of 362* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 363* GNU General Public License for more details. 364* 365* You should have received a copy of the GNU General Public License 366* along with this program (see SLA_CONDITIONS); if not, write to the 367* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 368* Boston, MA 02110-1301 USA 369* 370*- 371 372 IMPLICIT NONE 373 374 DOUBLE PRECISION V(6),A,B,R,AD,BD,RD 375 376 DOUBLE PRECISION X,Y,Z,XD,YD,ZD,RXY2,RXY,R2,XYP 377 378 379 380* Components of position/velocity vector 381 X=V(1) 382 Y=V(2) 383 Z=V(3) 384 XD=V(4) 385 YD=V(5) 386 ZD=V(6) 387 388* Component of R in XY plane squared 389 RXY2=X*X+Y*Y 390 391* Modulus squared 392 R2=RXY2+Z*Z 393 394* Protection against null vector 395 IF (R2.EQ.0D0) THEN 396 X=XD 397 Y=YD 398 Z=ZD 399 RXY2=X*X+Y*Y 400 R2=RXY2+Z*Z 401 END IF 402 403* Position and velocity in spherical coordinates 404 RXY=SQRT(RXY2) 405 XYP=X*XD+Y*YD 406 IF (RXY2.NE.0D0) THEN 407 A=ATAN2(Y,X) 408 B=ATAN2(Z,RXY) 409 AD=(X*YD-Y*XD)/RXY2 410 BD=(ZD*RXY2-Z*XYP)/(R2*RXY) 411 ELSE 412 A=0D0 413 IF (Z.NE.0D0) THEN 414 B=ATAN2(Z,RXY) 415 ELSE 416 B=0D0 417 END IF 418 AD=0D0 419 BD=0D0 420 END IF 421 R=SQRT(R2) 422 IF (R.NE.0D0) THEN 423 RD=(XYP+Z*ZD)/R 424 ELSE 425 RD=0D0 426 END IF 427 428 END 429 SUBROUTINE sla_DCC2S (V, A, B) 430*+ 431* - - - - - - 432* D C C 2 S 433* - - - - - - 434* 435* Cartesian to spherical coordinates (double precision) 436* 437* Given: 438* V d(3) x,y,z vector 439* 440* Returned: 441* A,B d spherical coordinates in radians 442* 443* The spherical coordinates are longitude (+ve anticlockwise looking 444* from the +ve latitude pole) and latitude. The Cartesian coordinates 445* are right handed, with the x axis at zero longitude and latitude, and 446* the z axis at the +ve latitude pole. 447* 448* If V is null, zero A and B are returned. At either pole, zero A is 449* returned. 450* 451* Last revision: 22 July 2004 452* 453* Copyright P.T.Wallace. All rights reserved. 454* 455* License: 456* This program is free software; you can redistribute it and/or modify 457* it under the terms of the GNU General Public License as published by 458* the Free Software Foundation; either version 2 of the License, or 459* (at your option) any later version. 460* 461* This program is distributed in the hope that it will be useful, 462* but WITHOUT ANY WARRANTY; without even the implied warranty of 463* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 464* GNU General Public License for more details. 465* 466* You should have received a copy of the GNU General Public License 467* along with this program (see SLA_CONDITIONS); if not, write to the 468* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 469* Boston, MA 02110-1301 USA 470* 471*- 472 473 IMPLICIT NONE 474 475 DOUBLE PRECISION V(3),A,B 476 477 DOUBLE PRECISION X,Y,Z,R 478 479 480 X = V(1) 481 Y = V(2) 482 Z = V(3) 483 R = SQRT(X*X+Y*Y) 484 485 IF (R.EQ.0D0) THEN 486 A = 0D0 487 ELSE 488 A = ATAN2(Y,X) 489 END IF 490 491 IF (Z.EQ.0D0) THEN 492 B = 0D0 493 ELSE 494 B = ATAN2(Z,R) 495 END IF 496 497 END 498 SUBROUTINE sla_DCS2C (A, B, V) 499*+ 500* - - - - - - 501* D C S 2 C 502* - - - - - - 503* 504* Spherical coordinates to direction cosines (double precision) 505* 506* Given: 507* A,B d spherical coordinates in radians 508* (RA,Dec), (long,lat) etc. 509* 510* Returned: 511* V d(3) x,y,z unit vector 512* 513* The spherical coordinates are longitude (+ve anticlockwise looking 514* from the +ve latitude pole) and latitude. The Cartesian coordinates 515* are right handed, with the x axis at zero longitude and latitude, and 516* the z axis at the +ve latitude pole. 517* 518* Last revision: 26 December 2004 519* 520* Copyright P.T.Wallace. All rights reserved. 521* 522* License: 523* This program is free software; you can redistribute it and/or modify 524* it under the terms of the GNU General Public License as published by 525* the Free Software Foundation; either version 2 of the License, or 526* (at your option) any later version. 527* 528* This program is distributed in the hope that it will be useful, 529* but WITHOUT ANY WARRANTY; without even the implied warranty of 530* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 531* GNU General Public License for more details. 532* 533* You should have received a copy of the GNU General Public License 534* along with this program (see SLA_CONDITIONS); if not, write to the 535* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 536* Boston, MA 02110-1301 USA 537* 538*- 539 540 IMPLICIT NONE 541 542 DOUBLE PRECISION A,B,V(3) 543 544 DOUBLE PRECISION COSB 545 546 547 COSB = COS(B) 548 549 V(1) = COS(A)*COSB 550 V(2) = SIN(A)*COSB 551 V(3) = SIN(B) 552 553 END 554 SUBROUTINE sla_DE2H (HA, DEC, PHI, AZ, EL) 555*+ 556* - - - - - 557* D E 2 H 558* - - - - - 559* 560* Equatorial to horizon coordinates: HA,Dec to Az,El 561* 562* (double precision) 563* 564* Given: 565* HA d hour angle 566* DEC d declination 567* PHI d observatory latitude 568* 569* Returned: 570* AZ d azimuth 571* EL d elevation 572* 573* Notes: 574* 575* 1) All the arguments are angles in radians. 576* 577* 2) Azimuth is returned in the range 0-2pi; north is zero, 578* and east is +pi/2. Elevation is returned in the range 579* +/-pi/2. 580* 581* 3) The latitude must be geodetic. In critical applications, 582* corrections for polar motion should be applied. 583* 584* 4) In some applications it will be important to specify the 585* correct type of hour angle and declination in order to 586* produce the required type of azimuth and elevation. In 587* particular, it may be important to distinguish between 588* elevation as affected by refraction, which would 589* require the "observed" HA,Dec, and the elevation 590* in vacuo, which would require the "topocentric" HA,Dec. 591* If the effects of diurnal aberration can be neglected, the 592* "apparent" HA,Dec may be used instead of the topocentric 593* HA,Dec. 594* 595* 5) No range checking of arguments is carried out. 596* 597* 6) In applications which involve many such calculations, rather 598* than calling the present routine it will be more efficient to 599* use inline code, having previously computed fixed terms such 600* as sine and cosine of latitude, and (for tracking a star) 601* sine and cosine of declination. 602* 603* P.T.Wallace Starlink 9 July 1994 604* 605* Copyright (C) 1995 Rutherford Appleton Laboratory 606* 607* License: 608* This program is free software; you can redistribute it and/or modify 609* it under the terms of the GNU General Public License as published by 610* the Free Software Foundation; either version 2 of the License, or 611* (at your option) any later version. 612* 613* This program is distributed in the hope that it will be useful, 614* but WITHOUT ANY WARRANTY; without even the implied warranty of 615* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 616* GNU General Public License for more details. 617* 618* You should have received a copy of the GNU General Public License 619* along with this program (see SLA_CONDITIONS); if not, write to the 620* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 621* Boston, MA 02110-1301 USA 622* 623*- 624 625 IMPLICIT NONE 626 627 DOUBLE PRECISION HA,DEC,PHI,AZ,EL 628 629 DOUBLE PRECISION D2PI 630 PARAMETER (D2PI=6.283185307179586476925286766559D0) 631 632 DOUBLE PRECISION SH,CH,SD,CD,SP,CP,X,Y,Z,R,A 633 634 635* Useful trig functions 636 SH=SIN(HA) 637 CH=COS(HA) 638 SD=SIN(DEC) 639 CD=COS(DEC) 640 SP=SIN(PHI) 641 CP=COS(PHI) 642 643* Az,El as x,y,z 644 X=-CH*CD*SP+SD*CP 645 Y=-SH*CD 646 Z=CH*CD*CP+SD*SP 647 648* To spherical 649 R=SQRT(X*X+Y*Y) 650 IF (R.EQ.0D0) THEN 651 A=0D0 652 ELSE 653 A=ATAN2(Y,X) 654 END IF 655 IF (A.LT.0D0) A=A+D2PI 656 AZ=A 657 EL=ATAN2(Z,R) 658 659 END 660 SUBROUTINE sla_DEULER (ORDER, PHI, THETA, PSI, RMAT) 661*+ 662* - - - - - - - 663* D E U L E R 664* - - - - - - - 665* 666* Form a rotation matrix from the Euler angles - three successive 667* rotations about specified Cartesian axes (double precision) 668* 669* Given: 670* ORDER c*(*) specifies about which axes the rotations occur 671* PHI d 1st rotation (radians) 672* THETA d 2nd rotation ( " ) 673* PSI d 3rd rotation ( " ) 674* 675* Returned: 676* RMAT d(3,3) rotation matrix 677* 678* A rotation is positive when the reference frame rotates 679* anticlockwise as seen looking towards the origin from the 680* positive region of the specified axis. 681* 682* The characters of ORDER define which axes the three successive 683* rotations are about. A typical value is 'ZXZ', indicating that 684* RMAT is to become the direction cosine matrix corresponding to 685* rotations of the reference frame through PHI radians about the 686* old Z-axis, followed by THETA radians about the resulting X-axis, 687* then PSI radians about the resulting Z-axis. 688* 689* The axis names can be any of the following, in any order or 690* combination: X, Y, Z, uppercase or lowercase, 1, 2, 3. Normal 691* axis labelling/numbering conventions apply; the xyz (=123) 692* triad is right-handed. Thus, the 'ZXZ' example given above 693* could be written 'zxz' or '313' (or even 'ZxZ' or '3xZ'). ORDER 694* is terminated by length or by the first unrecognized character. 695* 696* Fewer than three rotations are acceptable, in which case the later 697* angle arguments are ignored. If all rotations are zero, the 698* identity matrix is produced. 699* 700* P.T.Wallace Starlink 23 May 1997 701* 702* Copyright (C) 1997 Rutherford Appleton Laboratory 703* 704* License: 705* This program is free software; you can redistribute it and/or modify 706* it under the terms of the GNU General Public License as published by 707* the Free Software Foundation; either version 2 of the License, or 708* (at your option) any later version. 709* 710* This program is distributed in the hope that it will be useful, 711* but WITHOUT ANY WARRANTY; without even the implied warranty of 712* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 713* GNU General Public License for more details. 714* 715* You should have received a copy of the GNU General Public License 716* along with this program (see SLA_CONDITIONS); if not, write to the 717* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 718* Boston, MA 02110-1301 USA 719* 720*- 721 722 IMPLICIT NONE 723 724 CHARACTER*(*) ORDER 725 DOUBLE PRECISION PHI,THETA,PSI,RMAT(3,3) 726 727 INTEGER J,I,L,N,K 728 DOUBLE PRECISION RESULT(3,3),ROTN(3,3),ANGLE,S,C,W,WM(3,3) 729 CHARACTER AXIS 730 731 732 733* Initialize result matrix 734 DO J=1,3 735 DO I=1,3 736 IF (I.NE.J) THEN 737 RESULT(I,J) = 0D0 738 ELSE 739 RESULT(I,J) = 1D0 740 END IF 741 END DO 742 END DO 743 744* Establish length of axis string 745 L = LEN(ORDER) 746 747* Look at each character of axis string until finished 748 DO N=1,3 749 IF (N.LE.L) THEN 750 751* Initialize rotation matrix for the current rotation 752 DO J=1,3 753 DO I=1,3 754 IF (I.NE.J) THEN 755 ROTN(I,J) = 0D0 756 ELSE 757 ROTN(I,J) = 1D0 758 END IF 759 END DO 760 END DO 761 762* Pick up the appropriate Euler angle and take sine & cosine 763 IF (N.EQ.1) THEN 764 ANGLE = PHI 765 ELSE IF (N.EQ.2) THEN 766 ANGLE = THETA 767 ELSE 768 ANGLE = PSI 769 END IF 770 S = SIN(ANGLE) 771 C = COS(ANGLE) 772 773* Identify the axis 774 AXIS = ORDER(N:N) 775 IF (AXIS.EQ.'X'.OR. 776 : AXIS.EQ.'x'.OR. 777 : AXIS.EQ.'1') THEN 778 779* Matrix for x-rotation 780 ROTN(2,2) = C 781 ROTN(2,3) = S 782 ROTN(3,2) = -S 783 ROTN(3,3) = C 784 785 ELSE IF (AXIS.EQ.'Y'.OR. 786 : AXIS.EQ.'y'.OR. 787 : AXIS.EQ.'2') THEN 788 789* Matrix for y-rotation 790 ROTN(1,1) = C 791 ROTN(1,3) = -S 792 ROTN(3,1) = S 793 ROTN(3,3) = C 794 795 ELSE IF (AXIS.EQ.'Z'.OR. 796 : AXIS.EQ.'z'.OR. 797 : AXIS.EQ.'3') THEN 798 799* Matrix for z-rotation 800 ROTN(1,1) = C 801 ROTN(1,2) = S 802 ROTN(2,1) = -S 803 ROTN(2,2) = C 804 805 ELSE 806 807* Unrecognized character - fake end of string 808 L = 0 809 810 END IF 811 812* Apply the current rotation (matrix ROTN x matrix RESULT) 813 DO I=1,3 814 DO J=1,3 815 W = 0D0 816 DO K=1,3 817 W = W+ROTN(I,K)*RESULT(K,J) 818 END DO 819 WM(I,J) = W 820 END DO 821 END DO 822 DO J=1,3 823 DO I=1,3 824 RESULT(I,J) = WM(I,J) 825 END DO 826 END DO 827 828 END IF 829 830 END DO 831 832* Copy the result 833 DO J=1,3 834 DO I=1,3 835 RMAT(I,J) = RESULT(I,J) 836 END DO 837 END DO 838 839 END 840 SUBROUTINE sla_DMOON (DATE, PV) 841*+ 842* - - - - - - 843* D M O O N 844* - - - - - - 845* 846* Approximate geocentric position and velocity of the Moon 847* (double precision) 848* 849* Given: 850* DATE D TDB (loosely ET) as a Modified Julian Date 851* (JD-2400000.5) 852* 853* Returned: 854* PV D(6) Moon x,y,z,xdot,ydot,zdot, mean equator and 855* equinox of date (AU, AU/s) 856* 857* Notes: 858* 859* 1 This routine is a full implementation of the algorithm 860* published by Meeus (see reference). 861* 862* 2 Meeus quotes accuracies of 10 arcsec in longitude, 3 arcsec in 863* latitude and 0.2 arcsec in HP (equivalent to about 20 km in 864* distance). Comparison with JPL DE200 over the interval 865* 1960-2025 gives RMS errors of 3.7 arcsec and 83 mas/hour in 866* longitude, 2.3 arcsec and 48 mas/hour in latitude, 11 km 867* and 81 mm/s in distance. The maximum errors over the same 868* interval are 18 arcsec and 0.50 arcsec/hour in longitude, 869* 11 arcsec and 0.24 arcsec/hour in latitude, 40 km and 0.29 m/s 870* in distance. 871* 872* 3 The original algorithm is expressed in terms of the obsolete 873* timescale Ephemeris Time. Either TDB or TT can be used, but 874* not UT without incurring significant errors (30 arcsec at 875* the present time) due to the Moon's 0.5 arcsec/sec movement. 876* 877* 4 The algorithm is based on pre IAU 1976 standards. However, 878* the result has been moved onto the new (FK5) equinox, an 879* adjustment which is in any case much smaller than the 880* intrinsic accuracy of the procedure. 881* 882* 5 Velocity is obtained by a complete analytical differentiation 883* of the Meeus model. 884* 885* Reference: 886* Meeus, l'Astronomie, June 1984, p348. 887* 888* P.T.Wallace Starlink 22 January 1998 889* 890* Copyright (C) 1998 Rutherford Appleton Laboratory 891* 892* License: 893* This program is free software; you can redistribute it and/or modify 894* it under the terms of the GNU General Public License as published by 895* the Free Software Foundation; either version 2 of the License, or 896* (at your option) any later version. 897* 898* This program is distributed in the hope that it will be useful, 899* but WITHOUT ANY WARRANTY; without even the implied warranty of 900* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 901* GNU General Public License for more details. 902* 903* You should have received a copy of the GNU General Public License 904* along with this program (see SLA_CONDITIONS); if not, write to the 905* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 906* Boston, MA 02110-1301 USA 907* 908*- 909 910 IMPLICIT NONE 911 912 DOUBLE PRECISION DATE,PV(6) 913 914* Degrees, arcseconds and seconds of time to radians 915 DOUBLE PRECISION D2R,DAS2R,DS2R 916 PARAMETER (D2R=0.0174532925199432957692369D0, 917 : DAS2R=4.848136811095359935899141D-6, 918 : DS2R=7.272205216643039903848712D-5) 919 920* Seconds per Julian century (86400*36525) 921 DOUBLE PRECISION CJ 922 PARAMETER (CJ=3155760000D0) 923 924* Julian epoch of B1950 925 DOUBLE PRECISION B1950 926 PARAMETER (B1950=1949.9997904423D0) 927 928* Earth equatorial radius in AU ( = 6378.137 / 149597870 ) 929 DOUBLE PRECISION ERADAU 930 PARAMETER (ERADAU=4.2635212653763D-5) 931 932 DOUBLE PRECISION T,THETA,SINOM,COSOM,DOMCOM,WA,DWA,WB,DWB,WOM, 933 : DWOM,SINWOM,COSWOM,V,DV,COEFF,EMN,EMPN,DN,FN,EN, 934 : DEN,DTHETA,FTHETA,EL,DEL,B,DB,BF,DBF,P,DP,SP,R, 935 : DR,X,Y,Z,XD,YD,ZD,SEL,CEL,SB,CB,RCB,RBD,W,EPJ, 936 : EQCOR,EPS,SINEPS,COSEPS,ES,EC 937 INTEGER N,I 938 939* 940* Coefficients for fundamental arguments 941* 942* at J1900: T**0, T**1, T**2, T**3 943* at epoch: T**0, T**1 944* 945* Units are degrees for position and Julian centuries for time 946* 947 948* Moon's mean longitude 949 DOUBLE PRECISION ELP0,ELP1,ELP2,ELP3,ELP,DELP 950 PARAMETER (ELP0=270.434164D0, 951 : ELP1=481267.8831D0, 952 : ELP2=-0.001133D0, 953 : ELP3=0.0000019D0) 954 955* Sun's mean anomaly 956 DOUBLE PRECISION EM0,EM1,EM2,EM3,EM,DEM 957 PARAMETER (EM0=358.475833D0, 958 : EM1=35999.0498D0, 959 : EM2=-0.000150D0, 960 : EM3=-0.0000033D0) 961 962* Moon's mean anomaly 963 DOUBLE PRECISION EMP0,EMP1,EMP2,EMP3,EMP,DEMP 964 PARAMETER (EMP0=296.104608D0, 965 : EMP1=477198.8491D0, 966 : EMP2=0.009192D0, 967 : EMP3=0.0000144D0) 968 969* Moon's mean elongation 970 DOUBLE PRECISION D0,D1,D2,D3,D,DD 971 PARAMETER (D0=350.737486D0, 972 : D1=445267.1142D0, 973 : D2=-0.001436D0, 974 : D3=0.0000019D0) 975 976* Mean distance of the Moon from its ascending node 977 DOUBLE PRECISION F0,F1,F2,F3,F,DF 978 PARAMETER (F0=11.250889D0, 979 : F1=483202.0251D0, 980 : F2=-0.003211D0, 981 : F3=-0.0000003D0) 982 983* Longitude of the Moon's ascending node 984 DOUBLE PRECISION OM0,OM1,OM2,OM3,OM,DOM 985 PARAMETER (OM0=259.183275D0, 986 : OM1=-1934.1420D0, 987 : OM2=0.002078D0, 988 : OM3=0.0000022D0) 989 990* Coefficients for (dimensionless) E factor 991 DOUBLE PRECISION E1,E2,E,DE,ESQ,DESQ 992 PARAMETER (E1=-0.002495D0,E2=-0.00000752D0) 993 994* Coefficients for periodic variations etc 995 DOUBLE PRECISION PAC,PA0,PA1 996 PARAMETER (PAC=0.000233D0,PA0=51.2D0,PA1=20.2D0) 997 DOUBLE PRECISION PBC 998 PARAMETER (PBC=-0.001778D0) 999 DOUBLE PRECISION PCC 1000 PARAMETER (PCC=0.000817D0) 1001 DOUBLE PRECISION PDC 1002 PARAMETER (PDC=0.002011D0) 1003 DOUBLE PRECISION PEC,PE0,PE1,PE2 1004 PARAMETER (PEC=0.003964D0, 1005 : PE0=346.560D0,PE1=132.870D0,PE2=-0.0091731D0) 1006 DOUBLE PRECISION PFC 1007 PARAMETER (PFC=0.001964D0) 1008 DOUBLE PRECISION PGC 1009 PARAMETER (PGC=0.002541D0) 1010 DOUBLE PRECISION PHC 1011 PARAMETER (PHC=0.001964D0) 1012 DOUBLE PRECISION PIC 1013 PARAMETER (PIC=-0.024691D0) 1014 DOUBLE PRECISION PJC,PJ0,PJ1 1015 PARAMETER (PJC=-0.004328D0,PJ0=275.05D0,PJ1=-2.30D0) 1016 DOUBLE PRECISION CW1 1017 PARAMETER (CW1=0.0004664D0) 1018 DOUBLE PRECISION CW2 1019 PARAMETER (CW2=0.0000754D0) 1020 1021* 1022* Coefficients for Moon position 1023* 1024* Tx(N) = coefficient of L, B or P term (deg) 1025* ITx(N,1-5) = coefficients of M, M', D, F, E**n in argument 1026* 1027 INTEGER NL,NB,NP 1028 PARAMETER (NL=50,NB=45,NP=31) 1029 DOUBLE PRECISION TL(NL),TB(NB),TP(NP) 1030 INTEGER ITL(5,NL),ITB(5,NB),ITP(5,NP) 1031* 1032* Longitude 1033* M M' D F n 1034 DATA TL( 1)/ +6.288750D0 /, 1035 : (ITL(I, 1),I=1,5)/ +0, +1, +0, +0, 0 / 1036 DATA TL( 2)/ +1.274018D0 /, 1037 : (ITL(I, 2),I=1,5)/ +0, -1, +2, +0, 0 / 1038 DATA TL( 3)/ +0.658309D0 /, 1039 : (ITL(I, 3),I=1,5)/ +0, +0, +2, +0, 0 / 1040 DATA TL( 4)/ +0.213616D0 /, 1041 : (ITL(I, 4),I=1,5)/ +0, +2, +0, +0, 0 / 1042 DATA TL( 5)/ -0.185596D0 /, 1043 : (ITL(I, 5),I=1,5)/ +1, +0, +0, +0, 1 / 1044 DATA TL( 6)/ -0.114336D0 /, 1045 : (ITL(I, 6),I=1,5)/ +0, +0, +0, +2, 0 / 1046 DATA TL( 7)/ +0.058793D0 /, 1047 : (ITL(I, 7),I=1,5)/ +0, -2, +2, +0, 0 / 1048 DATA TL( 8)/ +0.057212D0 /, 1049 : (ITL(I, 8),I=1,5)/ -1, -1, +2, +0, 1 / 1050 DATA TL( 9)/ +0.053320D0 /, 1051 : (ITL(I, 9),I=1,5)/ +0, +1, +2, +0, 0 / 1052 DATA TL(10)/ +0.045874D0 /, 1053 : (ITL(I,10),I=1,5)/ -1, +0, +2, +0, 1 / 1054 DATA TL(11)/ +0.041024D0 /, 1055 : (ITL(I,11),I=1,5)/ -1, +1, +0, +0, 1 / 1056 DATA TL(12)/ -0.034718D0 /, 1057 : (ITL(I,12),I=1,5)/ +0, +0, +1, +0, 0 / 1058 DATA TL(13)/ -0.030465D0 /, 1059 : (ITL(I,13),I=1,5)/ +1, +1, +0, +0, 1 / 1060 DATA TL(14)/ +0.015326D0 /, 1061 : (ITL(I,14),I=1,5)/ +0, +0, +2, -2, 0 / 1062 DATA TL(15)/ -0.012528D0 /, 1063 : (ITL(I,15),I=1,5)/ +0, +1, +0, +2, 0 / 1064 DATA TL(16)/ -0.010980D0 /, 1065 : (ITL(I,16),I=1,5)/ +0, -1, +0, +2, 0 / 1066 DATA TL(17)/ +0.010674D0 /, 1067 : (ITL(I,17),I=1,5)/ +0, -1, +4, +0, 0 / 1068 DATA TL(18)/ +0.010034D0 /, 1069 : (ITL(I,18),I=1,5)/ +0, +3, +0, +0, 0 / 1070 DATA TL(19)/ +0.008548D0 /, 1071 : (ITL(I,19),I=1,5)/ +0, -2, +4, +0, 0 / 1072 DATA TL(20)/ -0.007910D0 /, 1073 : (ITL(I,20),I=1,5)/ +1, -1, +2, +0, 1 / 1074 DATA TL(21)/ -0.006783D0 /, 1075 : (ITL(I,21),I=1,5)/ +1, +0, +2, +0, 1 / 1076 DATA TL(22)/ +0.005162D0 /, 1077 : (ITL(I,22),I=1,5)/ +0, +1, -1, +0, 0 / 1078 DATA TL(23)/ +0.005000D0 /, 1079 : (ITL(I,23),I=1,5)/ +1, +0, +1, +0, 1 / 1080 DATA TL(24)/ +0.004049D0 /, 1081 : (ITL(I,24),I=1,5)/ -1, +1, +2, +0, 1 / 1082 DATA TL(25)/ +0.003996D0 /, 1083 : (ITL(I,25),I=1,5)/ +0, +2, +2, +0, 0 / 1084 DATA TL(26)/ +0.003862D0 /, 1085 : (ITL(I,26),I=1,5)/ +0, +0, +4, +0, 0 / 1086 DATA TL(27)/ +0.003665D0 /, 1087 : (ITL(I,27),I=1,5)/ +0, -3, +2, +0, 0 / 1088 DATA TL(28)/ +0.002695D0 /, 1089 : (ITL(I,28),I=1,5)/ -1, +2, +0, +0, 1 / 1090 DATA TL(29)/ +0.002602D0 /, 1091 : (ITL(I,29),I=1,5)/ +0, +1, -2, -2, 0 / 1092 DATA TL(30)/ +0.002396D0 /, 1093 : (ITL(I,30),I=1,5)/ -1, -2, +2, +0, 1 / 1094 DATA TL(31)/ -0.002349D0 /, 1095 : (ITL(I,31),I=1,5)/ +0, +1, +1, +0, 0 / 1096 DATA TL(32)/ +0.002249D0 /, 1097 : (ITL(I,32),I=1,5)/ -2, +0, +2, +0, 2 / 1098 DATA TL(33)/ -0.002125D0 /, 1099 : (ITL(I,33),I=1,5)/ +1, +2, +0, +0, 1 / 1100 DATA TL(34)/ -0.002079D0 /, 1101 : (ITL(I,34),I=1,5)/ +2, +0, +0, +0, 2 / 1102 DATA TL(35)/ +0.002059D0 /, 1103 : (ITL(I,35),I=1,5)/ -2, -1, +2, +0, 2 / 1104 DATA TL(36)/ -0.001773D0 /, 1105 : (ITL(I,36),I=1,5)/ +0, +1, +2, -2, 0 / 1106 DATA TL(37)/ -0.001595D0 /, 1107 : (ITL(I,37),I=1,5)/ +0, +0, +2, +2, 0 / 1108 DATA TL(38)/ +0.001220D0 /, 1109 : (ITL(I,38),I=1,5)/ -1, -1, +4, +0, 1 / 1110 DATA TL(39)/ -0.001110D0 /, 1111 : (ITL(I,39),I=1,5)/ +0, +2, +0, +2, 0 / 1112 DATA TL(40)/ +0.000892D0 /, 1113 : (ITL(I,40),I=1,5)/ +0, +1, -3, +0, 0 / 1114 DATA TL(41)/ -0.000811D0 /, 1115 : (ITL(I,41),I=1,5)/ +1, +1, +2, +0, 1 / 1116 DATA TL(42)/ +0.000761D0 /, 1117 : (ITL(I,42),I=1,5)/ -1, -2, +4, +0, 1 / 1118 DATA TL(43)/ +0.000717D0 /, 1119 : (ITL(I,43),I=1,5)/ -2, +1, +0, +0, 2 / 1120 DATA TL(44)/ +0.000704D0 /, 1121 : (ITL(I,44),I=1,5)/ -2, +1, -2, +0, 2 / 1122 DATA TL(45)/ +0.000693D0 /, 1123 : (ITL(I,45),I=1,5)/ +1, -2, +2, +0, 1 / 1124 DATA TL(46)/ +0.000598D0 /, 1125 : (ITL(I,46),I=1,5)/ -1, +0, +2, -2, 1 / 1126 DATA TL(47)/ +0.000550D0 /, 1127 : (ITL(I,47),I=1,5)/ +0, +1, +4, +0, 0 / 1128 DATA TL(48)/ +0.000538D0 /, 1129 : (ITL(I,48),I=1,5)/ +0, +4, +0, +0, 0 / 1130 DATA TL(49)/ +0.000521D0 /, 1131 : (ITL(I,49),I=1,5)/ -1, +0, +4, +0, 1 / 1132 DATA TL(50)/ +0.000486D0 /, 1133 : (ITL(I,50),I=1,5)/ +0, +2, -1, +0, 0 / 1134* 1135* Latitude 1136* M M' D F n 1137 DATA TB( 1)/ +5.128189D0 /, 1138 : (ITB(I, 1),I=1,5)/ +0, +0, +0, +1, 0 / 1139 DATA TB( 2)/ +0.280606D0 /, 1140 : (ITB(I, 2),I=1,5)/ +0, +1, +0, +1, 0 / 1141 DATA TB( 3)/ +0.277693D0 /, 1142 : (ITB(I, 3),I=1,5)/ +0, +1, +0, -1, 0 / 1143 DATA TB( 4)/ +0.173238D0 /, 1144 : (ITB(I, 4),I=1,5)/ +0, +0, +2, -1, 0 / 1145 DATA TB( 5)/ +0.055413D0 /, 1146 : (ITB(I, 5),I=1,5)/ +0, -1, +2, +1, 0 / 1147 DATA TB( 6)/ +0.046272D0 /, 1148 : (ITB(I, 6),I=1,5)/ +0, -1, +2, -1, 0 / 1149 DATA TB( 7)/ +0.032573D0 /, 1150 : (ITB(I, 7),I=1,5)/ +0, +0, +2, +1, 0 / 1151 DATA TB( 8)/ +0.017198D0 /, 1152 : (ITB(I, 8),I=1,5)/ +0, +2, +0, +1, 0 / 1153 DATA TB( 9)/ +0.009267D0 /, 1154 : (ITB(I, 9),I=1,5)/ +0, +1, +2, -1, 0 / 1155 DATA TB(10)/ +0.008823D0 /, 1156 : (ITB(I,10),I=1,5)/ +0, +2, +0, -1, 0 / 1157 DATA TB(11)/ +0.008247D0 /, 1158 : (ITB(I,11),I=1,5)/ -1, +0, +2, -1, 1 / 1159 DATA TB(12)/ +0.004323D0 /, 1160 : (ITB(I,12),I=1,5)/ +0, -2, +2, -1, 0 / 1161 DATA TB(13)/ +0.004200D0 /, 1162 : (ITB(I,13),I=1,5)/ +0, +1, +2, +1, 0 / 1163 DATA TB(14)/ +0.003372D0 /, 1164 : (ITB(I,14),I=1,5)/ -1, +0, -2, +1, 1 / 1165 DATA TB(15)/ +0.002472D0 /, 1166 : (ITB(I,15),I=1,5)/ -1, -1, +2, +1, 1 / 1167 DATA TB(16)/ +0.002222D0 /, 1168 : (ITB(I,16),I=1,5)/ -1, +0, +2, +1, 1 / 1169 DATA TB(17)/ +0.002072D0 /, 1170 : (ITB(I,17),I=1,5)/ -1, -1, +2, -1, 1 / 1171 DATA TB(18)/ +0.001877D0 /, 1172 : (ITB(I,18),I=1,5)/ -1, +1, +0, +1, 1 / 1173 DATA TB(19)/ +0.001828D0 /, 1174 : (ITB(I,19),I=1,5)/ +0, -1, +4, -1, 0 / 1175 DATA TB(20)/ -0.001803D0 /, 1176 : (ITB(I,20),I=1,5)/ +1, +0, +0, +1, 1 / 1177 DATA TB(21)/ -0.001750D0 /, 1178 : (ITB(I,21),I=1,5)/ +0, +0, +0, +3, 0 / 1179 DATA TB(22)/ +0.001570D0 /, 1180 : (ITB(I,22),I=1,5)/ -1, +1, +0, -1, 1 / 1181 DATA TB(23)/ -0.001487D0 /, 1182 : (ITB(I,23),I=1,5)/ +0, +0, +1, +1, 0 / 1183 DATA TB(24)/ -0.001481D0 /, 1184 : (ITB(I,24),I=1,5)/ +1, +1, +0, +1, 1 / 1185 DATA TB(25)/ +0.001417D0 /, 1186 : (ITB(I,25),I=1,5)/ -1, -1, +0, +1, 1 / 1187 DATA TB(26)/ +0.001350D0 /, 1188 : (ITB(I,26),I=1,5)/ -1, +0, +0, +1, 1 / 1189 DATA TB(27)/ +0.001330D0 /, 1190 : (ITB(I,27),I=1,5)/ +0, +0, -1, +1, 0 / 1191 DATA TB(28)/ +0.001106D0 /, 1192 : (ITB(I,28),I=1,5)/ +0, +3, +0, +1, 0 / 1193 DATA TB(29)/ +0.001020D0 /, 1194 : (ITB(I,29),I=1,5)/ +0, +0, +4, -1, 0 / 1195 DATA TB(30)/ +0.000833D0 /, 1196 : (ITB(I,30),I=1,5)/ +0, -1, +4, +1, 0 / 1197 DATA TB(31)/ +0.000781D0 /, 1198 : (ITB(I,31),I=1,5)/ +0, +1, +0, -3, 0 / 1199 DATA TB(32)/ +0.000670D0 /, 1200 : (ITB(I,32),I=1,5)/ +0, -2, +4, +1, 0 / 1201 DATA TB(33)/ +0.000606D0 /, 1202 : (ITB(I,33),I=1,5)/ +0, +0, +2, -3, 0 / 1203 DATA TB(34)/ +0.000597D0 /, 1204 : (ITB(I,34),I=1,5)/ +0, +2, +2, -1, 0 / 1205 DATA TB(35)/ +0.000492D0 /, 1206 : (ITB(I,35),I=1,5)/ -1, +1, +2, -1, 1 / 1207 DATA TB(36)/ +0.000450D0 /, 1208 : (ITB(I,36),I=1,5)/ +0, +2, -2, -1, 0 / 1209 DATA TB(37)/ +0.000439D0 /, 1210 : (ITB(I,37),I=1,5)/ +0, +3, +0, -1, 0 / 1211 DATA TB(38)/ +0.000423D0 /, 1212 : (ITB(I,38),I=1,5)/ +0, +2, +2, +1, 0 / 1213 DATA TB(39)/ +0.000422D0 /, 1214 : (ITB(I,39),I=1,5)/ +0, -3, +2, -1, 0 / 1215 DATA TB(40)/ -0.000367D0 /, 1216 : (ITB(I,40),I=1,5)/ +1, -1, +2, +1, 1 / 1217 DATA TB(41)/ -0.000353D0 /, 1218 : (ITB(I,41),I=1,5)/ +1, +0, +2, +1, 1 / 1219 DATA TB(42)/ +0.000331D0 /, 1220 : (ITB(I,42),I=1,5)/ +0, +0, +4, +1, 0 / 1221 DATA TB(43)/ +0.000317D0 /, 1222 : (ITB(I,43),I=1,5)/ -1, +1, +2, +1, 1 / 1223 DATA TB(44)/ +0.000306D0 /, 1224 : (ITB(I,44),I=1,5)/ -2, +0, +2, -1, 2 / 1225 DATA TB(45)/ -0.000283D0 /, 1226 : (ITB(I,45),I=1,5)/ +0, +1, +0, +3, 0 / 1227* 1228* Parallax 1229* M M' D F n 1230 DATA TP( 1)/ +0.950724D0 /, 1231 : (ITP(I, 1),I=1,5)/ +0, +0, +0, +0, 0 / 1232 DATA TP( 2)/ +0.051818D0 /, 1233 : (ITP(I, 2),I=1,5)/ +0, +1, +0, +0, 0 / 1234 DATA TP( 3)/ +0.009531D0 /, 1235 : (ITP(I, 3),I=1,5)/ +0, -1, +2, +0, 0 / 1236 DATA TP( 4)/ +0.007843D0 /, 1237 : (ITP(I, 4),I=1,5)/ +0, +0, +2, +0, 0 / 1238 DATA TP( 5)/ +0.002824D0 /, 1239 : (ITP(I, 5),I=1,5)/ +0, +2, +0, +0, 0 / 1240 DATA TP( 6)/ +0.000857D0 /, 1241 : (ITP(I, 6),I=1,5)/ +0, +1, +2, +0, 0 / 1242 DATA TP( 7)/ +0.000533D0 /, 1243 : (ITP(I, 7),I=1,5)/ -1, +0, +2, +0, 1 / 1244 DATA TP( 8)/ +0.000401D0 /, 1245 : (ITP(I, 8),I=1,5)/ -1, -1, +2, +0, 1 / 1246 DATA TP( 9)/ +0.000320D0 /, 1247 : (ITP(I, 9),I=1,5)/ -1, +1, +0, +0, 1 / 1248 DATA TP(10)/ -0.000271D0 /, 1249 : (ITP(I,10),I=1,5)/ +0, +0, +1, +0, 0 / 1250 DATA TP(11)/ -0.000264D0 /, 1251 : (ITP(I,11),I=1,5)/ +1, +1, +0, +0, 1 / 1252 DATA TP(12)/ -0.000198D0 /, 1253 : (ITP(I,12),I=1,5)/ +0, -1, +0, +2, 0 / 1254 DATA TP(13)/ +0.000173D0 /, 1255 : (ITP(I,13),I=1,5)/ +0, +3, +0, +0, 0 / 1256 DATA TP(14)/ +0.000167D0 /, 1257 : (ITP(I,14),I=1,5)/ +0, -1, +4, +0, 0 / 1258 DATA TP(15)/ -0.000111D0 /, 1259 : (ITP(I,15),I=1,5)/ +1, +0, +0, +0, 1 / 1260 DATA TP(16)/ +0.000103D0 /, 1261 : (ITP(I,16),I=1,5)/ +0, -2, +4, +0, 0 / 1262 DATA TP(17)/ -0.000084D0 /, 1263 : (ITP(I,17),I=1,5)/ +0, +2, -2, +0, 0 / 1264 DATA TP(18)/ -0.000083D0 /, 1265 : (ITP(I,18),I=1,5)/ +1, +0, +2, +0, 1 / 1266 DATA TP(19)/ +0.000079D0 /, 1267 : (ITP(I,19),I=1,5)/ +0, +2, +2, +0, 0 / 1268 DATA TP(20)/ +0.000072D0 /, 1269 : (ITP(I,20),I=1,5)/ +0, +0, +4, +0, 0 / 1270 DATA TP(21)/ +0.000064D0 /, 1271 : (ITP(I,21),I=1,5)/ -1, +1, +2, +0, 1 / 1272 DATA TP(22)/ -0.000063D0 /, 1273 : (ITP(I,22),I=1,5)/ +1, -1, +2, +0, 1 / 1274 DATA TP(23)/ +0.000041D0 /, 1275 : (ITP(I,23),I=1,5)/ +1, +0, +1, +0, 1 / 1276 DATA TP(24)/ +0.000035D0 /, 1277 : (ITP(I,24),I=1,5)/ -1, +2, +0, +0, 1 / 1278 DATA TP(25)/ -0.000033D0 /, 1279 : (ITP(I,25),I=1,5)/ +0, +3, -2, +0, 0 / 1280 DATA TP(26)/ -0.000030D0 /, 1281 : (ITP(I,26),I=1,5)/ +0, +1, +1, +0, 0 / 1282 DATA TP(27)/ -0.000029D0 /, 1283 : (ITP(I,27),I=1,5)/ +0, +0, -2, +2, 0 / 1284 DATA TP(28)/ -0.000029D0 /, 1285 : (ITP(I,28),I=1,5)/ +1, +2, +0, +0, 1 / 1286 DATA TP(29)/ +0.000026D0 /, 1287 : (ITP(I,29),I=1,5)/ -2, +0, +2, +0, 2 / 1288 DATA TP(30)/ -0.000023D0 /, 1289 : (ITP(I,30),I=1,5)/ +0, +1, -2, +2, 0 / 1290 DATA TP(31)/ +0.000019D0 /, 1291 : (ITP(I,31),I=1,5)/ -1, -1, +4, +0, 1 / 1292 1293 1294 1295* Centuries since J1900 1296 T=(DATE-15019.5D0)/36525D0 1297 1298* 1299* Fundamental arguments (radians) and derivatives (radians per 1300* Julian century) for the current epoch 1301* 1302 1303* Moon's mean longitude 1304 ELP=D2R*MOD(ELP0+(ELP1+(ELP2+ELP3*T)*T)*T,360D0) 1305 DELP=D2R*(ELP1+(2D0*ELP2+3D0*ELP3*T)*T) 1306 1307* Sun's mean anomaly 1308 EM=D2R*MOD(EM0+(EM1+(EM2+EM3*T)*T)*T,360D0) 1309 DEM=D2R*(EM1+(2D0*EM2+3D0*EM3*T)*T) 1310 1311* Moon's mean anomaly 1312 EMP=D2R*MOD(EMP0+(EMP1+(EMP2+EMP3*T)*T)*T,360D0) 1313 DEMP=D2R*(EMP1+(2D0*EMP2+3D0*EMP3*T)*T) 1314 1315* Moon's mean elongation 1316 D=D2R*MOD(D0+(D1+(D2+D3*T)*T)*T,360D0) 1317 DD=D2R*(D1+(2D0*D2+3D0*D3*T)*T) 1318 1319* Mean distance of the Moon from its ascending node 1320 F=D2R*MOD(F0+(F1+(F2+F3*T)*T)*T,360D0) 1321 DF=D2R*(F1+(2D0*F2+3D0*F3*T)*T) 1322 1323* Longitude of the Moon's ascending node 1324 OM=D2R*MOD(OM0+(OM1+(OM2+OM3*T)*T)*T,360D0) 1325 DOM=D2R*(OM1+(2D0*OM2+3D0*OM3*T)*T) 1326 SINOM=SIN(OM) 1327 COSOM=COS(OM) 1328 DOMCOM=DOM*COSOM 1329 1330* Add the periodic variations 1331 THETA=D2R*(PA0+PA1*T) 1332 WA=SIN(THETA) 1333 DWA=D2R*PA1*COS(THETA) 1334 THETA=D2R*(PE0+(PE1+PE2*T)*T) 1335 WB=PEC*SIN(THETA) 1336 DWB=D2R*PEC*(PE1+2D0*PE2*T)*COS(THETA) 1337 ELP=ELP+D2R*(PAC*WA+WB+PFC*SINOM) 1338 DELP=DELP+D2R*(PAC*DWA+DWB+PFC*DOMCOM) 1339 EM=EM+D2R*PBC*WA 1340 DEM=DEM+D2R*PBC*DWA 1341 EMP=EMP+D2R*(PCC*WA+WB+PGC*SINOM) 1342 DEMP=DEMP+D2R*(PCC*DWA+DWB+PGC*DOMCOM) 1343 D=D+D2R*(PDC*WA+WB+PHC*SINOM) 1344 DD=DD+D2R*(PDC*DWA+DWB+PHC*DOMCOM) 1345 WOM=OM+D2R*(PJ0+PJ1*T) 1346 DWOM=DOM+D2R*PJ1 1347 SINWOM=SIN(WOM) 1348 COSWOM=COS(WOM) 1349 F=F+D2R*(WB+PIC*SINOM+PJC*SINWOM) 1350 DF=DF+D2R*(DWB+PIC*DOMCOM+PJC*DWOM*COSWOM) 1351 1352* E-factor, and square 1353 E=1D0+(E1+E2*T)*T 1354 DE=E1+2D0*E2*T 1355 ESQ=E*E 1356 DESQ=2D0*E*DE 1357 1358* 1359* Series expansions 1360* 1361 1362* Longitude 1363 V=0D0 1364 DV=0D0 1365 DO N=NL,1,-1 1366 COEFF=TL(N) 1367 EMN=DBLE(ITL(1,N)) 1368 EMPN=DBLE(ITL(2,N)) 1369 DN=DBLE(ITL(3,N)) 1370 FN=DBLE(ITL(4,N)) 1371 I=ITL(5,N) 1372 IF (I.EQ.0) THEN 1373 EN=1D0 1374 DEN=0D0 1375 ELSE IF (I.EQ.1) THEN 1376 EN=E 1377 DEN=DE 1378 ELSE 1379 EN=ESQ 1380 DEN=DESQ 1381 END IF 1382 THETA=EMN*EM+EMPN*EMP+DN*D+FN*F 1383 DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF 1384 FTHETA=SIN(THETA) 1385 V=V+COEFF*FTHETA*EN 1386 DV=DV+COEFF*(COS(THETA)*DTHETA*EN+FTHETA*DEN) 1387 END DO 1388 EL=ELP+D2R*V 1389 DEL=(DELP+D2R*DV)/CJ 1390 1391* Latitude 1392 V=0D0 1393 DV=0D0 1394 DO N=NB,1,-1 1395 COEFF=TB(N) 1396 EMN=DBLE(ITB(1,N)) 1397 EMPN=DBLE(ITB(2,N)) 1398 DN=DBLE(ITB(3,N)) 1399 FN=DBLE(ITB(4,N)) 1400 I=ITB(5,N) 1401 IF (I.EQ.0) THEN 1402 EN=1D0 1403 DEN=0D0 1404 ELSE IF (I.EQ.1) THEN 1405 EN=E 1406 DEN=DE 1407 ELSE 1408 EN=ESQ 1409 DEN=DESQ 1410 END IF 1411 THETA=EMN*EM+EMPN*EMP+DN*D+FN*F 1412 DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF 1413 FTHETA=SIN(THETA) 1414 V=V+COEFF*FTHETA*EN 1415 DV=DV+COEFF*(COS(THETA)*DTHETA*EN+FTHETA*DEN) 1416 END DO 1417 BF=1D0-CW1*COSOM-CW2*COSWOM 1418 DBF=CW1*DOM*SINOM+CW2*DWOM*SINWOM 1419 B=D2R*V*BF 1420 DB=D2R*(DV*BF+V*DBF)/CJ 1421 1422* Parallax 1423 V=0D0 1424 DV=0D0 1425 DO N=NP,1,-1 1426 COEFF=TP(N) 1427 EMN=DBLE(ITP(1,N)) 1428 EMPN=DBLE(ITP(2,N)) 1429 DN=DBLE(ITP(3,N)) 1430 FN=DBLE(ITP(4,N)) 1431 I=ITP(5,N) 1432 IF (I.EQ.0) THEN 1433 EN=1D0 1434 DEN=0D0 1435 ELSE IF (I.EQ.1) THEN 1436 EN=E 1437 DEN=DE 1438 ELSE 1439 EN=ESQ 1440 DEN=DESQ 1441 END IF 1442 THETA=EMN*EM+EMPN*EMP+DN*D+FN*F 1443 DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF 1444 FTHETA=COS(THETA) 1445 V=V+COEFF*FTHETA*EN 1446 DV=DV+COEFF*(-SIN(THETA)*DTHETA*EN+FTHETA*DEN) 1447 END DO 1448 P=D2R*V 1449 DP=D2R*DV/CJ 1450 1451* 1452* Transformation into final form 1453* 1454 1455* Parallax to distance (AU, AU/sec) 1456 SP=SIN(P) 1457 R=ERADAU/SP 1458 DR=-R*DP*COS(P)/SP 1459 1460* Longitude, latitude to x,y,z (AU) 1461 SEL=SIN(EL) 1462 CEL=COS(EL) 1463 SB=SIN(B) 1464 CB=COS(B) 1465 RCB=R*CB 1466 RBD=R*DB 1467 W=RBD*SB-CB*DR 1468 X=RCB*CEL 1469 Y=RCB*SEL 1470 Z=R*SB 1471 XD=-Y*DEL-W*CEL 1472 YD=X*DEL-W*SEL 1473 ZD=RBD*CB+SB*DR 1474 1475* Julian centuries since J2000 1476 T=(DATE-51544.5D0)/36525D0 1477 1478* Fricke equinox correction 1479 EPJ=2000D0+T*100D0 1480 EQCOR=DS2R*(0.035D0+0.00085D0*(EPJ-B1950)) 1481 1482* Mean obliquity (IAU 1976) 1483 EPS=DAS2R*(84381.448D0+(-46.8150D0+(-0.00059D0+0.001813D0*T)*T)*T) 1484 1485* To the equatorial system, mean of date, FK5 system 1486 SINEPS=SIN(EPS) 1487 COSEPS=COS(EPS) 1488 ES=EQCOR*SINEPS 1489 EC=EQCOR*COSEPS 1490 PV(1)=X-EC*Y+ES*Z 1491 PV(2)=EQCOR*X+Y*COSEPS-Z*SINEPS 1492 PV(3)=Y*SINEPS+Z*COSEPS 1493 PV(4)=XD-EC*YD+ES*ZD 1494 PV(5)=EQCOR*XD+YD*COSEPS-ZD*SINEPS 1495 PV(6)=YD*SINEPS+ZD*COSEPS 1496 1497 END 1498 SUBROUTINE sla_DMXV (DM, VA, VB) 1499*+ 1500* - - - - - 1501* D M X V 1502* - - - - - 1503* 1504* Performs the 3-D forward unitary transformation: 1505* 1506* vector VB = matrix DM * vector VA 1507* 1508* (double precision) 1509* 1510* Given: 1511* DM dp(3,3) matrix 1512* VA dp(3) vector 1513* 1514* Returned: 1515* VB dp(3) result vector 1516* 1517* To comply with the ANSI Fortran 77 standard, VA and VB must be 1518* different arrays. However, the routine is coded so as to work 1519* properly on many platforms even if this rule is violated. 1520* 1521* Last revision: 26 December 2004 1522* 1523* Copyright P.T.Wallace. All rights reserved. 1524* 1525* License: 1526* This program is free software; you can redistribute it and/or modify 1527* it under the terms of the GNU General Public License as published by 1528* the Free Software Foundation; either version 2 of the License, or 1529* (at your option) any later version. 1530* 1531* This program is distributed in the hope that it will be useful, 1532* but WITHOUT ANY WARRANTY; without even the implied warranty of 1533* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1534* GNU General Public License for more details. 1535* 1536* You should have received a copy of the GNU General Public License 1537* along with this program (see SLA_CONDITIONS); if not, write to the 1538* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 1539* Boston, MA 02110-1301 USA 1540* 1541*- 1542 1543 IMPLICIT NONE 1544 1545 DOUBLE PRECISION DM(3,3),VA(3),VB(3) 1546 1547 INTEGER I,J 1548 DOUBLE PRECISION W,VW(3) 1549 1550 1551* Matrix DM * vector VA -> vector VW 1552 DO J=1,3 1553 W=0D0 1554 DO I=1,3 1555 W=W+DM(J,I)*VA(I) 1556 END DO 1557 VW(J)=W 1558 END DO 1559 1560* Vector VW -> vector VB 1561 DO J=1,3 1562 VB(J)=VW(J) 1563 END DO 1564 1565 END 1566 DOUBLE PRECISION FUNCTION sla_DRANGE (ANGLE) 1567*+ 1568* - - - - - - - 1569* D R A N G E 1570* - - - - - - - 1571* 1572* Normalize angle into range +/- pi (double precision) 1573* 1574* Given: 1575* ANGLE dp the angle in radians 1576* 1577* The result (double precision) is ANGLE expressed in the range +/- pi. 1578* 1579* P.T.Wallace Starlink 23 November 1995 1580* 1581* Copyright (C) 1995 Rutherford Appleton Laboratory 1582* 1583* License: 1584* This program is free software; you can redistribute it and/or modify 1585* it under the terms of the GNU General Public License as published by 1586* the Free Software Foundation; either version 2 of the License, or 1587* (at your option) any later version. 1588* 1589* This program is distributed in the hope that it will be useful, 1590* but WITHOUT ANY WARRANTY; without even the implied warranty of 1591* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1592* GNU General Public License for more details. 1593* 1594* You should have received a copy of the GNU General Public License 1595* along with this program (see SLA_CONDITIONS); if not, write to the 1596* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 1597* Boston, MA 02110-1301 USA 1598* 1599*- 1600 1601 IMPLICIT NONE 1602 1603 DOUBLE PRECISION ANGLE 1604 1605 DOUBLE PRECISION DPI,D2PI 1606 PARAMETER (DPI=3.141592653589793238462643D0) 1607 PARAMETER (D2PI=6.283185307179586476925287D0) 1608 1609 1610 sla_DRANGE=MOD(ANGLE,D2PI) 1611 IF (ABS(sla_DRANGE).GE.DPI) 1612 : sla_DRANGE=sla_DRANGE-SIGN(D2PI,ANGLE) 1613 1614 END 1615 DOUBLE PRECISION FUNCTION sla_DRANRM (ANGLE) 1616*+ 1617* - - - - - - - 1618* D R A N R M 1619* - - - - - - - 1620* 1621* Normalize angle into range 0-2 pi (double precision) 1622* 1623* Given: 1624* ANGLE dp the angle in radians 1625* 1626* The result is ANGLE expressed in the range 0-2 pi. 1627* 1628* Last revision: 22 July 2004 1629* 1630* Copyright P.T.Wallace. All rights reserved. 1631* 1632* License: 1633* This program is free software; you can redistribute it and/or modify 1634* it under the terms of the GNU General Public License as published by 1635* the Free Software Foundation; either version 2 of the License, or 1636* (at your option) any later version. 1637* 1638* This program is distributed in the hope that it will be useful, 1639* but WITHOUT ANY WARRANTY; without even the implied warranty of 1640* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1641* GNU General Public License for more details. 1642* 1643* You should have received a copy of the GNU General Public License 1644* along with this program (see SLA_CONDITIONS); if not, write to the 1645* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 1646* Boston, MA 02110-1301 USA 1647* 1648*- 1649 1650 IMPLICIT NONE 1651 1652 DOUBLE PRECISION ANGLE 1653 1654 DOUBLE PRECISION D2PI 1655 PARAMETER (D2PI=6.283185307179586476925286766559D0) 1656 1657 1658 sla_DRANRM = MOD(ANGLE,D2PI) 1659 IF (sla_DRANRM.LT.0D0) sla_DRANRM = sla_DRANRM+D2PI 1660 1661 END 1662 DOUBLE PRECISION FUNCTION sla_DTT (UTC) 1663*+ 1664* - - - - 1665* D T T 1666* - - - - 1667* 1668* Increment to be applied to Coordinated Universal Time UTC to give 1669* Terrestrial Time TT (formerly Ephemeris Time ET) 1670* 1671* (double precision) 1672* 1673* Given: 1674* UTC d UTC date as a modified JD (JD-2400000.5) 1675* 1676* Result: TT-UTC in seconds 1677* 1678* Notes: 1679* 1680* 1 The UTC is specified to be a date rather than a time to indicate 1681* that care needs to be taken not to specify an instant which lies 1682* within a leap second. Though in most cases UTC can include the 1683* fractional part, correct behaviour on the day of a leap second 1684* can only be guaranteed up to the end of the second 23:59:59. 1685* 1686* 2 Pre 1972 January 1 a fixed value of 10 + ET-TAI is returned. 1687* 1688* 3 See also the routine sla_DT, which roughly estimates ET-UT for 1689* historical epochs. 1690* 1691* Called: sla_DAT 1692* 1693* P.T.Wallace Starlink 6 December 1994 1694* 1695* Copyright (C) 1995 Rutherford Appleton Laboratory 1696* 1697* License: 1698* This program is free software; you can redistribute it and/or modify 1699* it under the terms of the GNU General Public License as published by 1700* the Free Software Foundation; either version 2 of the License, or 1701* (at your option) any later version. 1702* 1703* This program is distributed in the hope that it will be useful, 1704* but WITHOUT ANY WARRANTY; without even the implied warranty of 1705* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1706* GNU General Public License for more details. 1707* 1708* You should have received a copy of the GNU General Public License 1709* along with this program (see SLA_CONDITIONS); if not, write to the 1710* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 1711* Boston, MA 02110-1301 USA 1712* 1713*- 1714 1715 IMPLICIT NONE 1716 1717 DOUBLE PRECISION UTC 1718 1719 DOUBLE PRECISION sla_DAT 1720 1721 1722 sla_DTT=32.184D0+sla_DAT(UTC) 1723 1724 END 1725 SUBROUTINE sla_ECMAT (DATE, RMAT) 1726*+ 1727* - - - - - - 1728* E C M A T 1729* - - - - - - 1730* 1731* Form the equatorial to ecliptic rotation matrix - IAU 1980 theory 1732* (double precision) 1733* 1734* Given: 1735* DATE dp TDB (loosely ET) as Modified Julian Date 1736* (JD-2400000.5) 1737* Returned: 1738* RMAT dp(3,3) matrix 1739* 1740* Reference: 1741* Murray,C.A., Vectorial Astrometry, section 4.3. 1742* 1743* Note: 1744* The matrix is in the sense V(ecl) = RMAT * V(equ); the 1745* equator, equinox and ecliptic are mean of date. 1746* 1747* Called: sla_DEULER 1748* 1749* P.T.Wallace Starlink 23 August 1996 1750* 1751* Copyright (C) 1996 Rutherford Appleton Laboratory 1752* 1753* License: 1754* This program is free software; you can redistribute it and/or modify 1755* it under the terms of the GNU General Public License as published by 1756* the Free Software Foundation; either version 2 of the License, or 1757* (at your option) any later version. 1758* 1759* This program is distributed in the hope that it will be useful, 1760* but WITHOUT ANY WARRANTY; without even the implied warranty of 1761* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1762* GNU General Public License for more details. 1763* 1764* You should have received a copy of the GNU General Public License 1765* along with this program (see SLA_CONDITIONS); if not, write to the 1766* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 1767* Boston, MA 02110-1301 USA 1768* 1769*- 1770 1771 IMPLICIT NONE 1772 1773 DOUBLE PRECISION DATE,RMAT(3,3) 1774 1775* Arc seconds to radians 1776 DOUBLE PRECISION AS2R 1777 PARAMETER (AS2R=0.484813681109535994D-5) 1778 1779 DOUBLE PRECISION T,EPS0 1780 1781 1782 1783* Interval between basic epoch J2000.0 and current epoch (JC) 1784 T = (DATE-51544.5D0)/36525D0 1785 1786* Mean obliquity 1787 EPS0 = AS2R* 1788 : (84381.448D0+(-46.8150D0+(-0.00059D0+0.001813D0*T)*T)*T) 1789 1790* Matrix 1791 CALL sla_DEULER('X',EPS0,0D0,0D0,RMAT) 1792 1793 END 1794 DOUBLE PRECISION FUNCTION sla_EPJ (DATE) 1795*+ 1796* - - - - 1797* E P J 1798* - - - - 1799* 1800* Conversion of Modified Julian Date to Julian Epoch (double precision) 1801* 1802* Given: 1803* DATE dp Modified Julian Date (JD - 2400000.5) 1804* 1805* The result is the Julian Epoch. 1806* 1807* Reference: 1808* Lieske,J.H., 1979. Astron.Astrophys.,73,282. 1809* 1810* P.T.Wallace Starlink February 1984 1811* 1812* Copyright (C) 1995 Rutherford Appleton Laboratory 1813* 1814* License: 1815* This program is free software; you can redistribute it and/or modify 1816* it under the terms of the GNU General Public License as published by 1817* the Free Software Foundation; either version 2 of the License, or 1818* (at your option) any later version. 1819* 1820* This program is distributed in the hope that it will be useful, 1821* but WITHOUT ANY WARRANTY; without even the implied warranty of 1822* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1823* GNU General Public License for more details. 1824* 1825* You should have received a copy of the GNU General Public License 1826* along with this program (see SLA_CONDITIONS); if not, write to the 1827* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 1828* Boston, MA 02110-1301 USA 1829* 1830*- 1831 1832 IMPLICIT NONE 1833 1834 DOUBLE PRECISION DATE 1835 1836 1837 sla_EPJ = 2000D0 + (DATE-51544.5D0)/365.25D0 1838 1839 END 1840 SUBROUTINE sla_EQECL (DR, DD, DATE, DL, DB) 1841*+ 1842* - - - - - - 1843* E Q E C L 1844* - - - - - - 1845* 1846* Transformation from J2000.0 equatorial coordinates to 1847* ecliptic coordinates (double precision) 1848* 1849* Given: 1850* DR,DD dp J2000.0 mean RA,Dec (radians) 1851* DATE dp TDB (loosely ET) as Modified Julian Date 1852* (JD-2400000.5) 1853* Returned: 1854* DL,DB dp ecliptic longitude and latitude 1855* (mean of date, IAU 1980 theory, radians) 1856* 1857* Called: 1858* sla_DCS2C, sla_PREC, sla_EPJ, sla_DMXV, sla_ECMAT, sla_DCC2S, 1859* sla_DRANRM, sla_DRANGE 1860* 1861* P.T.Wallace Starlink March 1986 1862* 1863* Copyright (C) 1995 Rutherford Appleton Laboratory 1864* 1865* License: 1866* This program is free software; you can redistribute it and/or modify 1867* it under the terms of the GNU General Public License as published by 1868* the Free Software Foundation; either version 2 of the License, or 1869* (at your option) any later version. 1870* 1871* This program is distributed in the hope that it will be useful, 1872* but WITHOUT ANY WARRANTY; without even the implied warranty of 1873* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1874* GNU General Public License for more details. 1875* 1876* You should have received a copy of the GNU General Public License 1877* along with this program (see SLA_CONDITIONS); if not, write to the 1878* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 1879* Boston, MA 02110-1301 USA 1880* 1881*- 1882 1883 IMPLICIT NONE 1884 1885 DOUBLE PRECISION DR,DD,DATE,DL,DB 1886 1887 DOUBLE PRECISION sla_EPJ,sla_DRANRM,sla_DRANGE 1888 1889 DOUBLE PRECISION RMAT(3,3),V1(3),V2(3) 1890 1891 1892 1893* Spherical to Cartesian 1894 CALL sla_DCS2C(DR,DD,V1) 1895 1896* Mean J2000 to mean of date 1897 CALL sla_PREC(2000D0,sla_EPJ(DATE),RMAT) 1898 CALL sla_DMXV(RMAT,V1,V2) 1899 1900* Equatorial to ecliptic 1901 CALL sla_ECMAT(DATE,RMAT) 1902 CALL sla_DMXV(RMAT,V2,V1) 1903 1904* Cartesian to spherical 1905 CALL sla_DCC2S(V1,DL,DB) 1906 1907* Express in conventional ranges 1908 DL=sla_DRANRM(DL) 1909 DB=sla_DRANGE(DB) 1910 1911 END 1912 DOUBLE PRECISION FUNCTION sla_EQEQX (DATE) 1913*+ 1914* - - - - - - 1915* E Q E Q X 1916* - - - - - - 1917* 1918* Equation of the equinoxes (IAU 1994, double precision) 1919* 1920* Given: 1921* DATE dp TDB (loosely ET) as Modified Julian Date 1922* (JD-2400000.5) 1923* 1924* The result is the equation of the equinoxes (double precision) 1925* in radians: 1926* 1927* Greenwich apparent ST = GMST + sla_EQEQX 1928* 1929* References: IAU Resolution C7, Recommendation 3 (1994) 1930* Capitaine, N. & Gontier, A.-M., Astron. Astrophys., 1931* 275, 645-650 (1993) 1932* 1933* Called: sla_NUTC 1934* 1935* Patrick Wallace Starlink 23 August 1996 1936* 1937* Copyright (C) 1996 Rutherford Appleton Laboratory 1938* 1939* License: 1940* This program is free software; you can redistribute it and/or modify 1941* it under the terms of the GNU General Public License as published by 1942* the Free Software Foundation; either version 2 of the License, or 1943* (at your option) any later version. 1944* 1945* This program is distributed in the hope that it will be useful, 1946* but WITHOUT ANY WARRANTY; without even the implied warranty of 1947* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1948* GNU General Public License for more details. 1949* 1950* You should have received a copy of the GNU General Public License 1951* along with this program (see SLA_CONDITIONS); if not, write to the 1952* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 1953* Boston, MA 02110-1301 USA 1954* 1955*- 1956 1957 IMPLICIT NONE 1958 1959 DOUBLE PRECISION DATE 1960 1961* Turns to arc seconds and arc seconds to radians 1962 DOUBLE PRECISION T2AS,AS2R 1963 PARAMETER (T2AS=1296000D0, 1964 : AS2R=0.484813681109535994D-5) 1965 1966 DOUBLE PRECISION T,OM,DPSI,DEPS,EPS0 1967 1968 1969 1970* Interval between basic epoch J2000.0 and current epoch (JC) 1971 T=(DATE-51544.5D0)/36525D0 1972 1973* Longitude of the mean ascending node of the lunar orbit on the 1974* ecliptic, measured from the mean equinox of date 1975 OM=AS2R*(450160.280D0+(-5D0*T2AS-482890.539D0 1976 : +(7.455D0+0.008D0*T)*T)*T) 1977 1978* Nutation 1979 CALL sla_NUTC(DATE,DPSI,DEPS,EPS0) 1980 1981* Equation of the equinoxes 1982 sla_EQEQX=DPSI*COS(EPS0)+AS2R*(0.00264D0*SIN(OM)+ 1983 : 0.000063D0*SIN(OM+OM)) 1984 1985 END 1986 SUBROUTINE sla_GEOC (P, H, R, Z) 1987*+ 1988* - - - - - 1989* G E O C 1990* - - - - - 1991* 1992* Convert geodetic position to geocentric (double precision) 1993* 1994* Given: 1995* P dp latitude (geodetic, radians) 1996* H dp height above reference spheroid (geodetic, metres) 1997* 1998* Returned: 1999* R dp distance from Earth axis (AU) 2000* Z dp distance from plane of Earth equator (AU) 2001* 2002* Notes: 2003* 2004* 1 Geocentric latitude can be obtained by evaluating ATAN2(Z,R). 2005* 2006* 2 IAU 1976 constants are used. 2007* 2008* Reference: 2009* 2010* Green,R.M., Spherical Astronomy, CUP 1985, p98. 2011* 2012* Last revision: 22 July 2004 2013* 2014* Copyright P.T.Wallace. All rights reserved. 2015* 2016* License: 2017* This program is free software; you can redistribute it and/or modify 2018* it under the terms of the GNU General Public License as published by 2019* the Free Software Foundation; either version 2 of the License, or 2020* (at your option) any later version. 2021* 2022* This program is distributed in the hope that it will be useful, 2023* but WITHOUT ANY WARRANTY; without even the implied warranty of 2024* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 2025* GNU General Public License for more details. 2026* 2027* You should have received a copy of the GNU General Public License 2028* along with this program (see SLA_CONDITIONS); if not, write to the 2029* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 2030* Boston, MA 02110-1301 USA 2031* 2032*- 2033 2034 IMPLICIT NONE 2035 2036 DOUBLE PRECISION P,H,R,Z 2037 2038* Earth equatorial radius (metres) 2039 DOUBLE PRECISION A0 2040 PARAMETER (A0=6378140D0) 2041 2042* Reference spheroid flattening factor and useful function 2043 DOUBLE PRECISION F,B 2044 PARAMETER (F=1D0/298.257D0,B=(1D0-F)**2) 2045 2046* Astronomical unit in metres 2047 DOUBLE PRECISION AU 2048 PARAMETER (AU=1.49597870D11) 2049 2050 DOUBLE PRECISION SP,CP,C,S 2051 2052 2053 2054* Geodetic to geocentric conversion 2055 SP = SIN(P) 2056 CP = COS(P) 2057 C = 1D0/SQRT(CP*CP+B*SP*SP) 2058 S = B*C 2059 R = (A0*C+H)*CP/AU 2060 Z = (A0*S+H)*SP/AU 2061 2062 END 2063 DOUBLE PRECISION FUNCTION sla_GMST (UT1) 2064*+ 2065* - - - - - 2066* G M S T 2067* - - - - - 2068* 2069* Conversion from universal time to sidereal time (double precision) 2070* 2071* Given: 2072* UT1 dp universal time (strictly UT1) expressed as 2073* modified Julian Date (JD-2400000.5) 2074* 2075* The result is the Greenwich mean sidereal time (double 2076* precision, radians). 2077* 2078* The IAU 1982 expression (see page S15 of 1984 Astronomical Almanac) 2079* is used, but rearranged to reduce rounding errors. This expression 2080* is always described as giving the GMST at 0 hours UT. In fact, it 2081* gives the difference between the GMST and the UT, which happens to 2082* equal the GMST (modulo 24 hours) at 0 hours UT each day. In this 2083* routine, the entire UT is used directly as the argument for the 2084* standard formula, and the fractional part of the UT is added 2085* separately. Note that the factor 1.0027379... does not appear in the 2086* IAU 1982 expression explicitly but in the form of the coefficient 2087* 8640184.812866, which is 86400x36525x0.0027379... 2088* 2089* See also the routine sla_GMSTA, which delivers better numerical 2090* precision by accepting the UT date and time as separate arguments. 2091* 2092* Called: sla_DRANRM 2093* 2094* P.T.Wallace Starlink 14 October 2001 2095* 2096* Copyright (C) 2001 Rutherford Appleton Laboratory 2097* 2098* License: 2099* This program is free software; you can redistribute it and/or modify 2100* it under the terms of the GNU General Public License as published by 2101* the Free Software Foundation; either version 2 of the License, or 2102* (at your option) any later version. 2103* 2104* This program is distributed in the hope that it will be useful, 2105* but WITHOUT ANY WARRANTY; without even the implied warranty of 2106* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 2107* GNU General Public License for more details. 2108* 2109* You should have received a copy of the GNU General Public License 2110* along with this program (see SLA_CONDITIONS); if not, write to the 2111* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 2112* Boston, MA 02110-1301 USA 2113* 2114*- 2115 2116 IMPLICIT NONE 2117 2118 DOUBLE PRECISION UT1 2119 2120 DOUBLE PRECISION sla_DRANRM 2121 2122 DOUBLE PRECISION D2PI,S2R 2123 PARAMETER (D2PI=6.283185307179586476925286766559D0, 2124 : S2R=7.272205216643039903848711535369D-5) 2125 2126 DOUBLE PRECISION TU 2127 2128 2129 2130* Julian centuries from fundamental epoch J2000 to this UT 2131 TU=(UT1-51544.5D0)/36525D0 2132 2133* GMST at this UT 2134 sla_GMST=sla_DRANRM(MOD(UT1,1D0)*D2PI+ 2135 : (24110.54841D0+ 2136 : (8640184.812866D0+ 2137 : (0.093104D0-6.2D-6*TU)*TU)*TU)*S2R) 2138 2139 END 2140 SUBROUTINE sla_NUTC (DATE, DPSI, DEPS, EPS0) 2141*+ 2142* - - - - - 2143* N U T C 2144* - - - - - 2145* 2146* Nutation: longitude & obliquity components and mean obliquity, 2147* using the Shirai & Fukushima (2001) theory. 2148* 2149* Given: 2150* DATE d TDB (loosely ET) as Modified Julian Date 2151* (JD-2400000.5) 2152* Returned: 2153* DPSI,DEPS d nutation in longitude,obliquity 2154* EPS0 d mean obliquity 2155* 2156* Notes: 2157* 2158* 1 The routine predicts forced nutation (but not free core nutation) 2159* plus corrections to the IAU 1976 precession model. 2160* 2161* 2 Earth attitude predictions made by combining the present nutation 2162* model with IAU 1976 precession are accurate to 1 mas (with respect 2163* to the ICRF) for a few decades around 2000. 2164* 2165* 3 The sla_NUTC80 routine is the equivalent of the present routine 2166* but using the IAU 1980 nutation theory. The older theory is less 2167* accurate, leading to errors as large as 350 mas over the interval 2168* 1900-2100, mainly because of the error in the IAU 1976 precession. 2169* 2170* References: 2171* 2172* Shirai, T. & Fukushima, T., Astron.J. 121, 3270-3283 (2001). 2173* 2174* Fukushima, T., Astron.Astrophys. 244, L11 (1991). 2175* 2176* Simon, J. L., Bretagnon, P., Chapront, J., Chapront-Touze, M., 2177* Francou, G. & Laskar, J., Astron.Astrophys. 282, 663 (1994). 2178* 2179* This revision: 24 November 2005 2180* 2181* Copyright P.T.Wallace. All rights reserved. 2182* 2183* License: 2184* This program is free software; you can redistribute it and/or modify 2185* it under the terms of the GNU General Public License as published by 2186* the Free Software Foundation; either version 2 of the License, or 2187* (at your option) any later version. 2188* 2189* This program is distributed in the hope that it will be useful, 2190* but WITHOUT ANY WARRANTY; without even the implied warranty of 2191* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 2192* GNU General Public License for more details. 2193* 2194* You should have received a copy of the GNU General Public License 2195* along with this program (see SLA_CONDITIONS); if not, write to the 2196* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 2197* Boston, MA 02110-1301 USA 2198* 2199*- 2200 2201 IMPLICIT NONE 2202 2203 DOUBLE PRECISION DATE,DPSI,DEPS,EPS0 2204 2205* Degrees to radians 2206 DOUBLE PRECISION DD2R 2207 PARAMETER (DD2R=1.745329251994329576923691D-2) 2208 2209* Arc seconds to radians 2210 DOUBLE PRECISION DAS2R 2211 PARAMETER (DAS2R=4.848136811095359935899141D-6) 2212 2213* Arc seconds in a full circle 2214 DOUBLE PRECISION TURNAS 2215 PARAMETER (TURNAS=1296000D0) 2216 2217* Reference epoch (J2000), MJD 2218 DOUBLE PRECISION DJM0 2219 PARAMETER (DJM0=51544.5D0 ) 2220 2221* Days per Julian century 2222 DOUBLE PRECISION DJC 2223 PARAMETER (DJC=36525D0) 2224 2225 INTEGER I,J 2226 DOUBLE PRECISION T,EL,ELP,F,D,OM,VE,MA,JU,SA,THETA,C,S,DP,DE 2227 2228* Number of terms in the nutation model 2229 INTEGER NTERMS 2230 PARAMETER (NTERMS=194) 2231 2232* The SF2001 forced nutation model 2233 INTEGER NA(9,NTERMS) 2234 DOUBLE PRECISION PSI(4,NTERMS), EPS(4,NTERMS) 2235 2236* Coefficients of fundamental angles 2237 DATA ( ( NA(I,J), I=1,9 ), J=1,10 ) / 2238 : 0, 0, 0, 0, -1, 0, 0, 0, 0, 2239 : 0, 0, 2, -2, 2, 0, 0, 0, 0, 2240 : 0, 0, 2, 0, 2, 0, 0, 0, 0, 2241 : 0, 0, 0, 0, -2, 0, 0, 0, 0, 2242 : 0, 1, 0, 0, 0, 0, 0, 0, 0, 2243 : 0, 1, 2, -2, 2, 0, 0, 0, 0, 2244 : 1, 0, 0, 0, 0, 0, 0, 0, 0, 2245 : 0, 0, 2, 0, 1, 0, 0, 0, 0, 2246 : 1, 0, 2, 0, 2, 0, 0, 0, 0, 2247 : 0, -1, 2, -2, 2, 0, 0, 0, 0 / 2248 DATA ( ( NA(I,J), I=1,9 ), J=11,20 ) / 2249 : 0, 0, 2, -2, 1, 0, 0, 0, 0, 2250 : -1, 0, 2, 0, 2, 0, 0, 0, 0, 2251 : -1, 0, 0, 2, 0, 0, 0, 0, 0, 2252 : 1, 0, 0, 0, 1, 0, 0, 0, 0, 2253 : 1, 0, 0, 0, -1, 0, 0, 0, 0, 2254 : -1, 0, 2, 2, 2, 0, 0, 0, 0, 2255 : 1, 0, 2, 0, 1, 0, 0, 0, 0, 2256 : -2, 0, 2, 0, 1, 0, 0, 0, 0, 2257 : 0, 0, 0, 2, 0, 0, 0, 0, 0, 2258 : 0, 0, 2, 2, 2, 0, 0, 0, 0 / 2259 DATA ( ( NA(I,J), I=1,9 ), J=21,30 ) / 2260 : 2, 0, 0, -2, 0, 0, 0, 0, 0, 2261 : 2, 0, 2, 0, 2, 0, 0, 0, 0, 2262 : 1, 0, 2, -2, 2, 0, 0, 0, 0, 2263 : -1, 0, 2, 0, 1, 0, 0, 0, 0, 2264 : 2, 0, 0, 0, 0, 0, 0, 0, 0, 2265 : 0, 0, 2, 0, 0, 0, 0, 0, 0, 2266 : 0, 1, 0, 0, 1, 0, 0, 0, 0, 2267 : -1, 0, 0, 2, 1, 0, 0, 0, 0, 2268 : 0, 2, 2, -2, 2, 0, 0, 0, 0, 2269 : 0, 0, 2, -2, 0, 0, 0, 0, 0 / 2270 DATA ( ( NA(I,J), I=1,9 ), J=31,40 ) / 2271 : -1, 0, 0, 2, -1, 0, 0, 0, 0, 2272 : 0, 1, 0, 0, -1, 0, 0, 0, 0, 2273 : 0, 2, 0, 0, 0, 0, 0, 0, 0, 2274 : -1, 0, 2, 2, 1, 0, 0, 0, 0, 2275 : 1, 0, 2, 2, 2, 0, 0, 0, 0, 2276 : 0, 1, 2, 0, 2, 0, 0, 0, 0, 2277 : -2, 0, 2, 0, 0, 0, 0, 0, 0, 2278 : 0, 0, 2, 2, 1, 0, 0, 0, 0, 2279 : 0, -1, 2, 0, 2, 0, 0, 0, 0, 2280 : 0, 0, 0, 2, 1, 0, 0, 0, 0 / 2281 DATA ( ( NA(I,J), I=1,9 ), J=41,50 ) / 2282 : 1, 0, 2, -2, 1, 0, 0, 0, 0, 2283 : 2, 0, 0, -2, -1, 0, 0, 0, 0, 2284 : 2, 0, 2, -2, 2, 0, 0, 0, 0, 2285 : 2, 0, 2, 0, 1, 0, 0, 0, 0, 2286 : 0, 0, 0, 2, -1, 0, 0, 0, 0, 2287 : 0, -1, 2, -2, 1, 0, 0, 0, 0, 2288 : -1, -1, 0, 2, 0, 0, 0, 0, 0, 2289 : 2, 0, 0, -2, 1, 0, 0, 0, 0, 2290 : 1, 0, 0, 2, 0, 0, 0, 0, 0, 2291 : 0, 1, 2, -2, 1, 0, 0, 0, 0 / 2292 DATA ( ( NA(I,J), I=1,9 ), J=51,60 ) / 2293 : 1, -1, 0, 0, 0, 0, 0, 0, 0, 2294 : -2, 0, 2, 0, 2, 0, 0, 0, 0, 2295 : 0, -1, 0, 2, 0, 0, 0, 0, 0, 2296 : 3, 0, 2, 0, 2, 0, 0, 0, 0, 2297 : 0, 0, 0, 1, 0, 0, 0, 0, 0, 2298 : 1, -1, 2, 0, 2, 0, 0, 0, 0, 2299 : 1, 0, 0, -1, 0, 0, 0, 0, 0, 2300 : -1, -1, 2, 2, 2, 0, 0, 0, 0, 2301 : -1, 0, 2, 0, 0, 0, 0, 0, 0, 2302 : 2, 0, 0, 0, -1, 0, 0, 0, 0 / 2303 DATA ( ( NA(I,J), I=1,9 ), J=61,70 ) / 2304 : 0, -1, 2, 2, 2, 0, 0, 0, 0, 2305 : 1, 1, 2, 0, 2, 0, 0, 0, 0, 2306 : 2, 0, 0, 0, 1, 0, 0, 0, 0, 2307 : 1, 1, 0, 0, 0, 0, 0, 0, 0, 2308 : 1, 0, -2, 2, -1, 0, 0, 0, 0, 2309 : 1, 0, 2, 0, 0, 0, 0, 0, 0, 2310 : -1, 1, 0, 1, 0, 0, 0, 0, 0, 2311 : 1, 0, 0, 0, 2, 0, 0, 0, 0, 2312 : -1, 0, 1, 0, 1, 0, 0, 0, 0, 2313 : 0, 0, 2, 1, 2, 0, 0, 0, 0 / 2314 DATA ( ( NA(I,J), I=1,9 ), J=71,80 ) / 2315 : -1, 1, 0, 1, 1, 0, 0, 0, 0, 2316 : -1, 0, 2, 4, 2, 0, 0, 0, 0, 2317 : 0, -2, 2, -2, 1, 0, 0, 0, 0, 2318 : 1, 0, 2, 2, 1, 0, 0, 0, 0, 2319 : 1, 0, 0, 0, -2, 0, 0, 0, 0, 2320 : -2, 0, 2, 2, 2, 0, 0, 0, 0, 2321 : 1, 1, 2, -2, 2, 0, 0, 0, 0, 2322 : -2, 0, 2, 4, 2, 0, 0, 0, 0, 2323 : -1, 0, 4, 0, 2, 0, 0, 0, 0, 2324 : 2, 0, 2, -2, 1, 0, 0, 0, 0 / 2325 DATA ( ( NA(I,J), I=1,9 ), J=81,90 ) / 2326 : 1, 0, 0, -1, -1, 0, 0, 0, 0, 2327 : 2, 0, 2, 2, 2, 0, 0, 0, 0, 2328 : 1, 0, 0, 2, 1, 0, 0, 0, 0, 2329 : 3, 0, 0, 0, 0, 0, 0, 0, 0, 2330 : 0, 0, 2, -2, -1, 0, 0, 0, 0, 2331 : 3, 0, 2, -2, 2, 0, 0, 0, 0, 2332 : 0, 0, 4, -2, 2, 0, 0, 0, 0, 2333 : -1, 0, 0, 4, 0, 0, 0, 0, 0, 2334 : 0, 1, 2, 0, 1, 0, 0, 0, 0, 2335 : 0, 0, 2, -2, 3, 0, 0, 0, 0 / 2336 DATA ( ( NA(I,J), I=1,9 ), J=91,100 ) / 2337 : -2, 0, 0, 4, 0, 0, 0, 0, 0, 2338 : -1, -1, 0, 2, 1, 0, 0, 0, 0, 2339 : -2, 0, 2, 0, -1, 0, 0, 0, 0, 2340 : 0, 0, 2, 0, -1, 0, 0, 0, 0, 2341 : 0, -1, 2, 0, 1, 0, 0, 0, 0, 2342 : 0, 1, 0, 0, 2, 0, 0, 0, 0, 2343 : 0, 0, 2, -1, 2, 0, 0, 0, 0, 2344 : 2, 1, 0, -2, 0, 0, 0, 0, 0, 2345 : 0, 0, 2, 4, 2, 0, 0, 0, 0, 2346 : -1, -1, 0, 2, -1, 0, 0, 0, 0 / 2347 DATA ( ( NA(I,J), I=1,9 ), J=101,110 ) / 2348 : -1, 1, 0, 2, 0, 0, 0, 0, 0, 2349 : 1, -1, 0, 0, 1, 0, 0, 0, 0, 2350 : 0, -1, 2, -2, 0, 0, 0, 0, 0, 2351 : 0, 1, 0, 0, -2, 0, 0, 0, 0, 2352 : 1, -1, 2, 2, 2, 0, 0, 0, 0, 2353 : 1, 0, 0, 2, -1, 0, 0, 0, 0, 2354 : -1, 1, 2, 2, 2, 0, 0, 0, 0, 2355 : 3, 0, 2, 0, 1, 0, 0, 0, 0, 2356 : 0, 1, 2, 2, 2, 0, 0, 0, 0, 2357 : 1, 0, 2, -2, 0, 0, 0, 0, 0 / 2358 DATA ( ( NA(I,J), I=1,9 ), J=111,120 ) / 2359 : -1, 0, -2, 4, -1, 0, 0, 0, 0, 2360 : -1, -1, 2, 2, 1, 0, 0, 0, 0, 2361 : 0, -1, 2, 2, 1, 0, 0, 0, 0, 2362 : 2, -1, 2, 0, 2, 0, 0, 0, 0, 2363 : 0, 0, 0, 2, 2, 0, 0, 0, 0, 2364 : 1, -1, 2, 0, 1, 0, 0, 0, 0, 2365 : -1, 1, 2, 0, 2, 0, 0, 0, 0, 2366 : 0, 1, 0, 2, 0, 0, 0, 0, 0, 2367 : 0, 1, 2, -2, 0, 0, 0, 0, 0, 2368 : 0, 3, 2, -2, 2, 0, 0, 0, 0 / 2369 DATA ( ( NA(I,J), I=1,9 ), J=121,130 ) / 2370 : 0, 0, 0, 1, 1, 0, 0, 0, 0, 2371 : -1, 0, 2, 2, 0, 0, 0, 0, 0, 2372 : 2, 1, 2, 0, 2, 0, 0, 0, 0, 2373 : 1, 1, 0, 0, 1, 0, 0, 0, 0, 2374 : 2, 0, 0, 2, 0, 0, 0, 0, 0, 2375 : 1, 1, 2, 0, 1, 0, 0, 0, 0, 2376 : -1, 0, 0, 2, 2, 0, 0, 0, 0, 2377 : 1, 0, -2, 2, 0, 0, 0, 0, 0, 2378 : 0, -1, 0, 2, -1, 0, 0, 0, 0, 2379 : -1, 0, 1, 0, 2, 0, 0, 0, 0 / 2380 DATA ( ( NA(I,J), I=1,9 ), J=131,140 ) / 2381 : 0, 1, 0, 1, 0, 0, 0, 0, 0, 2382 : 1, 0, -2, 2, -2, 0, 0, 0, 0, 2383 : 0, 0, 0, 1, -1, 0, 0, 0, 0, 2384 : 1, -1, 0, 0, -1, 0, 0, 0, 0, 2385 : 0, 0, 0, 4, 0, 0, 0, 0, 0, 2386 : 1, -1, 0, 2, 0, 0, 0, 0, 0, 2387 : 1, 0, 2, 1, 2, 0, 0, 0, 0, 2388 : 1, 0, 2, -1, 2, 0, 0, 0, 0, 2389 : -1, 0, 0, 2, -2, 0, 0, 0, 0, 2390 : 0, 0, 2, 1, 1, 0, 0, 0, 0 / 2391 DATA ( ( NA(I,J), I=1,9 ), J=141,150 ) / 2392 : -1, 0, 2, 0, -1, 0, 0, 0, 0, 2393 : -1, 0, 2, 4, 1, 0, 0, 0, 0, 2394 : 0, 0, 2, 2, 0, 0, 0, 0, 0, 2395 : 1, 1, 2, -2, 1, 0, 0, 0, 0, 2396 : 0, 0, 1, 0, 1, 0, 0, 0, 0, 2397 : -1, 0, 2, -1, 1, 0, 0, 0, 0, 2398 : -2, 0, 2, 2, 1, 0, 0, 0, 0, 2399 : 2, -1, 0, 0, 0, 0, 0, 0, 0, 2400 : 4, 0, 2, 0, 2, 0, 0, 0, 0, 2401 : 2, 1, 2, -2, 2, 0, 0, 0, 0 / 2402 DATA ( ( NA(I,J), I=1,9 ), J=151,160 ) / 2403 : 0, 1, 2, 1, 2, 0, 0, 0, 0, 2404 : 1, 0, 4, -2, 2, 0, 0, 0, 0, 2405 : 1, 1, 0, 0, -1, 0, 0, 0, 0, 2406 : -2, 0, 2, 4, 1, 0, 0, 0, 0, 2407 : 2, 0, 2, 0, 0, 0, 0, 0, 0, 2408 : -1, 0, 1, 0, 0, 0, 0, 0, 0, 2409 : 1, 0, 0, 1, 0, 0, 0, 0, 0, 2410 : 0, 1, 0, 2, 1, 0, 0, 0, 0, 2411 : -1, 0, 4, 0, 1, 0, 0, 0, 0, 2412 : -1, 0, 0, 4, 1, 0, 0, 0, 0 / 2413 DATA ( ( NA(I,J), I=1,9 ), J=161,170 ) / 2414 : 2, 0, 2, 2, 1, 0, 0, 0, 0, 2415 : 2, 1, 0, 0, 0, 0, 0, 0, 0, 2416 : 0, 0, 5, -5, 5, -3, 0, 0, 0, 2417 : 0, 0, 0, 0, 0, 0, 0, 2, 0, 2418 : 0, 0, 1, -1, 1, 0, 0, -1, 0, 2419 : 0, 0, -1, 1, -1, 1, 0, 0, 0, 2420 : 0, 0, -1, 1, 0, 0, 2, 0, 0, 2421 : 0, 0, 3, -3, 3, 0, 0, -1, 0, 2422 : 0, 0, -8, 8, -7, 5, 0, 0, 0, 2423 : 0, 0, -1, 1, -1, 0, 2, 0, 0 / 2424 DATA ( ( NA(I,J), I=1,9 ), J=171,180 ) / 2425 : 0, 0, -2, 2, -2, 2, 0, 0, 0, 2426 : 0, 0, -6, 6, -6, 4, 0, 0, 0, 2427 : 0, 0, -2, 2, -2, 0, 8, -3, 0, 2428 : 0, 0, 6, -6, 6, 0, -8, 3, 0, 2429 : 0, 0, 4, -4, 4, -2, 0, 0, 0, 2430 : 0, 0, -3, 3, -3, 2, 0, 0, 0, 2431 : 0, 0, 4, -4, 3, 0, -8, 3, 0, 2432 : 0, 0, -4, 4, -5, 0, 8, -3, 0, 2433 : 0, 0, 0, 0, 0, 2, 0, 0, 0, 2434 : 0, 0, -4, 4, -4, 3, 0, 0, 0 / 2435 DATA ( ( NA(I,J), I=1,9 ), J=181,190 ) / 2436 : 0, 1, -1, 1, -1, 0, 0, 1, 0, 2437 : 0, 0, 0, 0, 0, 0, 0, 1, 0, 2438 : 0, 0, 1, -1, 1, 1, 0, 0, 0, 2439 : 0, 0, 2, -2, 2, 0, -2, 0, 0, 2440 : 0, -1, -7, 7, -7, 5, 0, 0, 0, 2441 : -2, 0, 2, 0, 2, 0, 0, -2, 0, 2442 : -2, 0, 2, 0, 1, 0, 0, -3, 0, 2443 : 0, 0, 2, -2, 2, 0, 0, -2, 0, 2444 : 0, 0, 1, -1, 1, 0, 0, 1, 0, 2445 : 0, 0, 0, 0, 0, 0, 0, 0, 2 / 2446 DATA ( ( NA(I,J), I=1,9 ), J=191,NTERMS ) / 2447 : 0, 0, 0, 0, 0, 0, 0, 0, 1, 2448 : 2, 0, -2, 0, -2, 0, 0, 3, 0, 2449 : 0, 0, 1, -1, 1, 0, 0, -2, 0, 2450 : 0, 0, -7, 7, -7, 5, 0, 0, 0 / 2451 2452* Nutation series: longitude 2453 DATA ( ( PSI(I,J), I=1,4 ), J=1,10 ) / 2454 : 3341.5D0, 17206241.8D0, 3.1D0, 17409.5D0, 2455 : -1716.8D0, -1317185.3D0, 1.4D0, -156.8D0, 2456 : 285.7D0, -227667.0D0, 0.3D0, -23.5D0, 2457 : -68.6D0, -207448.0D0, 0.0D0, -21.4D0, 2458 : 950.3D0, 147607.9D0, -2.3D0, -355.0D0, 2459 : -66.7D0, -51689.1D0, 0.2D0, 122.6D0, 2460 : -108.6D0, 71117.6D0, 0.0D0, 7.0D0, 2461 : 35.6D0, -38740.2D0, 0.1D0, -36.2D0, 2462 : 85.4D0, -30127.6D0, 0.0D0, -3.1D0, 2463 : 9.0D0, 21583.0D0, 0.1D0, -50.3D0 / 2464 DATA ( ( PSI(I,J), I=1,4 ), J=11,20 ) / 2465 : 22.1D0, 12822.8D0, 0.0D0, 13.3D0, 2466 : 3.4D0, 12350.8D0, 0.0D0, 1.3D0, 2467 : -21.1D0, 15699.4D0, 0.0D0, 1.6D0, 2468 : 4.2D0, 6313.8D0, 0.0D0, 6.2D0, 2469 : -22.8D0, 5796.9D0, 0.0D0, 6.1D0, 2470 : 15.7D0, -5961.1D0, 0.0D0, -0.6D0, 2471 : 13.1D0, -5159.1D0, 0.0D0, -4.6D0, 2472 : 1.8D0, 4592.7D0, 0.0D0, 4.5D0, 2473 : -17.5D0, 6336.0D0, 0.0D0, 0.7D0, 2474 : 16.3D0, -3851.1D0, 0.0D0, -0.4D0 / 2475 DATA ( ( PSI(I,J), I=1,4 ), J=21,30 ) / 2476 : -2.8D0, 4771.7D0, 0.0D0, 0.5D0, 2477 : 13.8D0, -3099.3D0, 0.0D0, -0.3D0, 2478 : 0.2D0, 2860.3D0, 0.0D0, 0.3D0, 2479 : 1.4D0, 2045.3D0, 0.0D0, 2.0D0, 2480 : -8.6D0, 2922.6D0, 0.0D0, 0.3D0, 2481 : -7.7D0, 2587.9D0, 0.0D0, 0.2D0, 2482 : 8.8D0, -1408.1D0, 0.0D0, 3.7D0, 2483 : 1.4D0, 1517.5D0, 0.0D0, 1.5D0, 2484 : -1.9D0, -1579.7D0, 0.0D0, 7.7D0, 2485 : 1.3D0, -2178.6D0, 0.0D0, -0.2D0 / 2486 DATA ( ( PSI(I,J), I=1,4 ), J=31,40 ) / 2487 : -4.8D0, 1286.8D0, 0.0D0, 1.3D0, 2488 : 6.3D0, 1267.2D0, 0.0D0, -4.0D0, 2489 : -1.0D0, 1669.3D0, 0.0D0, -8.3D0, 2490 : 2.4D0, -1020.0D0, 0.0D0, -0.9D0, 2491 : 4.5D0, -766.9D0, 0.0D0, 0.0D0, 2492 : -1.1D0, 756.5D0, 0.0D0, -1.7D0, 2493 : -1.4D0, -1097.3D0, 0.0D0, -0.5D0, 2494 : 2.6D0, -663.0D0, 0.0D0, -0.6D0, 2495 : 0.8D0, -714.1D0, 0.0D0, 1.6D0, 2496 : 0.4D0, -629.9D0, 0.0D0, -0.6D0 / 2497 DATA ( ( PSI(I,J), I=1,4 ), J=41,50 ) / 2498 : 0.3D0, 580.4D0, 0.0D0, 0.6D0, 2499 : -1.6D0, 577.3D0, 0.0D0, 0.5D0, 2500 : -0.9D0, 644.4D0, 0.0D0, 0.0D0, 2501 : 2.2D0, -534.0D0, 0.0D0, -0.5D0, 2502 : -2.5D0, 493.3D0, 0.0D0, 0.5D0, 2503 : -0.1D0, -477.3D0, 0.0D0, -2.4D0, 2504 : -0.9D0, 735.0D0, 0.0D0, -1.7D0, 2505 : 0.7D0, 406.2D0, 0.0D0, 0.4D0, 2506 : -2.8D0, 656.9D0, 0.0D0, 0.0D0, 2507 : 0.6D0, 358.0D0, 0.0D0, 2.0D0 / 2508 DATA ( ( PSI(I,J), I=1,4 ), J=51,60 ) / 2509 : -0.7D0, 472.5D0, 0.0D0, -1.1D0, 2510 : -0.1D0, -300.5D0, 0.0D0, 0.0D0, 2511 : -1.2D0, 435.1D0, 0.0D0, -1.0D0, 2512 : 1.8D0, -289.4D0, 0.0D0, 0.0D0, 2513 : 0.6D0, -422.6D0, 0.0D0, 0.0D0, 2514 : 0.8D0, -287.6D0, 0.0D0, 0.6D0, 2515 : -38.6D0, -392.3D0, 0.0D0, 0.0D0, 2516 : 0.7D0, -281.8D0, 0.0D0, 0.6D0, 2517 : 0.6D0, -405.7D0, 0.0D0, 0.0D0, 2518 : -1.2D0, 229.0D0, 0.0D0, 0.2D0 / 2519 DATA ( ( PSI(I,J), I=1,4 ), J=61,70 ) / 2520 : 1.1D0, -264.3D0, 0.0D0, 0.5D0, 2521 : -0.7D0, 247.9D0, 0.0D0, -0.5D0, 2522 : -0.2D0, 218.0D0, 0.0D0, 0.2D0, 2523 : 0.6D0, -339.0D0, 0.0D0, 0.8D0, 2524 : -0.7D0, 198.7D0, 0.0D0, 0.2D0, 2525 : -1.5D0, 334.0D0, 0.0D0, 0.0D0, 2526 : 0.1D0, 334.0D0, 0.0D0, 0.0D0, 2527 : -0.1D0, -198.1D0, 0.0D0, 0.0D0, 2528 : -106.6D0, 0.0D0, 0.0D0, 0.0D0, 2529 : -0.5D0, 165.8D0, 0.0D0, 0.0D0 / 2530 DATA ( ( PSI(I,J), I=1,4 ), J=71,80 ) / 2531 : 0.0D0, 134.8D0, 0.0D0, 0.0D0, 2532 : 0.9D0, -151.6D0, 0.0D0, 0.0D0, 2533 : 0.0D0, -129.7D0, 0.0D0, 0.0D0, 2534 : 0.8D0, -132.8D0, 0.0D0, -0.1D0, 2535 : 0.5D0, -140.7D0, 0.0D0, 0.0D0, 2536 : -0.1D0, 138.4D0, 0.0D0, 0.0D0, 2537 : 0.0D0, 129.0D0, 0.0D0, -0.3D0, 2538 : 0.5D0, -121.2D0, 0.0D0, 0.0D0, 2539 : -0.3D0, 114.5D0, 0.0D0, 0.0D0, 2540 : -0.1D0, 101.8D0, 0.0D0, 0.0D0 / 2541 DATA ( ( PSI(I,J), I=1,4 ), J=81,90 ) / 2542 : -3.6D0, -101.9D0, 0.0D0, 0.0D0, 2543 : 0.8D0, -109.4D0, 0.0D0, 0.0D0, 2544 : 0.2D0, -97.0D0, 0.0D0, 0.0D0, 2545 : -0.7D0, 157.3D0, 0.0D0, 0.0D0, 2546 : 0.2D0, -83.3D0, 0.0D0, 0.0D0, 2547 : -0.3D0, 93.3D0, 0.0D0, 0.0D0, 2548 : -0.1D0, 92.1D0, 0.0D0, 0.0D0, 2549 : -0.5D0, 133.6D0, 0.0D0, 0.0D0, 2550 : -0.1D0, 81.5D0, 0.0D0, 0.0D0, 2551 : 0.0D0, 123.9D0, 0.0D0, 0.0D0 / 2552 DATA ( ( PSI(I,J), I=1,4 ), J=91,100 ) / 2553 : -0.3D0, 128.1D0, 0.0D0, 0.0D0, 2554 : 0.1D0, 74.1D0, 0.0D0, -0.3D0, 2555 : -0.2D0, -70.3D0, 0.0D0, 0.0D0, 2556 : -0.4D0, 66.6D0, 0.0D0, 0.0D0, 2557 : 0.1D0, -66.7D0, 0.0D0, 0.0D0, 2558 : -0.7D0, 69.3D0, 0.0D0, -0.3D0, 2559 : 0.0D0, -70.4D0, 0.0D0, 0.0D0, 2560 : -0.1D0, 101.5D0, 0.0D0, 0.0D0, 2561 : 0.5D0, -69.1D0, 0.0D0, 0.0D0, 2562 : -0.2D0, 58.5D0, 0.0D0, 0.2D0 / 2563 DATA ( ( PSI(I,J), I=1,4 ), J=101,110 ) / 2564 : 0.1D0, -94.9D0, 0.0D0, 0.2D0, 2565 : 0.0D0, 52.9D0, 0.0D0, -0.2D0, 2566 : 0.1D0, 86.7D0, 0.0D0, -0.2D0, 2567 : -0.1D0, -59.2D0, 0.0D0, 0.2D0, 2568 : 0.3D0, -58.8D0, 0.0D0, 0.1D0, 2569 : -0.3D0, 49.0D0, 0.0D0, 0.0D0, 2570 : -0.2D0, 56.9D0, 0.0D0, -0.1D0, 2571 : 0.3D0, -50.2D0, 0.0D0, 0.0D0, 2572 : -0.2D0, 53.4D0, 0.0D0, -0.1D0, 2573 : 0.1D0, -76.5D0, 0.0D0, 0.0D0 / 2574 DATA ( ( PSI(I,J), I=1,4 ), J=111,120 ) / 2575 : -0.2D0, 45.3D0, 0.0D0, 0.0D0, 2576 : 0.1D0, -46.8D0, 0.0D0, 0.0D0, 2577 : 0.2D0, -44.6D0, 0.0D0, 0.0D0, 2578 : 0.2D0, -48.7D0, 0.0D0, 0.0D0, 2579 : 0.1D0, -46.8D0, 0.0D0, 0.0D0, 2580 : 0.1D0, -42.0D0, 0.0D0, 0.0D0, 2581 : 0.0D0, 46.4D0, 0.0D0, -0.1D0, 2582 : 0.2D0, -67.3D0, 0.0D0, 0.1D0, 2583 : 0.0D0, -65.8D0, 0.0D0, 0.2D0, 2584 : -0.1D0, -43.9D0, 0.0D0, 0.3D0 / 2585 DATA ( ( PSI(I,J), I=1,4 ), J=121,130 ) / 2586 : 0.0D0, -38.9D0, 0.0D0, 0.0D0, 2587 : -0.3D0, 63.9D0, 0.0D0, 0.0D0, 2588 : -0.2D0, 41.2D0, 0.0D0, 0.0D0, 2589 : 0.0D0, -36.1D0, 0.0D0, 0.2D0, 2590 : -0.3D0, 58.5D0, 0.0D0, 0.0D0, 2591 : -0.1D0, 36.1D0, 0.0D0, 0.0D0, 2592 : 0.0D0, -39.7D0, 0.0D0, 0.0D0, 2593 : 0.1D0, -57.7D0, 0.0D0, 0.0D0, 2594 : -0.2D0, 33.4D0, 0.0D0, 0.0D0, 2595 : 36.4D0, 0.0D0, 0.0D0, 0.0D0 / 2596 DATA ( ( PSI(I,J), I=1,4 ), J=131,140 ) / 2597 : -0.1D0, 55.7D0, 0.0D0, -0.1D0, 2598 : 0.1D0, -35.4D0, 0.0D0, 0.0D0, 2599 : 0.1D0, -31.0D0, 0.0D0, 0.0D0, 2600 : -0.1D0, 30.1D0, 0.0D0, 0.0D0, 2601 : -0.3D0, 49.2D0, 0.0D0, 0.0D0, 2602 : -0.2D0, 49.1D0, 0.0D0, 0.0D0, 2603 : -0.1D0, 33.6D0, 0.0D0, 0.0D0, 2604 : 0.1D0, -33.5D0, 0.0D0, 0.0D0, 2605 : 0.1D0, -31.0D0, 0.0D0, 0.0D0, 2606 : -0.1D0, 28.0D0, 0.0D0, 0.0D0 / 2607 DATA ( ( PSI(I,J), I=1,4 ), J=141,150 ) / 2608 : 0.1D0, -25.2D0, 0.0D0, 0.0D0, 2609 : 0.1D0, -26.2D0, 0.0D0, 0.0D0, 2610 : -0.2D0, 41.5D0, 0.0D0, 0.0D0, 2611 : 0.0D0, 24.5D0, 0.0D0, 0.1D0, 2612 : -16.2D0, 0.0D0, 0.0D0, 0.0D0, 2613 : 0.0D0, -22.3D0, 0.0D0, 0.0D0, 2614 : 0.0D0, 23.1D0, 0.0D0, 0.0D0, 2615 : -0.1D0, 37.5D0, 0.0D0, 0.0D0, 2616 : 0.2D0, -25.7D0, 0.0D0, 0.0D0, 2617 : 0.0D0, 25.2D0, 0.0D0, 0.0D0 / 2618 DATA ( ( PSI(I,J), I=1,4 ), J=151,160 ) / 2619 : 0.1D0, -24.5D0, 0.0D0, 0.0D0, 2620 : -0.1D0, 24.3D0, 0.0D0, 0.0D0, 2621 : 0.1D0, -20.7D0, 0.0D0, 0.0D0, 2622 : 0.1D0, -20.8D0, 0.0D0, 0.0D0, 2623 : -0.2D0, 33.4D0, 0.0D0, 0.0D0, 2624 : 32.9D0, 0.0D0, 0.0D0, 0.0D0, 2625 : 0.1D0, -32.6D0, 0.0D0, 0.0D0, 2626 : 0.0D0, 19.9D0, 0.0D0, 0.0D0, 2627 : -0.1D0, 19.6D0, 0.0D0, 0.0D0, 2628 : 0.0D0, -18.7D0, 0.0D0, 0.0D0 / 2629 DATA ( ( PSI(I,J), I=1,4 ), J=161,170 ) / 2630 : 0.1D0, -19.0D0, 0.0D0, 0.0D0, 2631 : 0.1D0, -28.6D0, 0.0D0, 0.0D0, 2632 : 4.0D0, 178.8D0,-11.8D0, 0.3D0, 2633 : 39.8D0, -107.3D0, -5.6D0, -1.0D0, 2634 : 9.9D0, 164.0D0, -4.1D0, 0.1D0, 2635 : -4.8D0, -135.3D0, -3.4D0, -0.1D0, 2636 : 50.5D0, 75.0D0, 1.4D0, -1.2D0, 2637 : -1.1D0, -53.5D0, 1.3D0, 0.0D0, 2638 : -45.0D0, -2.4D0, -0.4D0, 6.6D0, 2639 : -11.5D0, -61.0D0, -0.9D0, 0.4D0 / 2640 DATA ( ( PSI(I,J), I=1,4 ), J=171,180 ) / 2641 : 4.4D0, -68.4D0, -3.4D0, 0.0D0, 2642 : 7.7D0, -47.1D0, -4.7D0, -1.0D0, 2643 : -42.9D0, -12.6D0, -1.2D0, 4.2D0, 2644 : -42.8D0, 12.7D0, -1.2D0, -4.2D0, 2645 : -7.6D0, -44.1D0, 2.1D0, -0.5D0, 2646 : -64.1D0, 1.7D0, 0.2D0, 4.5D0, 2647 : 36.4D0, -10.4D0, 1.0D0, 3.5D0, 2648 : 35.6D0, 10.2D0, 1.0D0, -3.5D0, 2649 : -1.7D0, 39.5D0, 2.0D0, 0.0D0, 2650 : 50.9D0, -8.2D0, -0.8D0, -5.0D0 / 2651 DATA ( ( PSI(I,J), I=1,4 ), J=181,190 ) / 2652 : 0.0D0, 52.3D0, 1.2D0, 0.0D0, 2653 : -42.9D0, -17.8D0, 0.4D0, 0.0D0, 2654 : 2.6D0, 34.3D0, 0.8D0, 0.0D0, 2655 : -0.8D0, -48.6D0, 2.4D0, -0.1D0, 2656 : -4.9D0, 30.5D0, 3.7D0, 0.7D0, 2657 : 0.0D0, -43.6D0, 2.1D0, 0.0D0, 2658 : 0.0D0, -25.4D0, 1.2D0, 0.0D0, 2659 : 2.0D0, 40.9D0, -2.0D0, 0.0D0, 2660 : -2.1D0, 26.1D0, 0.6D0, 0.0D0, 2661 : 22.6D0, -3.2D0, -0.5D0, -0.5D0 / 2662 DATA ( ( PSI(I,J), I=1,4 ), J=191,NTERMS ) / 2663 : -7.6D0, 24.9D0, -0.4D0, -0.2D0, 2664 : -6.2D0, 34.9D0, 1.7D0, 0.3D0, 2665 : 2.0D0, 17.4D0, -0.4D0, 0.1D0, 2666 : -3.9D0, 20.5D0, 2.4D0, 0.6D0 / 2667 2668* Nutation series: obliquity 2669 DATA ( ( EPS(I,J), I=1,4 ), J=1,10 ) / 2670 : 9205365.8D0, -1506.2D0, 885.7D0, -0.2D0, 2671 : 573095.9D0, -570.2D0, -305.0D0, -0.3D0, 2672 : 97845.5D0, 147.8D0, -48.8D0, -0.2D0, 2673 : -89753.6D0, 28.0D0, 46.9D0, 0.0D0, 2674 : 7406.7D0, -327.1D0, -18.2D0, 0.8D0, 2675 : 22442.3D0, -22.3D0, -67.6D0, 0.0D0, 2676 : -683.6D0, 46.8D0, 0.0D0, 0.0D0, 2677 : 20070.7D0, 36.0D0, 1.6D0, 0.0D0, 2678 : 12893.8D0, 39.5D0, -6.2D0, 0.0D0, 2679 : -9593.2D0, 14.4D0, 30.2D0, -0.1D0 / 2680 DATA ( ( EPS(I,J), I=1,4 ), J=11,20 ) / 2681 : -6899.5D0, 4.8D0, -0.6D0, 0.0D0, 2682 : -5332.5D0, -0.1D0, 2.7D0, 0.0D0, 2683 : -125.2D0, 10.5D0, 0.0D0, 0.0D0, 2684 : -3323.4D0, -0.9D0, -0.3D0, 0.0D0, 2685 : 3142.3D0, 8.9D0, 0.3D0, 0.0D0, 2686 : 2552.5D0, 7.3D0, -1.2D0, 0.0D0, 2687 : 2634.4D0, 8.8D0, 0.2D0, 0.0D0, 2688 : -2424.4D0, 1.6D0, -0.4D0, 0.0D0, 2689 : -123.3D0, 3.9D0, 0.0D0, 0.0D0, 2690 : 1642.4D0, 7.3D0, -0.8D0, 0.0D0 / 2691 DATA ( ( EPS(I,J), I=1,4 ), J=21,30 ) / 2692 : 47.9D0, 3.2D0, 0.0D0, 0.0D0, 2693 : 1321.2D0, 6.2D0, -0.6D0, 0.0D0, 2694 : -1234.1D0, -0.3D0, 0.6D0, 0.0D0, 2695 : -1076.5D0, -0.3D0, 0.0D0, 0.0D0, 2696 : -61.6D0, 1.8D0, 0.0D0, 0.0D0, 2697 : -55.4D0, 1.6D0, 0.0D0, 0.0D0, 2698 : 856.9D0, -4.9D0, -2.1D0, 0.0D0, 2699 : -800.7D0, -0.1D0, 0.0D0, 0.0D0, 2700 : 685.1D0, -0.6D0, -3.8D0, 0.0D0, 2701 : -16.9D0, -1.5D0, 0.0D0, 0.0D0 / 2702 DATA ( ( EPS(I,J), I=1,4 ), J=31,40 ) / 2703 : 695.7D0, 1.8D0, 0.0D0, 0.0D0, 2704 : 642.2D0, -2.6D0, -1.6D0, 0.0D0, 2705 : 13.3D0, 1.1D0, -0.1D0, 0.0D0, 2706 : 521.9D0, 1.6D0, 0.0D0, 0.0D0, 2707 : 325.8D0, 2.0D0, -0.1D0, 0.0D0, 2708 : -325.1D0, -0.5D0, 0.9D0, 0.0D0, 2709 : 10.1D0, 0.3D0, 0.0D0, 0.0D0, 2710 : 334.5D0, 1.6D0, 0.0D0, 0.0D0, 2711 : 307.1D0, 0.4D0, -0.9D0, 0.0D0, 2712 : 327.2D0, 0.5D0, 0.0D0, 0.0D0 / 2713 DATA ( ( EPS(I,J), I=1,4 ), J=41,50 ) / 2714 : -304.6D0, -0.1D0, 0.0D0, 0.0D0, 2715 : 304.0D0, 0.6D0, 0.0D0, 0.0D0, 2716 : -276.8D0, -0.5D0, 0.1D0, 0.0D0, 2717 : 268.9D0, 1.3D0, 0.0D0, 0.0D0, 2718 : 271.8D0, 1.1D0, 0.0D0, 0.0D0, 2719 : 271.5D0, -0.4D0, -0.8D0, 0.0D0, 2720 : -5.2D0, 0.5D0, 0.0D0, 0.0D0, 2721 : -220.5D0, 0.1D0, 0.0D0, 0.0D0, 2722 : -20.1D0, 0.3D0, 0.0D0, 0.0D0, 2723 : -191.0D0, 0.1D0, 0.5D0, 0.0D0 / 2724 DATA ( ( EPS(I,J), I=1,4 ), J=51,60 ) / 2725 : -4.1D0, 0.3D0, 0.0D0, 0.0D0, 2726 : 130.6D0, -0.1D0, 0.0D0, 0.0D0, 2727 : 3.0D0, 0.3D0, 0.0D0, 0.0D0, 2728 : 122.9D0, 0.8D0, 0.0D0, 0.0D0, 2729 : 3.7D0, -0.3D0, 0.0D0, 0.0D0, 2730 : 123.1D0, 0.4D0, -0.3D0, 0.0D0, 2731 : -52.7D0, 15.3D0, 0.0D0, 0.0D0, 2732 : 120.7D0, 0.3D0, -0.3D0, 0.0D0, 2733 : 4.0D0, -0.3D0, 0.0D0, 0.0D0, 2734 : 126.5D0, 0.5D0, 0.0D0, 0.0D0 / 2735 DATA ( ( EPS(I,J), I=1,4 ), J=61,70 ) / 2736 : 112.7D0, 0.5D0, -0.3D0, 0.0D0, 2737 : -106.1D0, -0.3D0, 0.3D0, 0.0D0, 2738 : -112.9D0, -0.2D0, 0.0D0, 0.0D0, 2739 : 3.6D0, -0.2D0, 0.0D0, 0.0D0, 2740 : 107.4D0, 0.3D0, 0.0D0, 0.0D0, 2741 : -10.9D0, 0.2D0, 0.0D0, 0.0D0, 2742 : -0.9D0, 0.0D0, 0.0D0, 0.0D0, 2743 : 85.4D0, 0.0D0, 0.0D0, 0.0D0, 2744 : 0.0D0, -88.8D0, 0.0D0, 0.0D0, 2745 : -71.0D0, -0.2D0, 0.0D0, 0.0D0 / 2746 DATA ( ( EPS(I,J), I=1,4 ), J=71,80 ) / 2747 : -70.3D0, 0.0D0, 0.0D0, 0.0D0, 2748 : 64.5D0, 0.4D0, 0.0D0, 0.0D0, 2749 : 69.8D0, 0.0D0, 0.0D0, 0.0D0, 2750 : 66.1D0, 0.4D0, 0.0D0, 0.0D0, 2751 : -61.0D0, -0.2D0, 0.0D0, 0.0D0, 2752 : -59.5D0, -0.1D0, 0.0D0, 0.0D0, 2753 : -55.6D0, 0.0D0, 0.2D0, 0.0D0, 2754 : 51.7D0, 0.2D0, 0.0D0, 0.0D0, 2755 : -49.0D0, -0.1D0, 0.0D0, 0.0D0, 2756 : -52.7D0, -0.1D0, 0.0D0, 0.0D0 / 2757 DATA ( ( EPS(I,J), I=1,4 ), J=81,90 ) / 2758 : -49.6D0, 1.4D0, 0.0D0, 0.0D0, 2759 : 46.3D0, 0.4D0, 0.0D0, 0.0D0, 2760 : 49.6D0, 0.1D0, 0.0D0, 0.0D0, 2761 : -5.1D0, 0.1D0, 0.0D0, 0.0D0, 2762 : -44.0D0, -0.1D0, 0.0D0, 0.0D0, 2763 : -39.9D0, -0.1D0, 0.0D0, 0.0D0, 2764 : -39.5D0, -0.1D0, 0.0D0, 0.0D0, 2765 : -3.9D0, 0.1D0, 0.0D0, 0.0D0, 2766 : -42.1D0, -0.1D0, 0.0D0, 0.0D0, 2767 : -17.2D0, 0.1D0, 0.0D0, 0.0D0 / 2768 DATA ( ( EPS(I,J), I=1,4 ), J=91,100 ) / 2769 : -2.3D0, 0.1D0, 0.0D0, 0.0D0, 2770 : -39.2D0, 0.0D0, 0.0D0, 0.0D0, 2771 : -38.4D0, 0.1D0, 0.0D0, 0.0D0, 2772 : 36.8D0, 0.2D0, 0.0D0, 0.0D0, 2773 : 34.6D0, 0.1D0, 0.0D0, 0.0D0, 2774 : -32.7D0, 0.3D0, 0.0D0, 0.0D0, 2775 : 30.4D0, 0.0D0, 0.0D0, 0.0D0, 2776 : 0.4D0, 0.1D0, 0.0D0, 0.0D0, 2777 : 29.3D0, 0.2D0, 0.0D0, 0.0D0, 2778 : 31.6D0, 0.1D0, 0.0D0, 0.0D0 / 2779 DATA ( ( EPS(I,J), I=1,4 ), J=101,110 ) / 2780 : 0.8D0, -0.1D0, 0.0D0, 0.0D0, 2781 : -27.9D0, 0.0D0, 0.0D0, 0.0D0, 2782 : 2.9D0, 0.0D0, 0.0D0, 0.0D0, 2783 : -25.3D0, 0.0D0, 0.0D0, 0.0D0, 2784 : 25.0D0, 0.1D0, 0.0D0, 0.0D0, 2785 : 27.5D0, 0.1D0, 0.0D0, 0.0D0, 2786 : -24.4D0, -0.1D0, 0.0D0, 0.0D0, 2787 : 24.9D0, 0.2D0, 0.0D0, 0.0D0, 2788 : -22.8D0, -0.1D0, 0.0D0, 0.0D0, 2789 : 0.9D0, -0.1D0, 0.0D0, 0.0D0 / 2790 DATA ( ( EPS(I,J), I=1,4 ), J=111,120 ) / 2791 : 24.4D0, 0.1D0, 0.0D0, 0.0D0, 2792 : 23.9D0, 0.1D0, 0.0D0, 0.0D0, 2793 : 22.5D0, 0.1D0, 0.0D0, 0.0D0, 2794 : 20.8D0, 0.1D0, 0.0D0, 0.0D0, 2795 : 20.1D0, 0.0D0, 0.0D0, 0.0D0, 2796 : 21.5D0, 0.1D0, 0.0D0, 0.0D0, 2797 : -20.0D0, 0.0D0, 0.0D0, 0.0D0, 2798 : 1.4D0, 0.0D0, 0.0D0, 0.0D0, 2799 : -0.2D0, -0.1D0, 0.0D0, 0.0D0, 2800 : 19.0D0, 0.0D0, -0.1D0, 0.0D0 / 2801 DATA ( ( EPS(I,J), I=1,4 ), J=121,130 ) / 2802 : 20.5D0, 0.0D0, 0.0D0, 0.0D0, 2803 : -2.0D0, 0.0D0, 0.0D0, 0.0D0, 2804 : -17.6D0, -0.1D0, 0.0D0, 0.0D0, 2805 : 19.0D0, 0.0D0, 0.0D0, 0.0D0, 2806 : -2.4D0, 0.0D0, 0.0D0, 0.0D0, 2807 : -18.4D0, -0.1D0, 0.0D0, 0.0D0, 2808 : 17.1D0, 0.0D0, 0.0D0, 0.0D0, 2809 : 0.4D0, 0.0D0, 0.0D0, 0.0D0, 2810 : 18.4D0, 0.1D0, 0.0D0, 0.0D0, 2811 : 0.0D0, 17.4D0, 0.0D0, 0.0D0 / 2812 DATA ( ( EPS(I,J), I=1,4 ), J=131,140 ) / 2813 : -0.6D0, 0.0D0, 0.0D0, 0.0D0, 2814 : -15.4D0, 0.0D0, 0.0D0, 0.0D0, 2815 : -16.8D0, -0.1D0, 0.0D0, 0.0D0, 2816 : 16.3D0, 0.0D0, 0.0D0, 0.0D0, 2817 : -2.0D0, 0.0D0, 0.0D0, 0.0D0, 2818 : -1.5D0, 0.0D0, 0.0D0, 0.0D0, 2819 : -14.3D0, -0.1D0, 0.0D0, 0.0D0, 2820 : 14.4D0, 0.0D0, 0.0D0, 0.0D0, 2821 : -13.4D0, 0.0D0, 0.0D0, 0.0D0, 2822 : -14.3D0, -0.1D0, 0.0D0, 0.0D0 / 2823 DATA ( ( EPS(I,J), I=1,4 ), J=141,150 ) / 2824 : -13.7D0, 0.0D0, 0.0D0, 0.0D0, 2825 : 13.1D0, 0.1D0, 0.0D0, 0.0D0, 2826 : -1.7D0, 0.0D0, 0.0D0, 0.0D0, 2827 : -12.8D0, 0.0D0, 0.0D0, 0.0D0, 2828 : 0.0D0, -14.4D0, 0.0D0, 0.0D0, 2829 : 12.4D0, 0.0D0, 0.0D0, 0.0D0, 2830 : -12.0D0, 0.0D0, 0.0D0, 0.0D0, 2831 : -0.8D0, 0.0D0, 0.0D0, 0.0D0, 2832 : 10.9D0, 0.1D0, 0.0D0, 0.0D0, 2833 : -10.8D0, 0.0D0, 0.0D0, 0.0D0 / 2834 DATA ( ( EPS(I,J), I=1,4 ), J=151,160 ) / 2835 : 10.5D0, 0.0D0, 0.0D0, 0.0D0, 2836 : -10.4D0, 0.0D0, 0.0D0, 0.0D0, 2837 : -11.2D0, 0.0D0, 0.0D0, 0.0D0, 2838 : 10.5D0, 0.1D0, 0.0D0, 0.0D0, 2839 : -1.4D0, 0.0D0, 0.0D0, 0.0D0, 2840 : 0.0D0, 0.1D0, 0.0D0, 0.0D0, 2841 : 0.7D0, 0.0D0, 0.0D0, 0.0D0, 2842 : -10.3D0, 0.0D0, 0.0D0, 0.0D0, 2843 : -10.0D0, 0.0D0, 0.0D0, 0.0D0, 2844 : 9.6D0, 0.0D0, 0.0D0, 0.0D0 / 2845 DATA ( ( EPS(I,J), I=1,4 ), J=161,170 ) / 2846 : 9.4D0, 0.1D0, 0.0D0, 0.0D0, 2847 : 0.6D0, 0.0D0, 0.0D0, 0.0D0, 2848 : -87.7D0, 4.4D0, -0.4D0, -6.3D0, 2849 : 46.3D0, 22.4D0, 0.5D0, -2.4D0, 2850 : 15.6D0, -3.4D0, 0.1D0, 0.4D0, 2851 : 5.2D0, 5.8D0, 0.2D0, -0.1D0, 2852 : -30.1D0, 26.9D0, 0.7D0, 0.0D0, 2853 : 23.2D0, -0.5D0, 0.0D0, 0.6D0, 2854 : 1.0D0, 23.2D0, 3.4D0, 0.0D0, 2855 : -12.2D0, -4.3D0, 0.0D0, 0.0D0 / 2856 DATA ( ( EPS(I,J), I=1,4 ), J=171,180 ) / 2857 : -2.1D0, -3.7D0, -0.2D0, 0.1D0, 2858 : -18.6D0, -3.8D0, -0.4D0, 1.8D0, 2859 : 5.5D0, -18.7D0, -1.8D0, -0.5D0, 2860 : -5.5D0, -18.7D0, 1.8D0, -0.5D0, 2861 : 18.4D0, -3.6D0, 0.3D0, 0.9D0, 2862 : -0.6D0, 1.3D0, 0.0D0, 0.0D0, 2863 : -5.6D0, -19.5D0, 1.9D0, 0.0D0, 2864 : 5.5D0, -19.1D0, -1.9D0, 0.0D0, 2865 : -17.3D0, -0.8D0, 0.0D0, 0.9D0, 2866 : -3.2D0, -8.3D0, -0.8D0, 0.3D0 / 2867 DATA ( ( EPS(I,J), I=1,4 ), J=181,190 ) / 2868 : -0.1D0, 0.0D0, 0.0D0, 0.0D0, 2869 : -5.4D0, 7.8D0, -0.3D0, 0.0D0, 2870 : -14.8D0, 1.4D0, 0.0D0, 0.3D0, 2871 : -3.8D0, 0.4D0, 0.0D0, -0.2D0, 2872 : 12.6D0, 3.2D0, 0.5D0, -1.5D0, 2873 : 0.1D0, 0.0D0, 0.0D0, 0.0D0, 2874 : -13.6D0, 2.4D0, -0.1D0, 0.0D0, 2875 : 0.9D0, 1.2D0, 0.0D0, 0.0D0, 2876 : -11.9D0, -0.5D0, 0.0D0, 0.3D0, 2877 : 0.4D0, 12.0D0, 0.3D0, -0.2D0 / 2878 DATA ( ( EPS(I,J), I=1,4 ), J=191,NTERMS ) / 2879 : 8.3D0, 6.1D0, -0.1D0, 0.1D0, 2880 : 0.0D0, 0.0D0, 0.0D0, 0.0D0, 2881 : 0.4D0, -10.8D0, 0.3D0, 0.0D0, 2882 : 9.6D0, 2.2D0, 0.3D0, -1.2D0 / 2883 2884 2885 2886* Interval between fundamental epoch J2000.0 and given epoch (JC). 2887 T = (DATE-DJM0)/DJC 2888 2889* Mean anomaly of the Moon. 2890 EL = 134.96340251D0*DD2R+ 2891 : MOD(T*(1717915923.2178D0+ 2892 : T*( 31.8792D0+ 2893 : T*( 0.051635D0+ 2894 : T*( - 0.00024470D0)))),TURNAS)*DAS2R 2895 2896* Mean anomaly of the Sun. 2897 ELP = 357.52910918D0*DD2R+ 2898 : MOD(T*( 129596581.0481D0+ 2899 : T*( - 0.5532D0+ 2900 : T*( 0.000136D0+ 2901 : T*( - 0.00001149D0)))),TURNAS)*DAS2R 2902 2903* Mean argument of the latitude of the Moon. 2904 F = 93.27209062D0*DD2R+ 2905 : MOD(T*(1739527262.8478D0+ 2906 : T*( - 12.7512D0+ 2907 : T*( - 0.001037D0+ 2908 : T*( 0.00000417D0)))),TURNAS)*DAS2R 2909 2910* Mean elongation of the Moon from the Sun. 2911 D = 297.85019547D0*DD2R+ 2912 : MOD(T*(1602961601.2090D0+ 2913 : T*( - 6.3706D0+ 2914 : T*( 0.006539D0+ 2915 : T*( - 0.00003169D0)))),TURNAS)*DAS2R 2916 2917* Mean longitude of the ascending node of the Moon. 2918 OM = 125.04455501D0*DD2R+ 2919 : MOD(T*( - 6962890.5431D0+ 2920 : T*( 7.4722D0+ 2921 : T*( 0.007702D0+ 2922 : T*( - 0.00005939D0)))),TURNAS)*DAS2R 2923 2924* Mean longitude of Venus. 2925 VE = 181.97980085D0*DD2R+MOD(210664136.433548D0*T,TURNAS)*DAS2R 2926 2927* Mean longitude of Mars. 2928 MA = 355.43299958D0*DD2R+MOD( 68905077.493988D0*T,TURNAS)*DAS2R 2929 2930* Mean longitude of Jupiter. 2931 JU = 34.35151874D0*DD2R+MOD( 10925660.377991D0*T,TURNAS)*DAS2R 2932 2933* Mean longitude of Saturn. 2934 SA = 50.07744430D0*DD2R+MOD( 4399609.855732D0*T,TURNAS)*DAS2R 2935 2936* Geodesic nutation (Fukushima 1991) in microarcsec. 2937 DP = -153.1D0*SIN(ELP)-1.9D0*SIN(2D0*ELP) 2938 DE = 0D0 2939 2940* Shirai & Fukushima (2001) nutation series. 2941 DO J=NTERMS,1,-1 2942 THETA = DBLE(NA(1,J))*EL+ 2943 : DBLE(NA(2,J))*ELP+ 2944 : DBLE(NA(3,J))*F+ 2945 : DBLE(NA(4,J))*D+ 2946 : DBLE(NA(5,J))*OM+ 2947 : DBLE(NA(6,J))*VE+ 2948 : DBLE(NA(7,J))*MA+ 2949 : DBLE(NA(8,J))*JU+ 2950 : DBLE(NA(9,J))*SA 2951 C = COS(THETA) 2952 S = SIN(THETA) 2953 DP = DP+(PSI(1,J)+PSI(3,J)*T)*C+(PSI(2,J)+PSI(4,J)*T)*S 2954 DE = DE+(EPS(1,J)+EPS(3,J)*T)*C+(EPS(2,J)+EPS(4,J)*T)*S 2955 END DO 2956 2957* Change of units, and addition of the precession correction. 2958 DPSI = (DP*1D-6-0.042888D0-0.29856D0*T)*DAS2R 2959 DEPS = (DE*1D-6-0.005171D0-0.02408D0*T)*DAS2R 2960 2961* Mean obliquity of date (Simon et al. 1994). 2962 EPS0 = (84381.412D0+ 2963 : (-46.80927D0+ 2964 : (-0.000152D0+ 2965 : (0.0019989D0+ 2966 : (-0.00000051D0+ 2967 : (-0.000000025D0)*T)*T)*T)*T)*T)*DAS2R 2968 2969 END 2970 SUBROUTINE sla_NUT (DATE, RMATN) 2971*+ 2972* - - - - 2973* N U T 2974* - - - - 2975* 2976* Form the matrix of nutation for a given date - Shirai & Fukushima 2977* 2001 theory (double precision) 2978* 2979* Reference: 2980* Shirai, T. & Fukushima, T., Astron.J. 121, 3270-3283 (2001). 2981* 2982* Given: 2983* DATE d TDB (loosely ET) as Modified Julian Date 2984* (=JD-2400000.5) 2985* Returned: 2986* RMATN d(3,3) nutation matrix 2987* 2988* Notes: 2989* 2990* 1 The matrix is in the sense v(true) = rmatn * v(mean) . 2991* where v(true) is the star vector relative to the true equator and 2992* equinox of date and v(mean) is the star vector relative to the 2993* mean equator and equinox of date. 2994* 2995* 2 The matrix represents forced nutation (but not free core 2996* nutation) plus corrections to the IAU~1976 precession model. 2997* 2998* 3 Earth attitude predictions made by combining the present nutation 2999* matrix with IAU~1976 precession are accurate to 1~mas (with 3000* respect to the ICRS) for a few decades around 2000. 3001* 3002* 4 The distinction between the required TDB and TT is always 3003* negligible. Moreover, for all but the most critical applications 3004* UTC is adequate. 3005* 3006* Called: sla_NUTC, sla_DEULER 3007* 3008* Last revision: 1 December 2005 3009* 3010* Copyright P.T.Wallace. All rights reserved. 3011* 3012* License: 3013* This program is free software; you can redistribute it and/or modify 3014* it under the terms of the GNU General Public License as published by 3015* the Free Software Foundation; either version 2 of the License, or 3016* (at your option) any later version. 3017* 3018* This program is distributed in the hope that it will be useful, 3019* but WITHOUT ANY WARRANTY; without even the implied warranty of 3020* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 3021* GNU General Public License for more details. 3022* 3023* You should have received a copy of the GNU General Public License 3024* along with this program (see SLA_CONDITIONS); if not, write to the 3025* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 3026* Boston, MA 02110-1301 USA 3027* 3028*- 3029 3030 IMPLICIT NONE 3031 3032 DOUBLE PRECISION DATE,RMATN(3,3) 3033 3034 DOUBLE PRECISION DPSI,DEPS,EPS0 3035 3036 3037 3038* Nutation components and mean obliquity 3039 CALL sla_NUTC(DATE,DPSI,DEPS,EPS0) 3040 3041* Rotation matrix 3042 CALL sla_DEULER('XZX',EPS0,-DPSI,-(EPS0+DEPS),RMATN) 3043 3044 END 3045 SUBROUTINE sla_PREBN (BEP0, BEP1, RMATP) 3046*+ 3047* - - - - - - 3048* P R E B N 3049* - - - - - - 3050* 3051* Generate the matrix of precession between two epochs, 3052* using the old, pre-IAU1976, Bessel-Newcomb model, using 3053* Kinoshita's formulation (double precision) 3054* 3055* Given: 3056* BEP0 dp beginning Besselian epoch 3057* BEP1 dp ending Besselian epoch 3058* 3059* Returned: 3060* RMATP dp(3,3) precession matrix 3061* 3062* The matrix is in the sense V(BEP1) = RMATP * V(BEP0) 3063* 3064* Reference: 3065* Kinoshita, H. (1975) 'Formulas for precession', SAO Special 3066* Report No. 364, Smithsonian Institution Astrophysical 3067* Observatory, Cambridge, Massachusetts. 3068* 3069* Called: sla_DEULER 3070* 3071* P.T.Wallace Starlink 23 August 1996 3072* 3073* Copyright (C) 1996 Rutherford Appleton Laboratory 3074* 3075* License: 3076* This program is free software; you can redistribute it and/or modify 3077* it under the terms of the GNU General Public License as published by 3078* the Free Software Foundation; either version 2 of the License, or 3079* (at your option) any later version. 3080* 3081* This program is distributed in the hope that it will be useful, 3082* but WITHOUT ANY WARRANTY; without even the implied warranty of 3083* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 3084* GNU General Public License for more details. 3085* 3086* You should have received a copy of the GNU General Public License 3087* along with this program (see SLA_CONDITIONS); if not, write to the 3088* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 3089* Boston, MA 02110-1301 USA 3090* 3091*- 3092 3093 IMPLICIT NONE 3094 3095 DOUBLE PRECISION BEP0,BEP1,RMATP(3,3) 3096 3097* Arc seconds to radians 3098 DOUBLE PRECISION AS2R 3099 PARAMETER (AS2R=0.484813681109535994D-5) 3100 3101 DOUBLE PRECISION BIGT,T,TAS2R,W,ZETA,Z,THETA 3102 3103 3104 3105* Interval between basic epoch B1850.0 and beginning epoch in TC 3106 BIGT = (BEP0-1850D0)/100D0 3107 3108* Interval over which precession required, in tropical centuries 3109 T = (BEP1-BEP0)/100D0 3110 3111* Euler angles 3112 TAS2R = T*AS2R 3113 W = 2303.5548D0+(1.39720D0+0.000059D0*BIGT)*BIGT 3114 3115 ZETA = (W+(0.30242D0-0.000269D0*BIGT+0.017996D0*T)*T)*TAS2R 3116 Z = (W+(1.09478D0+0.000387D0*BIGT+0.018324D0*T)*T)*TAS2R 3117 THETA = (2005.1125D0+(-0.85294D0-0.000365D0*BIGT)*BIGT+ 3118 : (-0.42647D0-0.000365D0*BIGT-0.041802D0*T)*T)*TAS2R 3119 3120* Rotation matrix 3121 CALL sla_DEULER('ZYZ',-ZETA,THETA,-Z,RMATP) 3122 3123 END 3124 SUBROUTINE sla_PRECES (SYSTEM, EP0, EP1, RA, DC) 3125*+ 3126* - - - - - - - 3127* P R E C E S 3128* - - - - - - - 3129* 3130* Precession - either FK4 (Bessel-Newcomb, pre IAU 1976) or 3131* FK5 (Fricke, post IAU 1976) as required. 3132* 3133* Given: 3134* SYSTEM char precession to be applied: 'FK4' or 'FK5' 3135* EP0,EP1 dp starting and ending epoch 3136* RA,DC dp RA,Dec, mean equator & equinox of epoch EP0 3137* 3138* Returned: 3139* RA,DC dp RA,Dec, mean equator & equinox of epoch EP1 3140* 3141* Called: sla_DRANRM, sla_PREBN, sla_PREC, sla_DCS2C, 3142* sla_DMXV, sla_DCC2S 3143* 3144* Notes: 3145* 3146* 1) Lowercase characters in SYSTEM are acceptable. 3147* 3148* 2) The epochs are Besselian if SYSTEM='FK4' and Julian if 'FK5'. 3149* For example, to precess coordinates in the old system from 3150* equinox 1900.0 to 1950.0 the call would be: 3151* CALL sla_PRECES ('FK4', 1900D0, 1950D0, RA, DC) 3152* 3153* 3) This routine will NOT correctly convert between the old and 3154* the new systems - for example conversion from B1950 to J2000. 3155* For these purposes see sla_FK425, sla_FK524, sla_FK45Z and 3156* sla_FK54Z. 3157* 3158* 4) If an invalid SYSTEM is supplied, values of -99D0,-99D0 will 3159* be returned for both RA and DC. 3160* 3161* P.T.Wallace Starlink 20 April 1990 3162* 3163* Copyright (C) 1995 Rutherford Appleton Laboratory 3164* 3165* License: 3166* This program is free software; you can redistribute it and/or modify 3167* it under the terms of the GNU General Public License as published by 3168* the Free Software Foundation; either version 2 of the License, or 3169* (at your option) any later version. 3170* 3171* This program is distributed in the hope that it will be useful, 3172* but WITHOUT ANY WARRANTY; without even the implied warranty of 3173* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 3174* GNU General Public License for more details. 3175* 3176* You should have received a copy of the GNU General Public License 3177* along with this program (see SLA_CONDITIONS); if not, write to the 3178* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 3179* Boston, MA 02110-1301 USA 3180* 3181*- 3182 3183 IMPLICIT NONE 3184 3185 CHARACTER SYSTEM*(*) 3186 DOUBLE PRECISION EP0,EP1,RA,DC 3187 3188 DOUBLE PRECISION PM(3,3),V1(3),V2(3) 3189 CHARACTER SYSUC*3 3190 3191 DOUBLE PRECISION sla_DRANRM 3192 3193 3194 3195 3196* Convert to uppercase and validate SYSTEM 3197 SYSUC=SYSTEM 3198 IF (SYSUC(1:1).EQ.'f') SYSUC(1:1)='F' 3199 IF (SYSUC(2:2).EQ.'k') SYSUC(2:2)='K' 3200 IF (SYSUC.NE.'FK4'.AND.SYSUC.NE.'FK5') THEN 3201 RA=-99D0 3202 DC=-99D0 3203 ELSE 3204 3205* Generate appropriate precession matrix 3206 IF (SYSUC.EQ.'FK4') THEN 3207 CALL sla_PREBN(EP0,EP1,PM) 3208 ELSE 3209 CALL sla_PREC(EP0,EP1,PM) 3210 END IF 3211 3212* Convert RA,Dec to x,y,z 3213 CALL sla_DCS2C(RA,DC,V1) 3214 3215* Precess 3216 CALL sla_DMXV(PM,V1,V2) 3217 3218* Back to RA,Dec 3219 CALL sla_DCC2S(V2,RA,DC) 3220 RA=sla_DRANRM(RA) 3221 3222 END IF 3223 3224 END 3225 SUBROUTINE sla_PREC (EP0, EP1, RMATP) 3226*+ 3227* - - - - - 3228* P R E C 3229* - - - - - 3230* 3231* Form the matrix of precession between two epochs (IAU 1976, FK5) 3232* (double precision) 3233* 3234* Given: 3235* EP0 dp beginning epoch 3236* EP1 dp ending epoch 3237* 3238* Returned: 3239* RMATP dp(3,3) precession matrix 3240* 3241* Notes: 3242* 3243* 1) The epochs are TDB (loosely ET) Julian epochs. 3244* 3245* 2) The matrix is in the sense V(EP1) = RMATP * V(EP0) 3246* 3247* 3) Though the matrix method itself is rigorous, the precession 3248* angles are expressed through canonical polynomials which are 3249* valid only for a limited time span. There are also known 3250* errors in the IAU precession rate. The absolute accuracy 3251* of the present formulation is better than 0.1 arcsec from 3252* 1960AD to 2040AD, better than 1 arcsec from 1640AD to 2360AD, 3253* and remains below 3 arcsec for the whole of the period 3254* 500BC to 3000AD. The errors exceed 10 arcsec outside the 3255* range 1200BC to 3900AD, exceed 100 arcsec outside 4200BC to 3256* 5600AD and exceed 1000 arcsec outside 6800BC to 8200AD. 3257* The SLALIB routine sla_PRECL implements a more elaborate 3258* model which is suitable for problems spanning several 3259* thousand years. 3260* 3261* References: 3262* Lieske,J.H., 1979. Astron.Astrophys.,73,282. 3263* equations (6) & (7), p283. 3264* Kaplan,G.H., 1981. USNO circular no. 163, pA2. 3265* 3266* Called: sla_DEULER 3267* 3268* P.T.Wallace Starlink 23 August 1996 3269* 3270* Copyright (C) 1996 Rutherford Appleton Laboratory 3271* 3272* License: 3273* This program is free software; you can redistribute it and/or modify 3274* it under the terms of the GNU General Public License as published by 3275* the Free Software Foundation; either version 2 of the License, or 3276* (at your option) any later version. 3277* 3278* This program is distributed in the hope that it will be useful, 3279* but WITHOUT ANY WARRANTY; without even the implied warranty of 3280* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 3281* GNU General Public License for more details. 3282* 3283* You should have received a copy of the GNU General Public License 3284* along with this program (see SLA_CONDITIONS); if not, write to the 3285* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 3286* Boston, MA 02110-1301 USA 3287* 3288*- 3289 3290 IMPLICIT NONE 3291 3292 DOUBLE PRECISION EP0,EP1,RMATP(3,3) 3293 3294* Arc seconds to radians 3295 DOUBLE PRECISION AS2R 3296 PARAMETER (AS2R=0.484813681109535994D-5) 3297 3298 DOUBLE PRECISION T0,T,TAS2R,W,ZETA,Z,THETA 3299 3300 3301 3302* Interval between basic epoch J2000.0 and beginning epoch (JC) 3303 T0 = (EP0-2000D0)/100D0 3304 3305* Interval over which precession required (JC) 3306 T = (EP1-EP0)/100D0 3307 3308* Euler angles 3309 TAS2R = T*AS2R 3310 W = 2306.2181D0+(1.39656D0-0.000139D0*T0)*T0 3311 3312 ZETA = (W+((0.30188D0-0.000344D0*T0)+0.017998D0*T)*T)*TAS2R 3313 Z = (W+((1.09468D0+0.000066D0*T0)+0.018203D0*T)*T)*TAS2R 3314 THETA = ((2004.3109D0+(-0.85330D0-0.000217D0*T0)*T0) 3315 : +((-0.42665D0-0.000217D0*T0)-0.041833D0*T)*T)*TAS2R 3316 3317* Rotation matrix 3318 CALL sla_DEULER('ZYZ',-ZETA,THETA,-Z,RMATP) 3319 3320 END 3321 SUBROUTINE sla_PVOBS (P, H, STL, PV) 3322*+ 3323* - - - - - - 3324* P V O B S 3325* - - - - - - 3326* 3327* Position and velocity of an observing station (double precision) 3328* 3329* Given: 3330* P dp latitude (geodetic, radians) 3331* H dp height above reference spheroid (geodetic, metres) 3332* STL dp local apparent sidereal time (radians) 3333* 3334* Returned: 3335* PV dp(6) position/velocity 6-vector (AU, AU/s, true equator 3336* and equinox of date) 3337* 3338* Called: sla_GEOC 3339* 3340* IAU 1976 constants are used. 3341* 3342* P.T.Wallace Starlink 14 November 1994 3343* 3344* Copyright (C) 1995 Rutherford Appleton Laboratory 3345* 3346* License: 3347* This program is free software; you can redistribute it and/or modify 3348* it under the terms of the GNU General Public License as published by 3349* the Free Software Foundation; either version 2 of the License, or 3350* (at your option) any later version. 3351* 3352* This program is distributed in the hope that it will be useful, 3353* but WITHOUT ANY WARRANTY; without even the implied warranty of 3354* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 3355* GNU General Public License for more details. 3356* 3357* You should have received a copy of the GNU General Public License 3358* along with this program (see SLA_CONDITIONS); if not, write to the 3359* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 3360* Boston, MA 02110-1301 USA 3361* 3362*- 3363 3364 IMPLICIT NONE 3365 3366 DOUBLE PRECISION P,H,STL,PV(6) 3367 3368 DOUBLE PRECISION R,Z,S,C,V 3369 3370* Mean sidereal rate (at J2000) in radians per (UT1) second 3371 DOUBLE PRECISION SR 3372 PARAMETER (SR=7.292115855306589D-5) 3373 3374 3375 3376* Geodetic to geocentric conversion 3377 CALL sla_GEOC(P,H,R,Z) 3378 3379* Functions of ST 3380 S=SIN(STL) 3381 C=COS(STL) 3382 3383* Speed 3384 V=SR*R 3385 3386* Position 3387 PV(1)=R*C 3388 PV(2)=R*S 3389 PV(3)=Z 3390 3391* Velocity 3392 PV(4)=-V*S 3393 PV(5)=V*C 3394 PV(6)=0D0 3395 3396 END 3397