program koko c c Coordinate conversion program that uses command line inputs... c implicit none double precision rab,deb,raj,dej,l,b,pi,strad,sarad, & DR1950,DD1950,P1950,V1950,DR2000,DD2000,P2000,V2000, & rah,ram,ras,ded,dem,des,ra,de,dum,al,be integer maxpar,npar parameter(pi=3.1415927,maxpar=8) character*80 line,cpostn,chra,chde,pars(maxpar),fname integer narg, iargc, isys DR1950=0.0 DD1950=0.0 P1950=0.0 V1950=0.0 DR2000=0.0 DD2000=0.0 P2000=0.0 V2000=0.0 strad = 2.0*pi/86400.0 sarad = strad/15.0 fname = ' ' narg=iargc() if (narg.eq.2) then call getarg(2,fname) else if (narg.ne.3) then stop 'koko [J/B/G] [x] [y]' endif call getarg(1,line) if (line.eq.'G'.or.line.eq.'g') then write(*,*) 'Input system: IAU 1958 Galactic coords' isys=1 else if (line.eq.'B'.or.line.eq.'b') then write(*,*) 'Input system: FK4 B1950 Equatorial coords' isys=2 else if (line.eq.'J'.or.line.eq.'j') then write(*,*) 'Input system: FK5 J2000 Equatorial coords' isys=3 else stop 'Please specify J, B or G as the first argument!' endif if (fname.ne.' ') open(unit=10,file=fname,status='old') 5 if (fname.ne.' ') read(10,'(a)',end=10) line if (isys.eq.1) then if (fname.eq.' ') then call getarg(2,line) read(line,*) l call getarg(3,line) read(line,*) b else read(line,*) dum,l,b endif l=l*pi/180.0 b=b*pi/180.0 call sla_ge50(l,b,rab,deb) call sla_galeq(l,b,raj,dej) else if (fname.eq.' ') then call getarg(2,line) if (index(line,':').gt.0) then call parse(line,pars,npar,maxpar,':') read(pars(1),*) rah read(pars(2),*) ram read(pars(3),*) ras else read(line(1:2),*) rah read(line(3:4),*) ram read(line(5:),*) ras endif ra=rah+ram/60.0+ras/3600.0 ra=ra*pi/12.0 call getarg(3,line) call parse(line,pars,npar,maxpar,':') if (index(line,':').gt.0) then read(pars(1),*) ded read(pars(2),*) dem read(pars(3),*) des else if (line(1:1).eq.'+'.or.line(1:1).eq.'-') then read(line(1:3),*) ded read(line(4:5),*) dem read(line(6:),*) des else read(line(1:2),*) ded read(line(3:4),*) dem read(line(5:),*) des endif endif if (ded.lt.0.0) de=ded-dem/60.0-des/3600.0 if (ded.ge.0.0) de=ded+dem/60.0+des/3600.0 de=de*pi/180.0 else read(line,*) dum,ra,de ra=ra*pi/180.0 de=de*pi/180.0 endif if (isys.eq.2) then rab=ra deb=de call sla_eg50(rab,deb,l,b) call sla_fk425(rab,deb,DR1950,DD1950,P1950,V1950, & raj,dej,DR2000,DD2000,P2000,V2000) else raj=ra dej=de call sla_eqgal(raj,dej,l,b) call sla_fk524(raj,dej,DR2000,DD2000,P2000,V2000, & rab,deb,DR1950,DD1950,P1950,V1950) endif endif al=rab*180.0/3.1415927 be=deb*180.0/3.1415927 rab=rab/strad deb=deb/sarad chra=cpostn(rab,2,1) chde=cpostn(deb,2,1) if (fname.eq.' ') then write(*,*) 'B1950 FK4' write(*,'(a,x,a)') chra(1:index(chra,' ')-1), & chde(1:index(chde,' ')-1) write(*,*) al,be endif al=raj*180.0/3.1415927 be=dej*180.0/3.1415927 raj=raj/strad dej=dej/sarad chra=cpostn(raj,2,1) chde=cpostn(dej,2,1) if (fname.eq.' ') then write(*,*) 'J2000 FK5' write(*,'(a,x,a)') chra(1:index(chra,' ')-1), & chde(1:index(chde,' ')-1) write(*,*) al,be endif l=l*180.0/pi b=b*180.0/pi if (fname.eq.' ') then write(*,*) 'IAU 1958 Galactic coords' write(*,*) 'l:',float(l),' deg' write(*,*) 'b:',float(b),' deg' endif if (fname.eq.' ') stop write(*,*) int(dum),real(l),real(b) c & ,real(rab*strad),real(deb*sarad), c & real(raj*strad),real(dej*sarad) goto 5 10 close(unit=10) end SUBROUTINE sla_FK524 (R2000,D2000,DR2000,DD2000,P2000,V2000, : R1950,D1950,DR1950,DD1950,P1950,V1950) *+ * - - - - - - * F K 5 2 4 * - - - - - - * * Convert J2000.0 FK5 star data to B1950.0 FK4 (double precision) * * This routine converts stars from the new, IAU 1976, FK5, Fricke * system, to the old, Bessel-Newcomb, FK4 system, using Yallop's * implementation (see ref 2) of a matrix method due to Standish. * The numerical values of ref 2 are used canonically. * * Given: (all J2000.0,FK5) * R2000,D2000 dp J2000.0 RA,Dec (rad) * DR2000,DD2000 dp J2000.0 proper motions (rad/Jul.yr) * P2000 dp parallax (arcsec) * V2000 dp radial velocity (km/s, +ve = moving away) * * Returned: (all B1950.0,FK4) * R1950,D1950 dp B1950.0 RA,Dec (rad) * DR1950,DD1950 dp B1950.0 proper motions (rad/trop.yr) * P1950 dp parallax (arcsec) * V1950 dp radial velocity (km/s, +ve = moving away) * * Notes: * * 1) The proper motions in RA are dRA/dt rather than * cos(Dec)*dRA/dt, and are per year rather than per century. * * 2) Note that conversion from Julian epoch 2000.0 to Besselian * epoch 1950.0 only is provided for. Conversions involving * other epochs will require use of the appropriate precession, * proper motion, and E-terms routines before and/or after * FK524 is called. * * 3) In the FK4 catalogue the proper motions of stars within * 10 degrees of the poles do not embody the differential * E-term effect and should, strictly speaking, be handled * in a different manner from stars outside these regions. * However, given the general lack of homogeneity of the star * data available for routine astrometry, the difficulties of * handling positions that may have been determined from * astrometric fields spanning the polar and non-polar regions, * the likelihood that the differential E-terms effect was not * taken into account when allowing for proper motion in past * astrometry, and the undesirability of a discontinuity in * the algorithm, the decision has been made in this routine to * include the effect of differential E-terms on the proper * motions for all stars, whether polar or not. At epoch 2000, * and measuring on the sky rather than in terms of dRA, the * errors resulting from this simplification are less than * 1 milliarcsecond in position and 1 milliarcsecond per * century in proper motion. * * References: * * 1 Smith, C.A. et al, 1989. "The transformation of astrometric * catalog systems to the equinox J2000.0". Astron.J. 97, 265. * * 2 Yallop, B.D. et al, 1989. "Transformation of mean star places * from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space". * Astron.J. 97, 274. * * P.T.Wallace Starlink 22 November 1989 *- IMPLICIT NONE DOUBLE PRECISION R2000,D2000,DR2000,DD2000,P2000,V2000, : R1950,D1950,DR1950,DD1950,P1950,V1950 * Miscellaneous DOUBLE PRECISION R,D,UR,UD,PX,RV DOUBLE PRECISION SR,CR,SD,CD,X,Y,Z,W DOUBLE PRECISION V1(6),V2(6) DOUBLE PRECISION XD,YD,ZD DOUBLE PRECISION RXYZ,WD,RXYSQ,RXY INTEGER I,J * 2Pi DOUBLE PRECISION D2PI PARAMETER (D2PI=6.283185307179586476925287D0) * Radians per year to arcsec per century DOUBLE PRECISION PMF PARAMETER (PMF=100D0*60D0*60D0*360D0/D2PI) * Small number to avoid arithmetic problems DOUBLE PRECISION TINY PARAMETER (TINY=1D-30) * * CANONICAL CONSTANTS (see references) * * Km per sec to AU per tropical century * = 86400 * 36524.2198782 / 1.49597870E8 DOUBLE PRECISION VF PARAMETER (VF=21.0945D0) * Constant vector and matrix (by columns) DOUBLE PRECISION A(6),EMI(6,6) DATA A/ -1.62557D-6, -0.31919D-6, -0.13843D-6, : +1.245D-3, -1.580D-3, -0.659D-3/ DATA (EMI(I,1),I=1,6) / +0.999925679499910D0, : -0.011181482788805D0, : -0.004859004008828D0, : -0.000541640798032D0, : -0.237963047085011D0, : +0.436218238658637D0 / DATA (EMI(I,2),I=1,6) / +0.011181482840782D0, : +0.999937484898031D0, : -0.000027155744957D0, : +0.237912530551179D0, : -0.002660706488970D0, : -0.008537588719453D0 / DATA (EMI(I,3),I=1,6) / +0.004859003889183D0, : -0.000027177143501D0, : +0.999988194601879D0, : -0.436101961325347D0, : +0.012258830424865D0, : +0.002119065556992D0 / DATA (EMI(I,4),I=1,6) / -0.000002423898405D0, : +0.000000027105439D0, : +0.000000011777422D0, : +0.999904322043106D0, : -0.011181451601069D0, : -0.004858519608686D0 / DATA (EMI(I,5),I=1,6) / -0.000000027105439D0, : -0.000002423927017D0, : +0.000000000065851D0, : +0.011181451608968D0, : +0.999916125340107D0, : -0.000027162614355D0 / DATA (EMI(I,6),I=1,6) / -0.000000011777422D0, : +0.000000000065846D0, : -0.000002424049954D0, : +0.004858519590501D0, : -0.000027165866691D0, : +0.999966838131419D0 / * The above values were obtained by inverting C.Hohenkerk's * forward matrix (private communication), which agrees with * the one given in reference 2 but which has one additional * decimal place. * Pick up J2000 data (units radians and arcsec/JC) R=R2000 D=D2000 UR=DR2000*PMF UD=DD2000*PMF PX=P2000 RV=V2000 * Spherical to Cartesian SR=SIN(R) CR=COS(R) SD=SIN(D) CD=COS(D) X=CR*CD Y=SR*CD Z= SD W=VF*RV*PX V1(1)=X V1(2)=Y V1(3)=Z V1(4)=-UR*Y-CR*SD*UD+W*X V1(5)= UR*X-SR*SD*UD+W*Y V1(6)= CD*UD+W*Z * Convert position+velocity vector to BN system DO I=1,6 W=0D0 DO J=1,6 W=W+EMI(I,J)*V1(J) END DO V2(I)=W END DO * Vector components X=V2(1) Y=V2(2) Z=V2(3) XD=V2(4) YD=V2(5) ZD=V2(6) * Magnitude of position vector RXYZ=SQRT(X*X+Y*Y+Z*Z) * Radial velocity and parallax IF (PX.GT.TINY) THEN RV=(X*XD+Y*YD+Z*ZD)/(PX*VF*RXYZ) PX=PX/RXYZ END IF * Include E-terms W=X*A(1)+Y*A(2)+Z*A(3) WD=X*A(4)+Y*A(5)+Z*A(6) X=X+A(1)*RXYZ-W*X Y=Y+A(2)*RXYZ-W*Y Z=Z+A(3)*RXYZ-W*Z XD=XD+A(4)*RXYZ-WD*X YD=YD+A(5)*RXYZ-WD*Y ZD=ZD+A(6)*RXYZ-WD*Z * Convert to spherical RXYSQ=X*X+Y*Y RXY=SQRT(RXYSQ) IF (X.EQ.0D0.AND.Y.EQ.0D0) THEN R=0D0 ELSE R=ATAN2(Y,X) IF (R.LT.0.0D0) R=R+D2PI END IF D=ATAN2(Z,RXY) IF (RXY.GT.TINY) THEN UR=(X*YD-Y*XD)/RXYSQ UD=(ZD*RXYSQ-Z*(X*XD+Y*YD))/((RXYSQ+Z*Z)*RXY) END IF * Return results R1950=R D1950=D DR1950=UR/PMF DD1950=UD/PMF P1950=PX V1950=RV END SUBROUTINE sla_EQGAL (DR, DD, DL, DB) *+ * - - - - - - * E Q G A L * - - - - - - * * Transformation from J2000.0 equatorial coordinates to * IAU 1958 galactic coordinates (double precision) * * Given: * DR,DD dp J2000.0 RA,Dec * * Returned: * DL,DB dp galactic longitude and latitude L2,B2 * * (all arguments are radians) * * Called: * sla_DCS2C, sla_DMXV, sla_DCC2S, sla_DRANRM, sla_DRANGE * * Note: * The equatorial coordinates are J2000.0. Use the routine * sla_EG50 if conversion from B1950.0 'FK4' coordinates is * required. * * Reference: * Blaauw et al, Mon.Not.R.Astron.Soc.,121,123 (1960) * * P.T.Wallace Starlink November 1988 *- IMPLICIT NONE DOUBLE PRECISION DR,DD,DL,DB DOUBLE PRECISION sla_DRANRM,sla_DRANGE DOUBLE PRECISION V1(3),V2(3) * * L2,B2 system of galactic coordinates * * P = 192.25 RA of galactic north pole (mean B1950.0) * Q = 62.6 inclination of galactic to mean B1950.0 equator * R = 33 longitude of ascending node * * P,Q,R are degrees * * Equatorial to galactic rotation matrix (J2000.0), obtained by * applying the standard FK4 to FK5 transformation, for inertially * zero proper motion, to the columns of the B1950 equatorial to * galactic rotation matrix: * DOUBLE PRECISION RMAT(3,3) DATA RMAT(1,1),RMAT(1,2),RMAT(1,3), : RMAT(2,1),RMAT(2,2),RMAT(2,3), : RMAT(3,1),RMAT(3,2),RMAT(3,3)/ : -0.054875539726D0,-0.873437108010D0,-0.483834985808D0, : +0.494109453312D0,-0.444829589425D0,+0.746982251810D0, : -0.867666135858D0,-0.198076386122D0,+0.455983795705D0/ * Spherical to Cartesian CALL sla_DCS2C(DR,DD,V1) * Equatorial to galactic CALL sla_DMXV(RMAT,V1,V2) * Cartesian to spherical CALL sla_DCC2S(V2,DL,DB) * Express in conventional ranges DL=sla_DRANRM(DL) DB=sla_DRANGE(DB) END SUBROUTINE sla_GALEQ (DL, DB, DR, DD) *+ * - - - - - - * G A L E Q * - - - - - - * * Transformation from IAU 1958 galactic coordinates to * J2000.0 equatorial coordinates (double precision) * * Given: * DL,DB dp galactic longitude and latitude L2,B2 * * Returned: * DR,DD dp J2000.0 RA,Dec * * (all arguments are radians) * * Called: * sla_DCS2C, sla_DIMXV, sla_DCC2S, sla_DRANRM, sla_DRANGE * * Note: * The equatorial coordinates are J2000.0. Use the routine * sla_GE50 if conversion to B1950.0 'FK4' coordinates is * required. * * Reference: * Blaauw et al, Mon.Not.R.Astron.Soc.,121,123 (1960) * * P.T.Wallace Starlink November 1988 *- IMPLICIT NONE DOUBLE PRECISION DL,DB,DR,DD DOUBLE PRECISION sla_DRANRM,sla_DRANGE DOUBLE PRECISION V1(3),V2(3) * * L2,B2 system of galactic coordinates * * P = 192.25 RA of galactic north pole (mean B1950.0) * Q = 62.6 inclination of galactic to mean B1950.0 equator * R = 33 longitude of ascending node * * P,Q,R are degrees * * Equatorial to galactic rotation matrix (J2000.0), obtained by * applying the standard FK4 to FK5 transformation, for inertially * zero proper motion, to the columns of the B1950 equatorial to * galactic rotation matrix: * DOUBLE PRECISION RMAT(3,3) DATA RMAT(1,1),RMAT(1,2),RMAT(1,3), : RMAT(2,1),RMAT(2,2),RMAT(2,3), : RMAT(3,1),RMAT(3,2),RMAT(3,3)/ : -0.054875539726D0,-0.873437108010D0,-0.483834985808D0, : +0.494109453312D0,-0.444829589425D0,+0.746982251810D0, : -0.867666135858D0,-0.198076386122D0,+0.455983795705D0/ * Spherical to Cartesian CALL sla_DCS2C(DL,DB,V1) * Galactic to equatorial CALL sla_DIMXV(RMAT,V1,V2) * Cartesian to spherical CALL sla_DCC2S(V2,DR,DD) * Express in conventional ranges DR=sla_DRANRM(DR) DD=sla_DRANGE(DD) END SUBROUTINE sla_FK425 (R1950,D1950,DR1950,DD1950,P1950,V1950, : R2000,D2000,DR2000,DD2000,P2000,V2000) *+ * - - - - - - * F K 4 2 5 * - - - - - - * * Convert B1950.0 FK4 star data to J2000.0 FK5 (double precision) * * This routine converts stars from the old, Bessel-Newcomb, FK4 * system to the new, IAU 1976, FK5, Fricke system, using Yallop's * implementation (see ref 2) of a matrix method due to Standish. * The numerical values of ref 2 are used canonically. * * Given: (all B1950.0,FK4) * R1950,D1950 dp B1950.0 RA,Dec (rad) * DR1950,DD1950 dp B1950.0 proper motions (rad/trop.yr) * P1950 dp parallax (arcsec) * V1950 dp radial velocity (km/s, +ve = moving away) * * Returned: (all J2000.0,FK5) * R2000,D2000 dp J2000.0 RA,Dec (rad) * DR2000,DD2000 dp J2000.0 proper motions (rad/Jul.yr) * P2000 dp parallax (arcsec) * V2000 dp radial velocity (km/s, +ve = moving away) * * Notes: * * 1) The proper motions in RA are dRA/dt rather than * cos(Dec)*dRA/dt, and are per year rather than per century. * * 2) Conversion from Besselian epoch 1950.0 to Julian epoch * 2000.0 only is provided for. Conversions involving other * epochs will require use of the appropriate precession, * proper motion, and E-terms routines before and/or * after FK425 is called. * * 3) In the FK4 catalogue the proper motions of stars within * 10 degrees of the poles do not embody the differential * E-term effect and should, strictly speaking, be handled * in a different manner from stars outside these regions. * However, given the general lack of homogeneity of the star * data available for routine astrometry, the difficulties of * handling positions that may have been determined from * astrometric fields spanning the polar and non-polar regions, * the likelihood that the differential E-terms effect was not * taken into account when allowing for proper motion in past * astrometry, and the undesirability of a discontinuity in * the algorithm, the decision has been made in this routine to * include the effect of differential E-terms on the proper * motions for all stars, whether polar or not. At epoch 2000, * and measuring on the sky rather than in terms of dRA, the * errors resulting from this simplification are less than * 1 milliarcsecond in position and 1 milliarcsecond per * century in proper motion. * * References: * * 1 Smith, C.A. et al, 1989. "The transformation of astrometric * catalog systems to the equinox J2000.0". Astron.J. 97, 265. * * 2 Yallop, B.D. et al, 1989. "Transformation of mean star places * from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space". * Astron.J. 97, 274. * * P.T.Wallace Starlink 10 April 1990 *- IMPLICIT NONE DOUBLE PRECISION R1950,D1950,DR1950,DD1950,P1950,V1950, : R2000,D2000,DR2000,DD2000,P2000,V2000 * Miscellaneous DOUBLE PRECISION R,D,UR,UD,PX,RV,SR,CR,SD,CD,W,WD DOUBLE PRECISION X,Y,Z,XD,YD,ZD DOUBLE PRECISION RXYSQ,RXYZSQ,RXY,RXYZ,SPXY,SPXYZ INTEGER I,J * Star position and velocity vectors DOUBLE PRECISION R0(3),RD0(3) * Combined position and velocity vectors DOUBLE PRECISION V1(6),V2(6) * 2Pi DOUBLE PRECISION D2PI PARAMETER (D2PI=6.283185307179586476925287D0) * Radians per year to arcsec per century DOUBLE PRECISION PMF PARAMETER (PMF=100D0*60D0*60D0*360D0/D2PI) * Small number to avoid arithmetic problems DOUBLE PRECISION TINY PARAMETER (TINY=1D-30) * * CANONICAL CONSTANTS (see references) * * Km per sec to AU per tropical century * = 86400 * 36524.2198782 / 1.49597870E8 DOUBLE PRECISION VF PARAMETER (VF=21.0945D0) * Constant vector and matrix (by columns) DOUBLE PRECISION A(3),AD(3),EM(6,6) DATA A,AD/ -1.62557D-6, -0.31919D-6, -0.13843D-6, : +1.245D-3, -1.580D-3, -0.659D-3/ DATA (EM(I,1),I=1,6) / +0.999925678186902D0, : +0.011182059571766D0, : +0.004857946721186D0, : -0.000541652366951D0, : +0.237917612131583D0, : -0.436111276039270D0 / DATA (EM(I,2),I=1,6) / -0.011182059642247D0, : +0.999937478448132D0, : -0.000027147426498D0, : -0.237968129744288D0, : -0.002660763319071D0, : +0.012259092261564D0 / DATA (EM(I,3),I=1,6) / -0.004857946558960D0, : -0.000027176441185D0, : +0.999988199738770D0, : +0.436227555856097D0, : -0.008537771074048D0, : +0.002119110818172D0 / DATA (EM(I,4),I=1,6) / +0.000002423950176D0, : +0.000000027106627D0, : +0.000000011776559D0, : +0.999947035154614D0, : +0.011182506007242D0, : +0.004857669948650D0 / DATA (EM(I,5),I=1,6) / -0.000000027106627D0, : +0.000002423978783D0, : -0.000000000065816D0, : -0.011182506121805D0, : +0.999958833818833D0, : -0.000027137309539D0 / DATA (EM(I,6),I=1,6) / -0.000000011776558D0, : -0.000000000065874D0, : +0.000002424101735D0, : -0.004857669684959D0, : -0.000027184471371D0, : +1.000009560363559D0 / * Pick up B1950 data (units radians and arcsec/TC) R=R1950 D=D1950 UR=DR1950*PMF UD=DD1950*PMF PX=P1950 RV=V1950 * Spherical to Cartesian SR=SIN(R) CR=COS(R) SD=SIN(D) CD=COS(D) R0(1)=CR*CD R0(2)=SR*CD R0(3)= SD W=VF*RV*PX RD0(1)=-SR*CD*UR-CR*SD*UD+W*R0(1) RD0(2)= CR*CD*UR-SR*SD*UD+W*R0(2) RD0(3)= CD*UD+W*R0(3) * Allow for e-terms and express as position+velocity 6-vector W=R0(1)*A(1)+R0(2)*A(2)+R0(3)*A(3) WD=R0(1)*AD(1)+R0(2)*AD(2)+R0(3)*AD(3) DO I=1,3 V1(I)=R0(I)-A(I)+W*R0(I) V1(I+3)=RD0(I)-AD(I)+WD*R0(I) END DO * Convert position+velocity vector to Fricke system DO I=1,6 W=0D0 DO J=1,6 W=W+EM(I,J)*V1(J) END DO V2(I)=W END DO * Revert to spherical coordinates X=V2(1) Y=V2(2) Z=V2(3) XD=V2(4) YD=V2(5) ZD=V2(6) RXYSQ=X*X+Y*Y RXYZSQ=RXYSQ+Z*Z RXY=SQRT(RXYSQ) RXYZ=SQRT(RXYZSQ) SPXY=X*XD+Y*YD SPXYZ=SPXY+Z*ZD IF (X.EQ.0D0.AND.Y.EQ.0D0) THEN R=0D0 ELSE R=ATAN2(Y,X) IF (R.LT.0.0D0) R=R+D2PI END IF D=ATAN2(Z,RXY) IF (RXY.GT.TINY) THEN UR=(X*YD-Y*XD)/RXYSQ UD=(ZD*RXYSQ-Z*SPXY)/(RXYZSQ*RXY) END IF IF (PX.GT.TINY) THEN RV=SPXYZ/(PX*RXYZ*VF) PX=PX/RXYZ END IF * Return results R2000=R D2000=D DR2000=UR/PMF DD2000=UD/PMF V2000=RV P2000=PX END SUBROUTINE sla_EG50 (DR, DD, DL, DB) *+ * - - - - - * E G 5 0 * - - - - - * * Transformation from B1950.0 'FK4' equatorial coordinates to * IAU 1958 galactic coordinates (double precision) * * Given: * DR,DD dp B1950.0 'FK4' RA,Dec * * Returned: * DL,DB dp galactic longitude and latitude L2,B2 * * (all arguments are radians) * * Called: * sla_DCS2C, sla_DMXV, sla_DCC2S, sla_SUBET, sla_DRANRM, sla_DRANGE * * Note: * The equatorial coordinates are B1950.0 'FK4'. Use the * routine sla_EQGAL if conversion from J2000.0 coordinates * is required. * * Reference: * Blaauw et al, Mon.Not.R.Astron.Soc.,121,123 (1960) * * P.T.Wallace Starlink 10 May 1990 *- IMPLICIT NONE DOUBLE PRECISION DR,DD,DL,DB DOUBLE PRECISION sla_DRANRM,sla_DRANGE DOUBLE PRECISION V1(3),V2(3),R,D * * L2,B2 system of galactic coordinates * * P = 192.25 RA of galactic north pole (mean B1950.0) * Q = 62.6 inclination of galactic to mean B1950.0 equator * R = 33 longitude of ascending node * * P,Q,R are degrees * * * Equatorial to galactic rotation matrix * * The Euler angles are P, Q, 90-R, about the z then y then * z axes. * * +CP.CQ.SR-SP.CR +SP.CQ.SR+CP.CR -SQ.SR * * -CP.CQ.CR-SP.SR -SP.CQ.CR+CP.SR +SQ.CR * * +CP.SQ +SP.SQ +CQ * DOUBLE PRECISION RMAT(3,3) DATA RMAT(1,1),RMAT(1,2),RMAT(1,3), : RMAT(2,1),RMAT(2,2),RMAT(2,3), : RMAT(3,1),RMAT(3,2),RMAT(3,3)/ : -0.066988739415D0,-0.872755765852D0,-0.483538914632, : +0.492728466075D0,-0.450346958020D0,+0.744584633283, : -0.867600811151D0,-0.188374601723D0,+0.460199784784/ * Remove E-terms CALL sla_SUBET(DR,DD,1950D0,R,D) * Spherical to Cartesian CALL sla_DCS2C(R,D,V1) * Rotate to galactic CALL sla_DMXV(RMAT,V1,V2) * Cartesian to spherical CALL sla_DCC2S(V2,DL,DB) * Express angles in conventional ranges DL=sla_DRANRM(DL) DB=sla_DRANGE(DB) END SUBROUTINE sla_GE50 (DL, DB, DR, DD) *+ * - - - - - * G E 5 0 * - - - - - * * Transformation from IAU 1958 galactic coordinates to * B1950.0 'FK4' equatorial coordinates (double precision) * * Given: * DL,DB dp galactic longitude and latitude L2,B2 * * Returned: * DR,DD dp B1950.0 'FK4' RA,Dec * * (all arguments are radians) * * Called: * sla_DCS2C, sla_DIMXV, sla_DCC2S, sla_ADDET, sla_DRANRM, sla_DRANGE * * Note: * The equatorial coordinates are B1950.0 'FK4'. Use the * routine sla_GALEQ if conversion to J2000.0 coordinates * is required. * * Reference: * Blaauw et al, Mon.Not.R.Astron.Soc.,121,123 (1960) * * P.T.Wallace Starlink November 1988 *- IMPLICIT NONE DOUBLE PRECISION DL,DB,DR,DD DOUBLE PRECISION sla_DRANRM,sla_DRANGE DOUBLE PRECISION V1(3),V2(3),R,D,RE,DE * * L2,B2 system of galactic coordinates * * P = 192.25 RA of galactic north pole (mean B1950.0) * Q = 62.6 inclination of galactic to mean B1950.0 equator * R = 33 longitude of ascending node * * P,Q,R are degrees * * * Equatorial to galactic rotation matrix * * The Euler angles are P, Q, 90-R, about the z then y then * z axes. * * +CP.CQ.SR-SP.CR +SP.CQ.SR+CP.CR -SQ.SR * * -CP.CQ.CR-SP.SR -SP.CQ.CR+CP.SR +SQ.CR * * +CP.SQ +SP.SQ +CQ * DOUBLE PRECISION RMAT(3,3) DATA RMAT(1,1),RMAT(1,2),RMAT(1,3), : RMAT(2,1),RMAT(2,2),RMAT(2,3), : RMAT(3,1),RMAT(3,2),RMAT(3,3)/ : -0.066988739415D0,-0.872755765852D0,-0.483538914632, : +0.492728466075D0,-0.450346958020D0,+0.744584633283, : -0.867600811151D0,-0.188374601723D0,+0.460199784784/ * Spherical to Cartesian CALL sla_DCS2C(DL,DB,V1) * Rotate to mean B1950.0 CALL sla_DIMXV(RMAT,V1,V2) * Cartesian to spherical CALL sla_DCC2S(V2,R,D) * Introduce E-terms CALL sla_ADDET(R,D,1950D0,RE,DE) * Express in conventional ranges DR=sla_DRANRM(RE) DD=sla_DRANGE(DE) END DOUBLE PRECISION FUNCTION sla_DRANGE (ANGLE) *+ * - - - - - - - * D R A N G E * - - - - - - - * * Normalise angle into range +/- pi (double precision) * * Given: * ANGLE dp the angle in radians * * The result (double precision) is ANGLE expressed in the range +/- pi. * * P.T.Wallace Starlink 6 April 1990 *- IMPLICIT NONE DOUBLE PRECISION ANGLE DOUBLE PRECISION DPI,D2PI PARAMETER (DPI=3.141592653589793238462643D0) PARAMETER (D2PI=6.283185307179586476925287D0) sla_DRANGE=MOD(ANGLE,D2PI) IF (ABS(sla_DRANGE).GE.DPI) : sla_DRANGE=sla_DRANGE-SIGN(D2PI,ANGLE) END SUBROUTINE sla_DCC2S (V, A, B) *+ * - - - - - - * D C C 2 S * - - - - - - * * Direction cosines to spherical coordinates (double precision) * * Given: * V d(3) x,y,z vector * * Returned: * A,B d spherical coordinates in radians * * The spherical coordinates are longitude (+ve anticlockwise * looking from the +ve latitude pole) and latitude. The * Cartesian coordinates are right handed, with the x axis * at zero longitude and latitude, and the z axis at the * +ve latitude pole. * * If V is null, zero A and B are returned. * At either pole, zero A is returned. * * P.T.Wallace Starlink July 1989 *- IMPLICIT NONE DOUBLE PRECISION V(3),A,B DOUBLE PRECISION X,Y,Z,R X = V(1) Y = V(2) Z = V(3) R = SQRT(X*X+Y*Y) IF (R.EQ.0D0) THEN A = 0D0 ELSE A = ATAN2(Y,X) END IF IF (Z.EQ.0D0) THEN B = 0D0 ELSE B = ATAN2(Z,R) END IF END SUBROUTINE sla_DMXV (DM, VA, VB) *+ * - - - - - * D M X V * - - - - - * * Performs the 3-D forward unitary transformation: * * vector VB = matrix DM * vector VA * * (double precision) * * Given: * DM dp(3,3) matrix * VA dp(3) vector * * Returned: * VB dp(3) result vector * * P.T.Wallace Starlink March 1986 *- IMPLICIT NONE DOUBLE PRECISION DM(3,3),VA(3),VB(3) INTEGER I,J DOUBLE PRECISION W,VW(3) * Matrix DM * vector VA -> vector VW DO J=1,3 W=0D0 DO I=1,3 W=W+DM(J,I)*VA(I) END DO VW(J)=W END DO * Vector VW -> vector VB DO J=1,3 VB(J)=VW(J) END DO END DOUBLE PRECISION FUNCTION sla_DRANRM (ANGLE) *+ * - - - - - - - * D R A N R M * - - - - - - - * * Normalise angle into range 0-2 pi (double precision) * * Given: * ANGLE dp the angle in radians * * The result is ANGLE expressed in the range 0-2 pi (double * precision). * * P.T.Wallace Starlink December 1984 *- IMPLICIT NONE DOUBLE PRECISION ANGLE DOUBLE PRECISION D2PI PARAMETER (D2PI=6.283185307179586476925287D0) sla_DRANRM=MOD(ANGLE,D2PI) IF (sla_DRANRM.LT.0D0) sla_DRANRM=sla_DRANRM+D2PI END SUBROUTINE sla_ADDET (RM, DM, EQ, RC, DC) *+ * - - - - - - * A D D E T * - - - - - - * * Add the E-terms (elliptic component of annual aberration) * to a pre IAU 1976 mean place to conform to the old * catalogue convention (double precision) * * Given: * RM,DM dp RA,Dec (radians) without E-terms * EQ dp Besselian epoch of mean equator and equinox * * Returned: * RC,DC dp RA,Dec (radians) with E-terms included * * Called: * sla_ETRMS, sla_DCS2C, sla_DCC2S, sla_DRANRM, sla_DRANGE * * Explanation: * Most star positions from pre-1984 optical catalogues (or * derived from astrometry using such stars) embody the * E-terms. If it is necessary to convert a formal mean * place (for example a pulsar timing position) to one * consistent with such a star catalogue, then the RA,Dec * should be adjusted using this routine. * * Reference: * Explanatory Supplement to the Astronomical Ephemeris, * section 2D, page 48. * * P.T.Wallace Starlink July 1986 *- IMPLICIT NONE DOUBLE PRECISION RM,DM,EQ,RC,DC DOUBLE PRECISION sla_DRANRM DOUBLE PRECISION A(3),V(3) INTEGER I * E-terms vector CALL sla_ETRMS(EQ,A) * Spherical to Cartesian CALL sla_DCS2C(RM,DM,V) * Include the E-terms DO I=1,3 V(I)=V(I)+A(I) END DO * Cartesian to spherical CALL sla_DCC2S(V,RC,DC) * Bring RA into conventional range RC=sla_DRANRM(RC) END SUBROUTINE sla_DCS2C (A, B, V) *+ * - - - - - - * D C S 2 C * - - - - - - * * Spherical coordinates to direction cosines (double precision) * * Given: * A,B dp spherical coordinates in radians * (RA,Dec), (Long,Lat) etc * * Returned: * V dp(3) x,y,z unit vector * * The spherical coordinates are longitude (+ve anticlockwise * looking from the +ve latitude pole) and latitude. The * Cartesian coordinates are right handed, with the x axis * at zero longitude and latitude, and the z axis at the * +ve latitude pole. * * P.T.Wallace Starlink October 1984 *- IMPLICIT NONE DOUBLE PRECISION A,B,V(3) DOUBLE PRECISION COSB COSB=COS(B) V(1)=COS(A)*COSB V(2)=SIN(A)*COSB V(3)=SIN(B) END SUBROUTINE sla_SUBET (RC, DC, EQ, RM, DM) *+ * - - - - - - * S U B E T * - - - - - - * * Remove the E-terms (elliptic component of annual aberration) * from a pre IAU 1976 catalogue RA,Dec to give a mean place * (double precision) * * Given: * RC,DC dp RA,Dec (radians) with E-terms included * EQ dp Besselian epoch of mean equator and equinox * * Returned: * RM,DM dp RA,Dec (radians) without E-terms * * Called: * sla_ETRMS, sla_DCS2C, sla_,DVDV, sla_DCC2S, sla_DRANRM * * Explanation: * Most star positions from pre-1984 optical catalogues (or * derived from astrometry using such stars) embody the * E-terms. This routine converts such a position to a * formal mean place (allowing, for example, comparison with a * pulsar timing position). * * Reference: * Explanatory Supplement to the Astronomical Ephemeris, * section 2D, page 48. * * P.T.Wallace Starlink 10 May 1990 *- IMPLICIT NONE DOUBLE PRECISION RC,DC,EQ,RM,DM DOUBLE PRECISION sla_DRANRM,sla_DVDV DOUBLE PRECISION A(3),V(3),F INTEGER I * E-terms CALL sla_ETRMS(EQ,A) * Spherical to Cartesian CALL sla_DCS2C(RC,DC,V) * Include the E-terms F=1D0+sla_DVDV(V,A) DO I=1,3 V(I)=F*V(I)-A(I) END DO * Cartesian to spherical CALL sla_DCC2S(V,RM,DM) * Bring RA into conventional range RM=sla_DRANRM(RM) END SUBROUTINE sla_DIMXV (DM, VA, VB) *+ * - - - - - - * D I M X V * - - - - - - * * Performs the 3-D backward unitary transformation: * * vector VB = (inverse of matrix DM) * vector VA * * (double precision) * * (n.b. the matrix must be unitary, as this routine assumes that * the inverse and transpose are identical) * * Given: * DM dp(3,3) matrix * VA dp(3) vector * * Returned: * VB dp(3) result vector * * P.T.Wallace Starlink March 1986 *- IMPLICIT NONE DOUBLE PRECISION DM(3,3),VA(3),VB(3) INTEGER I,J DOUBLE PRECISION W,VW(3) * Inverse of matrix DM * vector VA -> vector VW DO J=1,3 W=0D0 DO I=1,3 W=W+DM(I,J)*VA(I) END DO VW(J)=W END DO * Vector VW -> vector VB DO J=1,3 VB(J)=VW(J) END DO END DOUBLE PRECISION FUNCTION sla_DVDV (VA, VB) *+ * - - - - - * D V D V * - - - - - * * Scalar product of two 3-vectors (double precision) * * Given: * VA dp(3) first vector * VB dp(3) second vector * * The result is the scalar product VA.VB (double precision) * * P.T.Wallace Starlink November 1984 *- IMPLICIT NONE DOUBLE PRECISION VA(3),VB(3) sla_DVDV=VA(1)*VB(1)+VA(2)*VB(2)+VA(3)*VB(3) END SUBROUTINE sla_ETRMS (EP, EV) *+ * - - - - - - * E T R M S * - - - - - - * * Compute the E-terms (elliptic component of annual aberration) * vector (double precision) * * Given: * EP dp Besselian epoch * * Returned: * EV dp(3) E-terms as (dx,dy,dz) * * Note the use of the J2000 aberration constant (20.49552 arcsec). * This is a reflection of the fact that the E-terms embodied in * existing star catalogues were computed from a variety of * aberration constants. Rather than adopting one of the old * constants the latest value is used here. * * References: * 1 Smith, C.A. et al., 1989. Astr.J. 97, 265. * 2 Yallop, B.D. et al., 1989. Astr.J. 97, 274. * * P.T.Wallace Starlink 10 April 1990 *- IMPLICIT NONE DOUBLE PRECISION EP,EV(3) * Arcseconds to radians DOUBLE PRECISION AS2R PARAMETER (AS2R=0.4848136811095359949D-5) DOUBLE PRECISION T,E,E0,P,EK,CP * Julian centuries since B1950 T=(EP-1950D0)*1.00002135903D-2 * Eccentricity E=0.01673011D0-(0.00004193D0+0.000000126D0*T)*T * Mean obliquity E0=(84404.836D0-(46.8495D0+(0.00319D0+0.00181D0*T)*T)*T)*AS2R * Mean longitude of perihelion P=(1015489.951D0+(6190.67D0+(1.65D0+0.012D0*T)*T)*T)*AS2R * E-terms EK=E*20.49552D0*AS2R CP=COS(P) EV(1)= EK*SIN(P) EV(2)=-EK*CP*COS(E0) EV(3)=-EK*CP*SIN(E0) END *DECK CPOSTN C C C CHARACTER*(*) FUNCTION CPOSTN(POS1,NDEC,ISIGN) C C FORM A CHARACTER STRING FOR A POSITION POS1 IN SECONDS WITH NDEC C DECIMAL PLACES. IF ISIGN.EQ.0 SUPRESS + SIGN. C CHARACTER FMT*13 DOUBLE PRECISION POS,POS1,DROUND,X C C ROUND THE POSITION TO THE REQUESTED NUMBER OF DECIMAL PLACES. C ND=MIN(27,MAX(0,NDEC)) POS=DROUND(POS1,ND) C C DEAL WITH SIGN. C CPOSTN=' ' IF(POS.LT.0.0) THEN CPOSTN(1:1)='-' IC=1 ELSEIF(ISIGN.NE.0) THEN CPOSTN(1:1)='+' IC=1 ELSE IC=0 ENDIF C C DELIMITERS. C CPOSTN(3+IC:3+IC)=':' CPOSTN(6+IC:6+IC)=':' C C HOURS OR DEGREES. C X=ABS(POS) I=X/3600.0 WRITE(CPOSTN(1+IC:2+IC),100) I C C MINUTES. C X=X-3600.0*I I=X/60.0 WRITE(CPOSTN(4+IC:5+IC),100) I C C SECONDS WITHOUT FRACTIONAL PART. C X=X-60.0*I IF(ND.LE.0) THEN I=NINT(X) WRITE(CPOSTN(7+IC:8+IC),100) I C C SECONDS WITH FRACTIONAL PART. C ELSE WRITE(FMT,101) ND+3,NDEC WRITE(CPOSTN(7+IC:ND+9+IC),FMT) X C C REPLACE SPACES WITH ZEROS. C DO 1 I=7+IC,8+IC IF(CPOSTN(I:I).EQ.' ') CPOSTN(I:I)='0' 1 CONTINUE ENDIF RETURN C C FORMATS. C 100 FORMAT(SS,I2.2) 101 FORMAT(SS,'(SS,F',I3,'.',I3,')') C C END OF CHARACTER*(*) FUNCTION CPOSTN. C END *DECK PARSE C C ********************************************************************* SUBROUTINE PARSE ( STRING, PARS, NPAR, MAXPAR, DELIM ) C ********************************************************************* C C SPLITS THE CHARACTER ARRAY STRING INTO PARAMETERS SEPARATED BY C THE DELIMITERS DELIM C THE NTH CHARACTER OF DELIM IS USED TO TERMINATE THE NTH C PARAMETER, THE FINAL CHARACTER BEING REPEATED IF NECESSARY C THE PARAMETERS ARE RETURNED IN PARS(MAXPAR) AND THEIR NUMBER IS NPAR C CHARACTER*(*) STRING,PARS(MAXPAR),DELIM C C L HOLDS THE POSITION OF THE LAST DELIMITER FOUND C IDEL POINTS TO THE CURRENT DELIMITER CHARACTER C LENDEL = LEN(DELIM) IDEL = 1 NPAR=0 L=0 DO 10 I=1,LEN(STRING) IF ( STRING(I:I).EQ.DELIM(IDEL:IDEL) ) THEN IF ( I-L.GT.1 ) THEN C C NON-ZERO LENGTH STRING BETWEEN DELIMITERS FOUND C CHECK NUMBER OF PARAMETERS AGAINST MAXPAR C STORE CURRENT PARAMETER C UPDATE DELIMITER CHARACTER C IF ( NPAR+1.GT.MAXPAR ) GOTO 100 NPAR = NPAR+1 PARS(NPAR) = STRING(L+1:I-1) IDEL = MIN(IDEL+1,LENDEL) ENDIF L=I ENDIF 10 CONTINUE C C IF A STRING IS REMAINING AT THE END, ADD IT TO PARS C IF ( L.LT.LEN(STRING) ) THEN IF ( NPAR+1.GT.MAXPAR ) GOTO 100 NPAR = NPAR+1 PARS(NPAR) = STRING(L+1:LEN(STRING)) ENDIF C 100 CONTINUE RETURN C C END OF SUBROUTINE PARSE C END *DECK DROUND C C C DOUBLE PRECISION FUNCTION DROUND(DVAL,NDEC) C C THIS ROUTINE RETURNS THE VALUE OF DVAL ROUNDED TO NDEC DECIMAL PLACES C IN THE SAME WAY AS THE FORTRAN OUTPUT ROUTINES. DVAL AND THE C FUNCTION VALUE ARE BOTH DOUBLE PRECISION. C CHARACTER BUFFER*30,FMT*8 DOUBLE PRECISION DVAL,D SAVE FMT DATA FMT/'(F30.??)'/ C C ADJUST THE FORMAT TO THE REQUESTED NUMBER OF DECIMAL PLACES. C ND=MIN(27,MAX(0,NDEC)) WRITE(FMT(6:7),800) ND C C WRITE DVAL INTO BUFFER WITH THIS NUMBER OF DECIMAL PLACES. C WRITE(BUFFER,FMT,ERR=100,IOSTAT=I) DVAL C C WRITE SUCCESSFUL, SET DROUND TO THIS VALUE. C READ(BUFFER,FMT,ERR=100,IOSTAT=I) D DROUND=D RETURN C C WRITE OR READ FAILED, RETURN AN ARITHMETICALLY ROUNDED VALUE. C 100 CONTINUE D=NINT(DVAL*10**ND)*10**(-ND) RETURN C C FORMAT STATEMENT. C 800 FORMAT(I2) C C END OF DOUBLE PRECISION FUNCTION DROUND. C END