subroutine  preces (t1,ra1,dec1,t2,ra2,dec2)
c
c     calcul de la precession
c     d'apres un sous-programme de p.lucke
c     utilise la formule exact de s.newcomb
c
c     on donne...
c
c     t1          annee des cordonnees ra1 et dec1
c     ra1         ascension droite a l'epoque t1, en radians
c     dec1        declinaison a l'epoque t2, en radians
c     t2          annee des cordonnees ra2 et dec2
c
c     on recoit...
c
c     ra2         ascension droite a l'epoque t2, en radians
c     dec2        declinaison a l'epoque t2, en radians
c
      double precision ra1,dec1,ra2,dec2,twopi
c
      data twopi/6.283185307/
c
      to=(amax1(t1,t2)-1900.)/100.
      tt=abs(t2-t1)/100.
      zt=tt*((0.0111713+0.67679e-5*to)+tt*(1.45e-6+tt*8.7e-8))
      zz=zt+3.83e-6*tt*tt
      th=tt*((0.009718963-0.41355e-5*to)-tt*(2.08e-6+tt*2.0e-7))
      if (t2.gt.t1) goto 10
      save=-zt
      zt=-zz
      zz=save
      th=-th
   10 tth2=tan(th/2.)
      sinth=sin(th)
      aa=ra1+zt
      csa=cos(aa)
      pp=sinth*(tan(dec1)+tth2*csa)
      da=atan(pp*sin(aa)/(1.-pp*csa))
      ra2=ra1+da+zt+zz
      if (ra2.ge.0.) goto 20
      ra2=ra2+twopi
      goto 30
   20 if (ra2.lt.twopi) goto 30
      ra2=ra2-twopi
   30 dec2=dec1+2.0*atan(tth2*cos(da/2.+aa)/cos(da/2.))
      return
      end