program b50j20 c Transforms B1950 positions (i.e. on the equator, equinox and epoch of c B1950.0) to the equator, equinox and epoch of J2000. c References: c Smith, C.A., et. al., "Mean and Apparent Place Computations in the c New IAU System. I. The Transformation of Astrometric Catalog c Systems of the Equinox of J2000.0" Astron. J. 97, 265 c Seidelmann, P.K. editor, "Explanatory Supplement to the Astronomical c Almanac. Ch 3.5 Transformation to the FK5 System and Epoch J2000.0 c This program is based on a PL1 program written by Dr. C. Smith. c The output from this program has been compared with that from c his original PL1 and they are identical. The original program c was written for use on an IBM4381 using the IBM PL1 compiler. implicit none character namei*60 character nameo*60 integer*4 fk4no character misc1*21,misc2*1,decsgn*1 real*8 rah,ram,ras real*8 ddeg,dmin,dsec real*8 ra,dec,parlx,rapm,decpm,radvel real*8 dint real*8 > X50E,Y50E,Z50E,X2000,Y2000,Z2000,MUX2000,MUY2000,MUZ2000,RA50, > RA50E,DEC50,DEC50E,MURA,MUDEC,MUX,MUY,MUZ,SINRA50,SINRA50E, > COSRA50, > COSRA50E,SINDEC50,SINDEC50E,COSDEC50,COSDEC50E,MUSQ, > MUXDOT,MUYDOT,MUZDOT,SINZETA0,COSZETA0,SINZ,COSZ,SINTHETA, > COSTHETA, > ZETA0,Z,THETA,XX,YX,ZX,XY,YY,ZY,XZ,YZ,ZZ,XT,YT,ZT,RA2000, > DEC2000, > MURA2000,MUDEC2000,TRPTOJUL,SECPERRAD,MUX2,MUY2,MUZ2, > SINRA2000,COSRA2000,COSDEC2000,FORSHORT,NEWMURA84,NEWMUDEC84, > X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,R2000SQ,T0,T,DALPHA,DDELTA real*8 > R1984SQ,RA84,DEC84,MUX84,MUY84,MUZ84,MURA84,MUDEC84, > SINRA84,COSRA84,SINDEC84,COSDEC84,X84,Y84,Z84,NEWRA84,COSNEWRA84, > SINNEWRA84,NEWX84,NEWY84,NEWZ84 real*8 rass,dsecs c/**************************************************** c/* File names c/**************************************************** data namei /'/users/ftp/pub/fk4b1950.input'/ ! input data nameo /'/users/ftp/pub/b50j20.out '/ ! output c/*********************************************************************/ c/* SECPERRAD IS THE NUMBER OF SECONDS OF ARC PER RADIAN */ c/*********************************************************************/ SECPERRAD=57.29577951d0*3600.d0 c/*********************************************************************/ c/* TRPTOJUL IS THE NUMBER OF JULIAN CENTURIES PER TROPICAL CENTURY */ c/* AND IS TAKEN AS THE RATIO 36525/36524.2198781 */ c/*********************************************************************/ TRPTOJUL=36525.d0/36524.2198781d0 c/*********************************************************************/ c/* BY INTERNATIONAL AGREEMENT, 1984 JAN. 1.0 IS THE DATE ON WHICH */ c/* THE NEW SYSTEM OF IAU CONSTANTS IS TO GO INTO USE BY TIMEKEEPING */ c/* LABORATORIES. TO PREVENT A DISCONTINUITY IN UT, THE DEFINITION OF */ c/* A NEW RELATIONSHIP BETWEEN UNIVERSAL AND SIDEREAL TIME WENT INTO */ c/* EFFECT AT THE SAME TIME. IT IS TO PRESERVE THIS RELATIONSHIP THAT */ c/* 1984 JAN. 1.0 HAS SPECIAL SIGNIFICANCE AS THE DATE AT WHICH THE */ c/* DISCONTINUITY IN THE RIGHT ASCENSIONS AND THEIR PROPER MOTIONS ARE*/ c/* INTRODUCED, AND THE PRECESSION BASIS FOR THE PROPER MOTIONS IS */ c/* CHANGED FROM NEWCOMB'S TO THAT OF LIESKE, ET AL. */ c/*********************************************************************/ c/* PRECESSION FROM B1950.0 TO 1984 JAN 1.0, I.E., */ c/* JD2433282.4235 TO JD2445700.5 */ c/*********** *************/ c/* B1950.0 MINUS J2000 IN UNITS OF JULIAN CENTURIES */ c/*********************************************************************/ T0= (2433282.4235d0-2451545.0d0)/36525.d0 c/*********************************************************************/ c/* 1984 JAN 1.0 MINUS B1950.0 IN UNITS OF JULIAN CENTURIES */ c/*********************************************************************/ T=(2445700.5d0-2433282.4235d0)/36525.d0 c/*********************************************************************/ c/* THE FOLLOWING EQUATIONS ARE BASED ON NEWCOMB'S PRECESSION AND */ c/* HAVE BEEN ADAPTED BY ARI FROM ANDOYER'S DEVELOPMENT GIVEN IN */ c/* WOOLARD AND CLEMENCE - SPHERICAL ASTRONOMY, PG. 262 */ c/* THE UNIT OF TIME IS THE JULIAN CENTURY. THE CONSTANTS IN THE EQUA-*/ c/* TIONS ARE REFERRED TO THE EPOCH OF J2000. SEE THE ARTICLE BY LED- */ c/* ERLE AND SCHWAN (ASTRON. AND ASTROPH., VOL. 134, NO. 1 MAY 1984, */ c/* PG. 1 */ c/*********************************************************************/ ZETA0=(2305.6997d0+1.39744d0*T0+0.000060d0*T0**2)*T+ > (0.30201d0-0.000270d0*T0)*T**2+0.017996d0*T**3 Z=(2305.6997d0+1.39744d0*T0+0.000060d0*T0**2)*T+ > (1.09543d0+0.000390d0*T0)*T**2+0.018326d0*T**3 THETA=(2003.8746d0-0.85405d0*T0-0.000370d0*T0**2)*T+ > (-0.42707d0-0.000370d0*T0)*T**2-0.041803d0*T**3 write(6,*) T0,T,ZETA0,Z,THETA ZETA0=ZETA0/3600.d0 Z=Z/3600.d0 THETA=THETA/3600.d0 SINZETA0=dSIND(ZETA0) COSZETA0=dCOSD(ZETA0) SINZ=dSIND(Z) COSZ=dCOSD(Z) SINTHETA=dSIND(THETA) COSTHETA=dCOSD(THETA) X1=COSZETA0*COSTHETA*COSZ-SINZETA0*SINZ Y1=-SINZETA0*COSTHETA*COSZ-COSZETA0*SINZ Z1=-SINTHETA*COSZ X2=COSZETA0*COSTHETA*SINZ+SINZETA0*COSZ Y2=-SINZETA0*COSTHETA*SINZ+COSZETA0*COSZ Z2=-SINTHETA*SINZ X3=COSZETA0*SINTHETA Y3=-SINZETA0*SINTHETA Z3=COSTHETA c/*********************************************************************/ c/* PRECESSION FROM 1984 JAN 1.0 TO J2000 (JD2445700.5 TO */ c/* JD2451545.0) */ c/* 1984 JAN 1.0 MINUS J2000 IN UNITS OF JULIAN CENTURIES */ c/*********************************************************************/ T0=(2445700.5d0-2451545.0d0)/36525.0d0 c/*********************************************************************/ c/* 2000 JAN 1.5 (J2000) MINUS 1984 JAN 1.0 IN JULIAN CENTURIES */ c/*********************************************************************/ T=(2451545.0d0-2445700.5d0)/36525.0d0 ZETA0=(2306.2181d0+1.39656d0*T0-0.000139d0*T0**2)*T+ > (0.30188d0-0.000344d0*T0)*T**2+0.017998d0*T**3 Z=(2306.2181d0+1.39656d0*T0-0.000139d0*T0**2)*T+ > (1.09468d0+0.000066d0*T0)*T**2+0.018203d0*T**3 THETA=(2004.3109d0-0.85330d0*T0-0.000217d0*T0**2)*T+ > (-0.42665d0-0.000217d0*T0)*T**2-0.041833d0*T**3 write(6,*) T0,T,ZETA0,Z,THETA ZETA0=ZETA0/3600.d0 Z=Z/3600.d0 THETA=THETA/3600.d0 SINZETA0=dSIND(ZETA0) COSZETA0=dCOSD(ZETA0) SINZ=dSIND(Z) COSZ=dCOSD(Z) SINTHETA=dSIND(THETA) COSTHETA=dCOSD(THETA) XX=COSZETA0*COSTHETA*COSZ-SINZETA0*SINZ YX=-SINZETA0*COSTHETA*COSZ-COSZETA0*SINZ ZX=-SINTHETA*COSZ XY=COSZETA0*COSTHETA*SINZ+SINZETA0*COSZ YY=-SINZETA0*COSTHETA*SINZ+COSZETA0*COSZ ZY=-SINTHETA*SINZ XZ=COSZETA0*SINTHETA YZ=-SINZETA0*SINTHETA ZZ=COSTHETA write(6,*) X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3 write(6,*) XX,YX,ZX,XY,YY,ZY,XZ,YZ,ZZ open(501,file=namei,status='OLD') open(502,file=nameo,status='UNKNOWN') 100 continue read(501,200,end=900) fk4no,misc1,rah,ram,ras,misc2,decsgn, > ddeg,dmin,dsec,ra,dec,parlx,rapm,decpm,radvel 200 format (i4,bz,a21,i2,i2,i5,a1,a1,i2,i2,i4,f13.10,f14.10,i4, > 2i7,f6.1) parlx=parlx*1.d-3 RA50=RAS RA50=RA50/1000. RA50=RA50+RAM*60.0+RAH*3600.0 DEC50=DSEC dec50=dec50/100. DEC50=DEC50+DMIN*60.0+DDEG*3600.0 IF(DECSGN.eq.'-') DEC50=-DEC50 c/*********************************************************************/ c/* CALCULATE DIRECTION COSINES */ c/*********************************************************************/ SINRA50=dSIND(RA50/240.d0) COSRA50=dCOSD(RA50/240.d0) SINDEC50=dSIND(DEC50/3600.d0) COSDEC50=dCOSD(DEC50/3600.d0) c/*********************************************************************/ c/* CALCULATE THE CORRECTION FOR THE ELLIPTIC TERMS IN ABERRATION */ c/*********************************************************************/ DALPHA=(0.065838d0*COSRA50-0.335299d0*SINRA50)/(15d0*COSDEC50) DDELTA=(-0.335299d0*COSRA50-0.065838d0*SINRA50)*SINDEC50+ > 0.028553d0*COSDEC50 RA50E=RA50+DALPHA DEC50E=DEC50+DDELTA c/*********************************************************************/ c/* CALCULATE DIRECTION COSINES AND ADVANCE THEM TO THE EPOCH OF 1984 */ c/*********************************************************************/ SINRA50E=dSIND(RA50E/240.d0) COSRA50E=dCOSD(RA50E/240.d0) SINDEC50E=dSIND(DEC50E/3600.d0) COSDEC50E=dCOSD(DEC50E/3600.d0) X50E=COSRA50E*COSDEC50E Y50E=SINRA50E*COSDEC50E Z50E=SINDEC50E MURA=RAPM MURA=MURA*1.d-3 MURA=MURA*15.d0/SECPERRAD MUDEC=DECPM MUDEC=MUDEC*1.d-2/SECPERRAD MUX=-MURA*Y50E-MUDEC*COSRA50E*SINDEC50E MUY=+MURA*X50E-MUDEC*SINRA50E*SINDEC50E MUZ=+MUDEC*COSDEC50E MUSQ=MUX**2+MUY**2+MUZ**2 FORSHORT=PARLX FORSHORT=FORSHORT*RADVEL*0.0002045381d0 MUXDOT=-X50E*MUSQ-FORSHORT*MUX MUYDOT=-Y50E*MUSQ-FORSHORT*MUY MUZDOT=-Z50E*MUSQ-FORSHORT*MUZ T=0.339995667d0 ! /************** UNIT = TROPICAL CENTURIES **********/ XT=X50E+(MUX+MUXDOT*T/2.d0)*T YT=Y50E+(MUY+MUYDOT*T/2.d0)*T ZT=Z50E+(MUZ+MUZDOT*T/2.d0)*T X84=X1*XT+Y1*YT+Z1*ZT Y84=X2*XT+Y2*YT+Z2*ZT Z84=X3*XT+Y3*YT+Z3*ZT R1984SQ=X84**2+Y84**2+Z84**2 RA84=dATAN2(Y84,X84) RA84=RA84*SECPERRAD/15.d0 DEC84=dATAN2(Z84,dSQRT(X84**2+Y84**2)) DEC84=DEC84*SECPERRAD IF(dABS(DEC84).ge.324000.d0) write(6,'("ABS(DEC1984) EXCEEDS 90" > " DEGREES - FK4NO IS",i5)') FK4NO c/*********************************************************************/ c/* ADVANCE THE PROPER MOTION DIRECTION COSINES TO THE EPOCH OF 1984 */ c/*********************************************************************/ MUX2=MUX+MUXDOT*T MUY2=MUY+MUYDOT*T MUZ2=MUZ+MUZDOT*T MUX84=X1*MUX2+Y1*MUY2+Z1*MUZ2 MUY84=X2*MUX2+Y2*MUY2+Z2*MUZ2 MUZ84=X3*MUX2+Y3*MUY2+Z3*MUZ2 MURA84=(MUY84*X84-MUX84*Y84)/(X84**2+Y84**2) MUDEC84=(MUZ84*R1984SQ-Z84*(X84*MUX84+Y84*MUY84+Z84*MUZ84)) > /(R1984SQ*SQRT(R1984SQ-Z84**2)) c/*********************************************************************/ c/* CORRECT THE 1984 RIGHT ASCENSION FOR ZERO POINT ERROR */ c/*********************************************************************/ NEWRA84=RA84+0.0639d0 c/*********************************************************************/ c/* CALCULATE DIRECTION COSINES OF THE CORRECTED RIGHT ASCENSION AND */ c/* ORIGINAL DECLINATION */ c/*********************************************************************/ SINRA84=dSIND(RA84/240.d0) SINNEWRA84=dSIND(NEWRA84/240.d0) COSRA84=dCOSD(RA84/240.d0) COSNEWRA84=dCOSD(NEWRA84/240.d0) SINDEC84=dSIND(DEC84/3600.d0) COSDEC84=dCOSD(DEC84/3600.d0) NEWX84=COSNEWRA84*COSDEC84 NEWY84=SINNEWRA84*COSDEC84 NEWZ84=SINDEC84 c/*********************************************************************/ c/* CONVERT UNIT OF TIME OF PROPER MOTION FROM TROPICAL TO JULIAN */ c/* CENTURIES , CHANGE FROM THE NEWCOMB TO THE LIESKE, ET AL. BASIS, */ c/* AND APPLY THE ZERO POINT CORRECTION TO THE RT. ASC. PROPER MOTION */ c/*********************************************************************/ MURA84=MURA84*TRPTOJUL*SECPERRAD/15.d0 NEWMURA84=MURA84+0.085d0-0.069138d0- > (133.629829d0*SINNEWRA84-133.600750d0*SINRA84)*SINDEC84/COSDEC84 MUDEC84=MUDEC84*TRPTOJUL*SECPERRAD NEWMUDEC84=MUDEC84- > (2004.447434d0*COSNEWRA84-2004.011250d0*COSRA84) MURA=NEWMURA84*15.d0/SECPERRAD MUDEC=NEWMUDEC84/SECPERRAD MUX=-MURA*NEWY84-MUDEC*COSNEWRA84*SINDEC84 MUY=+MURA*NEWX84-MUDEC*SINNEWRA84*SINDEC84 MUZ=MUDEC*COSDEC84 MUSQ=MUX**2+MUY**2+MUZ**2 FORSHORT=PARLX FORSHORT=FORSHORT*RADVEL*0.0002045381d0 MUXDOT=-NEWX84*MUSQ-FORSHORT*MUX MUYDOT=-NEWY84*MUSQ-FORSHORT*MUY MUZDOT=-NEWZ84*MUSQ-FORSHORT*MUZ c/*********************************************************************/ c/* ADVANCE THE DIRECTION COSINES FROM 1984 TO THE EPOCH AND EQUINOX */ c/* OF J2000 */ c/*********************************************************************/ T=0.1600136893d0 !/***************** UNIT = JULIAN CENTURIES ********/ XT=NEWX84+(MUX+MUXDOT*T/2.d0)*T YT=NEWY84+(MUY+MUYDOT*T/2.d0)*T ZT=NEWZ84+(MUZ+MUZDOT*T/2.d0)*T X2000=XX*XT+YX*YT+ZX*ZT Y2000=XY*XT+YY*YT+ZY*ZT Z2000=XZ*XT+YZ*YT+ZZ*ZT R2000SQ=X2000**2+Y2000**2+Z2000**2 RA2000=dATAN2(Y2000,X2000) COSRA2000=dCOS(RA2000) SINRA2000=dSIN(RA2000) RA2000=RA2000*SECPERRAD/15.d0 DEC2000=dATAN2(Z2000,dSQRT(X2000**2+Y2000**2)) COSDEC2000=dCOS(DEC2000) DEC2000=DEC2000*SECPERRAD IF(dABS(DEC2000).ge.324000.d0) write(6,'("ABS(DEC2000) EXCEEDS 90" > " DEGREES - FK4NO IS",i5)') FK4NO c/*********************************************************************/ c/* ADVANCE THE PROPER MOTIONS TO THE EQUINOX AND EPOCH OF J2000 */ c/*********************************************************************/ MUX2=MUX+MUXDOT*T MUY2=MUY+MUYDOT*T MUZ2=MUZ+MUZDOT*T MUX2000=XX*MUX2+YX*MUY2+ZX*MUZ2 MUY2000=XY*MUX2+YY*MUY2+ZY*MUZ2 MUZ2000=XZ*MUX2+YZ*MUY2+ZZ*MUZ2 MURA2000= > (((MUY2000*X2000-MUX2000*Y2000)/(X2000**2+Y2000**2)) > *SECPERRAD)/15.d0 MUDEC2000= > ((MUZ2000*R2000SQ-Z2000*(X2000*MUX2000+Y2000*MUY2000+ > Z2000*MUZ2000)) > /(R2000SQ*SQRT(R2000SQ-Z2000**2)))*SECPERRAD IF( RA2000 .lt. 0.d0) RA2000=RA2000+86400.d0 RA2000=RA2000+0.0005d0 RAH=dint(RA2000/3600.0) RAM=dint((RA2000-RAH*3600.0)/60.0) rass=RA2000-RAH*3600.0-RAM*60.0 IF(RAH.ge.24) RAH=RAH-24 RA2000=rass RA2000=RA2000+RAM*60.0+RAH*3600.0 RA=RA2000/SECPERRAD*15.d0+5.d-11 IF(DEC2000.lt.0.d0) THEN DECSGN='-' ELSE DECSGN='+' ENDIF DEC2000=ABS(DEC2000)+0.005d0 DDEG=dint(DEC2000/3600.0) DMIN=dint((DEC2000-DDEG*3600.0)/60.0) dsecs=DEC2000-DDEG*3600.0-DMIN*60.0 DEC2000=dsecs DEC2000=DEC2000+DDEG*3600.0+DMIN*60.0 DEC=DEC2000/SECPERRAD+5.d-11 IF(DECSGN.eq.'-') DEC=-DEC dsec= dsecs*100. ras = rass*1000. RAPM=(MURA2000 + dSIGN(5.d-4,MURA2000) )*1000. DECPM=(MUDEC2000 + dSIGN(5.d-3,MUDEC2000) )*100. parlx=parlx*1.d+3 write(502,200) fk4no,misc1,rah,ram,ras ,misc2,decsgn, > ddeg,dmin,dsec,ra,dec,parlx,rapm,decpm,radvel GO TO 100 900 stop 'finished' e n d