* Last processed by NICE on 12-Jun-2000 15:53:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
C nic_astro.f
      SUBROUTINE PRECMAT(EPOCHFR,EPOCHTO,AMAT)
      INCLUDE 'const.inc'
      REAL*8    AMAT(9)
      REAL*8    EPOCHFR,EPOCHTO
      REAL*8    T,T2,T3
      REAL*8    ZA,ZT,TH
      REAL*8    S1,S2,S3,C1,C2,C3
      REAL*8    A,B,C,D
      T = (EPOCHTO-EPOCHFR)/100D0
      T2 = T*T
      T3 = T*T2
      IF (ABS(EPOCHFR-1950D0).LT.ABS(EPOCHFR-2000D0)) THEN
*
* J1950 corrections from Almanach 1981  (p B18)
*
         ZA = (0.6402633D0*T + 0.0000839D0*T2 + 0.0000050D0*T3)
     $   *DEG_TO_RAD
         ZT = (0.6402633D0*T + 0.0003036D0*T2 + 0.0000050D0*T3)
     $   *DEG_TO_RAD
         TH = (0.5567376D0*T - 0.0001183D0*T2 - 0.0000117D0*T3)
     $   *DEG_TO_RAD
      ELSE
*
* J2000 corrections from Almanach 1995  (p B18)
*
         ZA = (0.6406161D0*T + 0.0000839D0*T2 + 0.0000050D0*T3)
     $   *DEG_TO_RAD
         ZT = (0.6406161D0*T + 0.0003041D0*T2 + 0.0000051D0*T3)
     $   *DEG_TO_RAD
         TH = (0.5567530D0*T - 0.0001185D0*T2 - 0.0000116D0*T3)
     $   *DEG_TO_RAD
      ENDIF
      S1 = DSIN(ZA)
      C1 = DCOS(ZA)
      S2 = DSIN(ZT)
      C2 = DCOS(ZT)
      S3 = DSIN(TH)
      C3 = DCOS(TH)
      A = S1 * S2
      B = C1 * S2
      C = S1 * C2
      D = C1 * C2
      AMAT(1) = D * C3 - A
      AMAT(2) = -C * C3 - B
      AMAT(3) = -S3 * C2
      AMAT(4) = B * C3 + C
      AMAT(5) = -A * C3 + D
      AMAT(6) = -S3 * S2
      AMAT(7) = C1 * S3
      AMAT(8) = -S1 * S3
      AMAT(9) = C3
      RETURN
      END
*
      SUBROUTINE AMSET(OBSLAT,AMAT)
c**************************************************************************
c   subroutine to provide transformation matrix from hour angle - dec     *
c       into azimuth - elevation.                                         *
c**************************************************************************
*
      INCLUDE 'const.inc'
*
      REAL*8  AMAT(9)
      REAL*8  OBSLAT
      REAL*8  X1,Y1,Z1
*
      X1=DEG_TO_RAD*OBSLAT
      Y1=DSIN(X1)
      Z1=DCOS(X1)
      AMAT(1)=-Y1
      AMAT(2)=0.D0
      AMAT(3)=Z1
      AMAT(4)=0.D0
      AMAT(5)=-1.D0
      AMAT(6)=0.D0
      AMAT(7)=Z1
      AMAT(8)=0.D0
      AMAT(9)=Y1
      RETURN
      END
*
*
      SUBROUTINE MOVE1(NYRI,NYRF,MO,NDA,RA,D,DELR,DELD,DC)
c----------------------------------------------------------------------
c this is a simple precession routine.
c    nyri = initial epoch in yrs a.d. (i.e. 1950).
c    nyrf = the year to which we are precessing.
c    mo = month in year nyrf (1 to 12)
c    nda = day in month mo.
c         if day no. is known put in nda with mo set to 1.
c    ra = right ascension of source in radians.
c    d = declination of source in radians.
c         best results if ra & d are for the mean epoch.
c
c    delr = corrn to be added to mean r.a. for year nyri.
c    deld = corrn to be added to mean dec for year nyri.
c 10 feb 88. add duffett-smith nutation algorithm. should do to 1".
c----------------------------------------------------------------------
*
      REAL*8  RA,D,DELR,DELD,DC
      INTEGER NYRI,NYRF,MO,NDA
*
      REAL*8  SND,CSD,TND
      REAL*8  CSR,SNR
      REAL*8  AL,TO,T
      REAL*8  ZETAO,Z
      REAL*8  THETA
      REAL*8  AM,AN,ALAM
      REAL*8  SNL,CSL
      REAL*8  DLONG,DOBLQ
*
      REAL*8  LS,LD,MS,MD,NM
*
      REAL*8  CONST,T1900,A,B,TLS,TNM
      REAL*8  TLD,DPSI,DEPS
*
      SND = DSIN(D)
      CSD = DCOS(D)
      TND = SND/CSD
      CSR = DCOS(RA)
      SNR = DSIN(RA)
      AL = 30.5D0 * (MO - 1) + NDA
      TO = DBLE(NYRI - 1900)/100.0D0
      T = (DBLE(NYRF - NYRI) + AL / 365.2421988D0) / 100.0D0
      ZETAO=(2304.25D0+1.396D0*TO)*T+0.302D0*T**2
     $+0.018D0*T**3
      Z=ZETAO+0.791D0*T**2
      THETA=(2004.682D0-0.853D0*TO)*T-0.426D0*T**2
     $-0.042D0*T**3
      AM=(ZETAO+Z)*4.848136811D-6
      AN=THETA*4.848136811D-6
      ALAM=(0.985647D0*AL+278.5D0)*0.0174532925D0
      SNL=DSIN(ALAM)
      CSL=DCOS(ALAM)
      DELR=-9.92413605D-5*(SNL*SNR+0.91745051D0*CSL*CSR)/CSD
     $+AM+AN*SNR*TND
      DELD=-9.92413605D-5*(SNL*CSR*SND-
     $0.91745051D0*CSL*SNR*SND
     $+0.39784993D0*CSL*CSD)+AN*CSR
c	omega=259.183275d0-1934.142d0*(to+t)
c	arg=omega*0.0174532925d0
c	dlong=-8.3597d-5*dsin(arg)
c	doblq=4.4678d-5*dcos(arg)
c      delr=delr+dlong*(0.91745051d0+0.39784993d0*snr*tnd)-csr*tnd*doblq
c      deld=deld+0.39784993d0*csr*dlong+snr*doblq
      CONST = 1.745329252D-2
      T1900=TO+T
c      print *,'t1900 = ',t1900
      A=1.000021358D2*T1900
      B=360.D0*(A-INT(A))
      LS=2.796967D2+3.03D-4*T1900*T1900+B
      A=1.336855231D3*T1900
      B=360.D0*(A-INT(A))
      LD=2.704342D2-1.133D-3*T1900*T1900+B
      A=9.999736056D1*T1900
      B=360.D0*(A-INT(A))
      MS=3.584758D2-1.5D-4*T1900*T1900+B
      A=1.325552359D3*T1900
      B=360.D0*(A-INT(A))
      MD=2.961046D2+9.192D-3*T1900*T1900+B
      A=5.372616667D0*T1900
      B=360.D0*(A-INT(A))
      NM=2.591833D2+2.078D-3*T1900*T1900-B
c      print *,' nm = ',nm
      TLS=2*CONST*LS
      NM=CONST*NM
      TNM=2*NM
      MS=CONST*MS
      TLD=2*CONST*LD
      MD=CONST*MD
      DPSI=(-17.2327D0-1.737D-2*T1900)*SIN(NM)
     $+(-1.2729D0-1.3D-4*T1900)*SIN(TLS)
     $+2.088D-1*SIN(TNM)-2.037D-1*SIN(TLD)
     $+(1.261D-1-3.1D-4*T1900)*SIN(MS)
     $+6.75D-2*SIN(MD)
     $-(4.97D-2-1.2D-4*T1900)*SIN(TLS+MS)
     $-3.42D-2*SIN(TLD-NM)-2.61D-2*SIN(TLD+MD)
     $+2.14D-2*SIN(TLS-MS)
     $-1.49D-2*SIN(TLS-TLD+MD)
     $+1.24D-2*SIN(TLS-NM)+1.14D-2*SIN(TLD-MD)
      DEPS=(9.21D0+9.1D-4*T1900)*COS(NM)
     $+(5.522D-1-2.9D-4*T1900)*COS(TLS)
     $-9.04D-2*COS(TNM)+8.84D-2*COS(TLD)
     $+2.16D-2*COS(TLS+MS)+1.83D-2*COS(TLD-NM)
     $+1.13D-2*COS(TLD+MD)-9.3D-3*COS(TLS-MS)
     $-6.6D-3*COS(TLS-NM)
c      print *,' dlong = ',dpsi,'doblq = ',deps
      DLONG=DPSI*4.848136812D-6
      DOBLQ=DEPS*4.848136812D-6
      DELR=DELR+DLONG*(0.91745051D0+0.39784993D0*SNR*TND)
     $-CSR*TND*DOBLQ
      DELD=DELD+0.39784993D0*CSR*DLONG+SNR*DOBLQ
      RETURN
      END
