c
      SUBROUTINE DR2AF(NDP,ANGLE,SIGN,IDMSF)
*   
*   
*  CONVERT ANGLE FOR OUTPUT 
*  ------------------------ 
*   
*   
*  GIVEN:   
*     NDP           NUMBER OF DECIMAL PLACES OF ARC SECONDS 
*     ANGLE         ANGLE IN RADIANS (DOUBLE PRECISION) 
*   
*  RETURNED:
*     SIGN          '+' OR '-'  
*     IDMSF         DEGREES, ARC MINUTES, ARC SECONDS, FRACTION 
*   
*   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
*   
      INTEGER NDP   
      DOUBLE PRECISION ANGLE
      CHARACTER*1 SIGN  
      INTEGER IDMSF(4)  
    
*  RADIANS TO ARC SECONDS   
      DOUBLE PRECISION R2AS 
      PARAMETER (R2AS=206264.80624709636D0) 
    
*  HANDLE SIGN  
      SIGN='+'  
      IF(ANGLE.LT.0D0) SIGN='-' 
    
*  FIELD UNITS IN TERMS OF LEAST SIGNIFICANT FIGURE 
      RAS=REAL(10**MAX(NDP,0))  
      RAM=RAS*60D0  
      RD=RAM*60D0   
    
*  ROUND ANGLE AND EXPRESS IN SMALLEST UNITS REQUIRED   
      A=NINT(RAS*R2AS*ABS(ANGLE))   
    
*  SEPARATE INTO FIELDS 
      AD=INT(A/RD)  
      A=A-AD*RD 
      AM=INT(A/RAM) 
      A=A-AM*RAM
      AS=INT(A/RAS) 
      AF=A-AS*RAS   
    
*  RETURN RESULTS   
      IDMSF(1)=INT(AD)  
      IDMSF(2)=INT(AM)  
      IDMSF(3)=INT(AS)  
      IDMSF(4)=INT(AF)  
    
      END   
      SUBROUTINE DR2TF(NDP,ANGLE,SIGN,IHMSF)
*   
*   
*  CONVERT ANGLE FOR OUTPUT, EXPRESSED AS TIME  
*  -------------------------------------------  
*   
*   
*  GIVEN:   
*     NDP           NUMBER OF DECIMAL PLACES OF SECONDS 
*     ANGLE         ANGLE IN RADIANS (DOUBLE PRECISION) 
*   
*  RETURNED:
*     SIGN          '+' OR '-'  
*     IHMSF         HOURS, MINUTES, SECONDS, FRACTION   
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      INTEGER NDP   
      DOUBLE PRECISION ANGLE
      CHARACTER*1 SIGN  
      INTEGER IHMSF(4)  
    
*  RADIANS TO SECONDS   
      DOUBLE PRECISION R2S  
      PARAMETER (R2S=13750.987083139757D0)  
    
*  HANDLE SIGN  
      SIGN='+'  
      IF(ANGLE.LT.0D0) SIGN='-' 
    
*  FIELD UNITS IN TERMS OF LEAST SIGNIFICANT FIGURE 
      RS=REAL(10**MAX(NDP,0))   
      RM=RS*60D0
      RH=RM*60D0
    
*  ROUND ANGLE AND EXPRESS IN SMALLEST UNITS REQUIRED   
      A=NINT(RS*R2S*ABS(ANGLE)) 
    
*  SEPARATE INTO FIELDS 
      AH=INT(A/RH)  
      A=A-AH*RH 
      AM=INT(A/RM)  
      A=A-AM*RM 
      AS=INT(A/RS)  
      AF=A-AS*RS
    
*  RETURN RESULTS   
      IHMSF(1)=INT(AH)  
      IHMSF(2)=INT(AM)  
      IHMSF(3)=INT(AS)  
      IHMSF(4)=INT(AF)  
    
      END   
      DOUBLE PRECISION FUNCTION DRANGE(ANGLE)   
*   
*   
*  NORMALISE ANGLE INTO RANGE +/- PI
*  ---------------------------------
*   
*   
*  ANGLE IS THE ANGLE IN RADIANS
*  BOTH ANGLE AND THE RESULT ARE DOUBLE PRECISION   
*   
    
      DOUBLE PRECISION ANGLE
    
      DOUBLE PRECISION DPI,D2PI 
      PARAMETER (DPI=3.141592653589793238462643D0)  
      PARAMETER (D2PI=6.283185307179586476925287D0) 
    
    
      DRANGE=DMOD(ANGLE,D2PI)   
      IF(DABS(DRANGE).GE.DPI) DRANGE=DRANGE-DSIGN(D2PI,ANGLE)   
    
      END   
      DOUBLE PRECISION FUNCTION DRANRM(ANGLE)   
*   
*   
*  NORMALISE ANGLE INTO RANGE 0-2PI 
*  -------------------------------- 
*   
*   
*  ANGLE IS THE ANGLE IN RADIANS
*  BOTH ANGLE AND THE RESULT ARE DOUBLE PRECISION   
    
      DOUBLE PRECISION ANGLE
    
      DOUBLE PRECISION D2PI 
      PARAMETER (D2PI=6.283185307179586476925287D0) 
    
    
      DRANRM=DMOD(ANGLE,D2PI)   
      IF(DRANRM.LT.0D0) DRANRM=DRANRM+D2PI  
    
      END   
      SUBROUTINE DS2TP(RA,DEC,RAZ,DECZ,XI,ETA)  
*   
*   
*  PROJECTION OF SPHERICAL COORDINATES ONTO TANGENT PLANE   
*  ('GNOMONIC' PROJECTION - 'STANDARD COORDINATES') 
*   
*   
*  GIVEN:   
*     RA,DEC      SPHERICAL COORDINATES OF POINT TO BE PROJECTED
*     RAZ,DECZ    SPHERICAL COORDINATES OF TANGENT POINT
*   
*  RETURNED:
*     XI,ETA      RECTANGULAR COORDINATES ON TANGENT PLANE  
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION RA,DEC,RAZ,DECZ,XI,ETA   
    
    
      SDECZ=SIN(DECZ)   
      SDEC=SIN(DEC) 
      CDECZ=COS(DECZ)   
      CDEC=COS(DEC) 
      RADIF=RA-RAZ  
      SRADIF=SIN(RADIF) 
      CRADIF=COS(RADIF) 
    
      DENOM=SDEC*SDECZ+CDEC*CDECZ*CRADIF
    
      XI=CDEC*SRADIF/DENOM  
      ETA=(SDEC*CDECZ-CDEC*SDECZ*CRADIF)/DENOM  
    
      END   
      SUBROUTINE DTP2S(XI,ETA,RAZ,DECZ,RA,DEC)  
*   
*   
*  TRANSFORM TANGENT PLANE COORDINATES INTO SPHERICAL   
*   
*   
*  GIVEN:   
*     XI,ETA      TANGENT PLANE RECTANGULAR COORDINATES 
*     RAZ,DECZ    SPHERICAL COORDINATES OF TANGENT POINT
*   
*  RETURNED:
*     RA,DEC      SPHERICAL COORDINATES 
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION XI,ETA,RAZ,DECZ,RA,DEC   
    
    
      SDECZ=SIN(DECZ)   
      CDECZ=COS(DECZ)   
    
      DENOM=CDECZ-ETA*SDECZ 
      RADIF=ATAN2(XI,DENOM) 
    
      RA=DRANRM(RADIF+RAZ)  
      DEC=ATAN2(SDECZ+ETA*CDECZ,SQRT(XI*XI+DENOM*DENOM))
    
      END   
      SUBROUTINE MMP(RM0, DM0, PR, PD, EPOCH0, EPOCH1, RM1, DM1)
*   
*  GIVEN:   
*     RM0, DM0           MEAN RA, DEC AT 1ST EPOCH  
*     PR, PD             PROPER MOTIONS IN RA, DEC AT SPECIFIED EPOCH   
*                          (RADIANS PER YEAR)   
*     EPOCH0             1ST EPOCH  
*     EPOCH1             2ND EPOCH  
*   
*  RETURNED:
*     RM1, DM1           MEAN RA, DEC AT 2ND EPOCH  
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION RM0, DM0, PR, PD, EPOCH0, EPOCH1, RM1, DM1   
    
      DOUBLE PRECISION RMATP(3,3), V(3) 
    
    
*  PROPER MOTIONS   
    
      TAU= EPOCH1 - EPOCH0  
      R= RM0 + PR*TAU   
      D= DM0 + PD*TAU   
    
    
*  PRECESSION MATRIX
    
      CALL PREC(EPOCH0, EPOCH1, RMATP)  
    
    
*  MEAN TO MEAN 
    
      CALL DCS2C(R, D, V)   
      CALL RMATMU(V, RMATP) 
      CALL DCC2S(R, DM1, V) 
      RM1= DRANRM(R)
    
      END   
      SUBROUTINE PREC(EPOCH0, EPOCH1, S)
*   
*  GIVEN:   
*     EPOCH0             STARTING EPOCH (YEARS AD)  
*     EPOCH1             ENDING EPOCH (YEARS AD)
*   
*  RETURNED:
*     S                  PRECESSION MATRIX  
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION EPOCH0, EPOCH1, S(3,3)   
    
    
      TSER(C0,C1,C2)=DRANGE(C0+(C1+C2*T)*T) 
    
    
*  TIME FROM 1900.0 TO ORIGINAL EPOCH (TROPICAL CENTURIES)  
    
      T0=(EPOCH0 - 1900D0)/100D0
    
    
*  TIME SINCE ORIGINAL EPOCH (TROPICAL CENTURIES)   
    
      T= (EPOCH1 - EPOCH0)/100D0
      TSQ=T*T   
    
    
*  PRECESSION PARAMETERS
    
      Z0=TSER(0., 6.768D-6*T0 + 1.11713192D-2, 1.464D-6)
      Z=Z0 + 3.835D-6*TSQ   
      TH=TSER(0D0, -4.135D-6*T0 + 9.7189726D-3, -2.065D-6)  
    
      SZ0=SIN(Z0)   
      CZ0=COS(Z0)   
    
      SZ=SIN(Z) 
      CZ=COS(Z) 
    
      STH=SIN(TH)   
      CTH=COS(TH)   
    
      CTHSZ=CTH*SZ  
      CTHCZ=CTH*CZ  
    
    
*  ROTATION MATRIX  
    
      S(1,1)= CZ0*CTHCZ - SZ0*SZ
      S(1,2)=-SZ0*CTHCZ - CZ0*SZ
      S(1,3)=-STH*CZ
      S(2,1)= CZ0*CTHSZ + SZ0*CZ
      S(2,2)=-SZ0*CTHSZ + CZ0*CZ
      S(2,3)=-STH*SZ
      S(3,1)= CZ0*STH   
      S(3,2)=-SZ0*STH   
      S(3,3)= CTH   
    
      END   
      SUBROUTINE PMC(PR, PD, EPOCH, DATE, DR, DD)   
*   
*  GIVEN:   
*     PR, PD             PROPER MOTIONS IN RA, DEC AT SPECIFIED EPOCH   
*                          (RADIANS PER YEAR)   
*     EPOCH              YEAR AD OF MEAN EQUATOR & EQUINOX  
*     DATE               DAYS SINCE JD 2400000.0 FOR APPARENT PLACE 
*   
*  RETURNED:
*     DR, DD             TOTAL PROPER MOTION CORRECTIONS (RADIANS)  
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION PR, PD, EPOCH, DATE, DR, DD  
    
    
*  TIME SINCE ORIGINAL EPOCH (TROPICAL YEARS)   
    
      TAU= (DATE - 15020.313D0)/365.2421957D0 - (EPOCH - 1900D0)
    
    
*  PROPER MOTION CORRECTIONS (RA,DEC VECTOR)
    
      DR= PR*TAU
      DD= PD*TAU
    
      END   
      SUBROUTINE DCS2C(A,B,V)   
*   
*     A IS RA/HA/LONG/AZ
*     B IS DEC/LAT/ZD   
*     V IS X,Y,Z
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION A, B, V(3)   
    
    
      AR= DRANGE(A) 
      BR= DRANGE(B) 
    
      COSBR= COS(BR)
    
      V(1)= COS(AR)*COSBR   
      V(2)= SIN(AR)*COSBR   
      V(3)= SIN(BR) 
    
      END   
      SUBROUTINE DCC2S(A,B,V)   
*   
*     A IS RA/HA/LONG/AZ
*     B IS DEC/LAT/ZD   
*     V IS X,Y,Z
*   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION A, B, V(3)   
    
    
      X=V(1)
      Y=V(2)
    
      A=ATAN2(Y,X)  
      B=ATAN2(V(3), SQRT(X*X + Y*Y))
    
      END   
      SUBROUTINE RMATMU(A,B)
*   
*  PERFORMS THE 3-D FORWARD UNITARY TRANSFORMATION: 
*     VECTOR A = MATRIX B * VECTOR A
*   
    
      DOUBLE PRECISION A(3), B(3,3) 
    
      DOUBLE PRECISION C(3) 
    
    
      DO J=1,3  
         C(J)=0D0   
         DO I=1,3   
            C(J)=C(J)+A(I)*B(J,I)   
         END DO 
      END DO
    
      DO J=1,3  
         A(J)=C(J)  
      END DO
    
      END   
      FUNCTION RJDATE ( NYEAR, NDIY )   
*   
*   
*  JULIAN DATE FROM YEAR & DAY-IN-YEAR  
*   
*   
*  ARGUMENTS:   
*   
*        NYEAR GREGORIAN YEAR AD (INTEGER; MUST BE 1901 THROUGH 2100)   
*        NDIY  DAY IN YEAR NUMBER (INTEGER; JAN 1ST = DAY 1)
*   
*        RESULT IS THE JULIAN DATE MINUS 2400000.0 FOR 0 HRS UT 
*   
    
      RJDATE=REAL(NDIY+(NYEAR-1)/4)+REAL(NYEAR-1901)*365.0+14909.5  
    
      END   
    
