C --- Subroutine generated on Mon Aug 27 09:17:59 CEST 2007

        subroutine fcnnut(mjd,X,Y,dX,dY)
        
C-------------------------------------------------------------------------C
C                                                                         C
C       Empirical model for FCN                                           C
C                                                                         C
C-------------------------------------------------------------------------C
C                                                                         C
C       Input: mjd, modified julian date                                  C
C                                                                         C
C       Output: X, Y, CIP offsets (microas)                               C
C               dX, dY, uncertainties on dX, dY (microas)                 C
C                                                                         C
C-------------------------------------------------------------------------C
C                                                                         C
C       Subroutine interface: none                                        C
C                                                                         C
C-------------------------------------------------------------------------C
C                                                                         C
C       Author: S. Lambert (s.lambert@oma.be)                             C
C                                                                         C
C-------------------------------------------------------------------------C

        implicit none
        
C       Input variables

        double precision mjd

C       Output variables
        
        double precision X,Y    ! FCN contributions
        double precision dX,dY  ! uncertainties on X, Y
        
C       Internal variables

        integer i,j,N
        parameter(N=24)
        double precision pi,per,phi,mpe
        double precision axc,axs,ayc,ays
        double precision daxc,daxs,dayc,days,dt,t
        double precision TABLE(4,N)
        double precision DATE(N),XC(N),XS(N),YC(N),YS(N),SX(N),SY(N)

        pi=3.14159265358979323846D0
        
C       Mean prediction error

        mpe=0.1325D0                       ! microas per day

C       FCN parameters

        per=-430.21D0                      ! period in days
        phi=(2.D0*pi/per)*(mjd-51544.5D0)  ! phase in rad
                
C       Block data of amplitudes for X (microas)
        
      data ((TABLE(i,j),i=1,4),j=1,N) /
     . 45700.D0,    78.96D0,  -125.26D0,    18.92D0, ! 1984.0
     . 46066.D0,   -93.86D0,  -128.36D0,    12.31D0, ! 1985.0
     . 46431.D0,  -237.90D0,  -154.16D0,    14.25D0, ! 1986.0
     . 46796.D0,  -262.72D0,  -152.24D0,    13.65D0, ! 1987.0
     . 47161.D0,  -201.58D0,   -27.63D0,    14.20D0, ! 1988.0
     . 47527.D0,  -191.86D0,   -87.75D0,    11.94D0, ! 1989.0
     . 47892.D0,  -201.03D0,   -60.17D0,     9.18D0, ! 1990.0
     . 48257.D0,  -147.84D0,    26.27D0,     7.93D0, ! 1991.0
     . 48622.D0,  -118.84D0,    40.64D0,     6.51D0, ! 1992.0
     . 48988.D0,  -111.84D0,     9.53D0,     7.50D0, ! 1993.0
     . 49353.D0,  -120.47D0,    34.42D0,     6.99D0, ! 1994.0
     . 49718.D0,  -123.23D0,    39.25D0,     7.67D0, ! 1995.0
     . 50083.D0,  -143.29D0,    17.16D0,     7.40D0, ! 1996.0
     . 50449.D0,  -127.16D0,    49.15D0,     7.50D0, ! 1997.0
     . 50814.D0,   -82.01D0,    30.81D0,     7.63D0, ! 1998.0
     . 51179.D0,   -20.04D0,   -35.49D0,     6.33D0, ! 1999.0
     . 51544.D0,    29.10D0,  -110.93D0,     5.95D0, ! 2000.0
     . 51910.D0,    87.13D0,  -123.81D0,     7.36D0, ! 2001.0
     . 52275.D0,    80.14D0,   -97.29D0,     7.13D0, ! 2002.0
     . 52640.D0,   109.91D0,   -54.17D0,     4.65D0, ! 2003.0
     . 53005.D0,   131.50D0,   -10.21D0,     4.55D0, ! 2004.0
     . 53371.D0,   157.75D0,   -15.35D0,     3.23D0, ! 2005.0
     . 53736.D0,   160.36D0,   -14.18D0,     3.09D0, ! 2006.0
     . 54101.D0,   141.55D0,    35.95D0,     4.72D0/ ! 2007.0
          
C       Amplitudes extracted from the table

        do i=1,N
           DATE(i)=TABLE(1,i)
           XC(i)=TABLE(2,i)
           XS(i)=TABLE(3,i)
           SX(i)=TABLE(4,i)
           YC(i)=XS(i)
           YS(i)=-XC(i)
           SY(i)=SX(i)
        end do

C       Prediction of the amplitude at the input date

        if (mjd.le.DATE(1)) then
           axc=XC(1)
           axs=XS(1)
           ayc=YC(1)
           ays=YS(1)
        else if (mjd.ge.DATE(N)) then
           axc=XC(N)
           axs=XS(N)
           ayc=YC(N)
           ays=YS(N)
        else
           do i=1,N-1
              if (mjd.ge.DATE(i).and.mjd.lt.DATE(i+1)) then
                 t=mjd-DATE(i)
                 dt=DATE(i+1)-DATE(i)
                 daxc=XC(i+1)-XC(i)
                 daxs=XS(i+1)-XS(i)
                 dayc=YC(i+1)-YC(i)
                 days=YS(i+1)-YS(i)
                 axc=XC(i)+(daxc/dt)*t
                 axs=XS(i)+(daxs/dt)*t
                 ayc=YC(i)+(dayc/dt)*t
                 ays=YS(i)+(days/dt)*t
              end if
           end do
        end if
        
C       Computation of X and Y

        X=axc*dcos(phi)-axs*dsin(phi)  ! microas
        Y=ayc*dcos(phi)-ays*dsin(phi)  ! microas
        
C       Prediction of the uncertainty at the input date

        if (mjd.le.DATE(1)) then
           axc=SX(1)+mpe*(DATE(1)-mjd)
           axs=SX(1)+mpe*(DATE(1)-mjd)
           ayc=SY(1)+mpe*(DATE(1)-mjd)
           ays=SY(1)+mpe*(DATE(1)-mjd)
        else if (mjd.ge.DATE(N)) then
           axc=SX(N)+mpe*(mjd-DATE(N))
           axs=SX(N)+mpe*(mjd-DATE(N))
           ayc=SY(N)+mpe*(mjd-DATE(N))
           ays=SY(N)+mpe*(mjd-DATE(N))
        else
           do i=1,N-1
              if (mjd.ge.DATE(i).and.mjd.lt.DATE(i+1)) then
                 t=mjd-DATE(i)
                 dt=DATE(i+1)-DATE(i)
                 daxc=SX(i+1)-SX(i)
                 daxs=SX(i+1)-SX(i)
                 dayc=SY(i+1)-SY(i)
                 days=SY(i+1)-SY(i)
                 axc=dabs(SX(i)+(daxc/dt)*t)
                 axs=dabs(SX(i)+(daxs/dt)*t)
                 ayc=dabs(SY(i)+(dayc/dt)*t)
                 ays=dabs(SY(i)+(days/dt)*t)
              end if
           end do
        end if
        
C       Computation of the uncertainties

        dX=axc+axs  ! microas
        dY=ayc+ays  ! microas
        
        end
