> Scott Daniel of the LSST project noticed two oddities in the Fortran
> SLALIB. One is harmless and the other is a bug. The latter was
> fixed long ago but I suspect is still in the Starlink version.
> Attached are revised GPL'ed copies of the two affected subroutines
> to replace the current Starlink ones.
>
> The harmless oddity was in sla_MAPQK. In the aberration code, a
> redundant normalization was applied, that had no effect on the
> spherical coordinates, just a minuscule efficiency hit.
>
> The bug was in sla_MAPQKZ, which neglected light deflection. This
> would have affected the returned position by typically a few mas
> for nighttime stars. And as I said earlier, it was caught years
> ago.
Hello,
Thank you for sending the updated files! I tried to apply them in the
Starlink libraries/sla/ directory, but I'm confused as the differences
(see git diff below) don't appear to be quite as you described them.
Aside from changes to the comments and copyright lines, the changes seem
to be:
mapqk.f:
* parameter VF altered from the 6th d.p.
mapqkz.f:
* normalization removed from abberation calculation
(removing variable P1DVP1)
Maybe we already had the fix for the light deflection bug and
normalization was removed from sla_MAPQKZ instead of sla_MAPQK?
Best regards,
Graham
diff --git a/libraries/sla/mapqk.f b/libraries/sla/mapqk.f
index 0f30a78..7586adb 100644
--- a/libraries/sla/mapqk.f
+++ b/libraries/sla/mapqk.f
@@ -30,9 +30,9 @@
* (1) time interval for proper motion (Julian years)
* (2-4) barycentric position of the Earth (AU)
* (5-7) heliocentric direction of the Earth (unit vector)
-* (8) (grav rad Sun)*2/(Sun-Earth distance)
+* (8) (Schwarzschild radius of Sun)/(Sun-Earth distance)
* (9-11) barycentric Earth velocity in units of c
-* (12) sqrt(1-v**2) where v=modulus(ABV)
+* (12) sqrt(1-v^2) where v=modulus(ABV)
* (13-21) precession/nutation (3,3) matrix
*
* Returned:
@@ -64,9 +64,9 @@
* sla_DCC2S Cartesian to spherical
* sla_DRANRM normalize angle 0-2Pi
*
-* P.T.Wallace Starlink 15 January 2000
+* Last revision: 23 November 2016
*
-* Copyright (C) 2000 Rutherford Appleton Laboratory
+* Copyright P.T.Wallace. All rights reserved.
*
* License:
* This program is free software; you can redistribute it and/or modify
@@ -81,8 +81,8 @@
*
* You should have received a copy of the GNU General Public License
* along with this program (see SLA_CONDITIONS); if not, write to the
-* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
-* Boston, MA 02110-1301 USA
+* Free Software Foundation, Inc., 59 Temple Place, Suite 330,
+* Boston, MA 02111-1307 USA
*
*-
@@ -96,7 +96,7 @@
* Km/s to AU/year
DOUBLE PRECISION VF
- PARAMETER (VF=0.21094502D0)
+ PARAMETER (VF=0.210949528D0)
INTEGER I
@@ -108,7 +108,7 @@
-* Unpack scalar and vector parameters
+* Unpack scalar and vector parameters.
PMT = AMPRMS(1)
GR2E = AMPRMS(8)
AB1 = AMPRMS(12)
@@ -118,17 +118,17 @@
ABV(I) = AMPRMS(I+8)
END DO
-* Spherical to x,y,z
+* Spherical to x,y,z.
CALL sla_DCS2C(RM,DM,Q)
-* Space motion (radians per year)
+* Space motion (radians per year).
PXR = PX*AS2R
W = VF*RV*PXR
EM(1) = -PR*Q(2)-PD*COS(RM)*SIN(DM)+W*Q(1)
EM(2) = PR*Q(1)-PD*SIN(RM)*SIN(DM)+W*Q(2)
EM(3) = PD*COS(DM) +W*Q(3)
-* Geocentric direction of star (normalized)
+* Geocentric direction of star (normalized).
DO I=1,3
P(I) = Q(I)+PMT*EM(I)-PXR*EB(I)
END DO
@@ -142,17 +142,17 @@
P1(I) = PN(I)+W*(EHN(I)-PDE*PN(I))
END DO
-* Aberration (normalization omitted)
+* Aberration (normalization omitted).
P1DV = sla_DVDV(P1,ABV)
W = 1D0+P1DV/(AB1+1D0)
DO I=1,3
P2(I) = AB1*P1(I)+W*ABV(I)
END DO
-* Precession and nutation
+* Precession and nutation.
CALL sla_DMXV(AMPRMS(13),P2,P3)
-* Geocentric apparent RA,Dec
+* Geocentric apparent RA,Dec.
CALL sla_DCC2S(P3,RA,DA)
RA = sla_DRANRM(RA)
diff --git a/libraries/sla/mapqkz.f b/libraries/sla/mapqkz.f
index 2e26bc8..a308adf 100644
--- a/libraries/sla/mapqkz.f
+++ b/libraries/sla/mapqkz.f
@@ -27,9 +27,9 @@
*
* (1-4) not used
* (5-7) heliocentric direction of the Earth (unit vector)
-* (8) (grav rad Sun)*2/(Sun-Earth distance)
+* (8) (Schwarzschild radius of Sun)/(Sun-Earth distance)
* (9-11) ABV: barycentric Earth velocity in units of c
-* (12) sqrt(1-v**2) where v=modulus(ABV)
+* (12) sqrt(1-v^2) where v=modulus(ABV)
* (13-21) precession/nutation (3,3) matrix
*
* Returned:
@@ -42,8 +42,8 @@
*
* Notes:
*
-* 1) The vectors AMPRMS(2-4) and AMPRMS(5-7) are referred to the
-* mean equinox and equator of epoch EQ.
+* 1) The vector AMPRMS(5-7) is referred to the mean equinox and
+* equator of epoch EQ.
*
* 2) Strictly speaking, the routine is not valid for solar-system
* sources, though the error will usually be extremely small.
@@ -56,9 +56,9 @@
*
* Called: sla_DCS2C, sla_DVDV, sla_DMXV, sla_DCC2S, sla_DRANRM
*
-* P.T.Wallace Starlink 18 March 1999
+* Last revision: 23 November 2016
*
-* Copyright (C) 1999 Rutherford Appleton Laboratory
+* Copyright P.T.Wallace. All rights reserved.
*
* License:
* This program is free software; you can redistribute it and/or modify
@@ -73,8 +73,8 @@
*
* You should have received a copy of the GNU General Public License
* along with this program (see SLA_CONDITIONS); if not, write to the
-* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
-* Boston, MA 02110-1301 USA
+* Free Software Foundation, Inc., 59 Temple Place, Suite 330,
+* Boston, MA 02111-1307 USA
*
*-
@@ -86,14 +86,13 @@
DOUBLE PRECISION GR2E,AB1,EHN(3),ABV(3),
: P(3),PDE,PDEP1,W,P1(3),P1DV,
- : P1DVP1,P2(3),P3(3)
+ : P2(3),P3(3)
DOUBLE PRECISION sla_DVDV,sla_DRANRM
-
-* Unpack scalar and vector parameters
+* Unpack scalar and vector parameters.
GR2E = AMPRMS(8)
AB1 = AMPRMS(12)
DO I=1,3
@@ -101,10 +100,10 @@
ABV(I) = AMPRMS(I+8)
END DO
-* Spherical to x,y,z
+* Spherical to x,y,z.
CALL sla_DCS2C(RM,DM,P)
-* Light deflection
+* Light deflection (restrained within the Sun's disc).
PDE = sla_DVDV(P,EHN)
PDEP1 = PDE+1D0
W = GR2E/MAX(PDEP1,1D-5)
@@ -112,18 +111,17 @@
P1(I) = P(I)+W*(EHN(I)-PDE*P(I))
END DO
-* Aberration
+* Aberration (normalization omitted).
P1DV = sla_DVDV(P1,ABV)
- P1DVP1 = P1DV+1D0
W = 1D0+P1DV/(AB1+1D0)
DO I=1,3
- P2(I) = (AB1*P1(I)+W*ABV(I))/P1DVP1
+ P2(I) = AB1*P1(I)+W*ABV(I)
END DO
-* Precession and nutation
+* Precession and nutation.
CALL sla_DMXV(AMPRMS(13),P2,P3)
-* Geocentric apparent RA,Dec
+* Geocentric apparent RA,Dec.
CALL sla_DCC2S(P3,RA,DA)
RA = sla_DRANRM(RA)
----
Starlink User Support list
For list configuration, including subscribing to and unsubscribing from the list, see
https://www.jiscmail.ac.uk/cgi-bin/webadmin?A0=STARLINK
