subroutine fourier(t1,rl,P0,Npts,ian,DISP,A0,A,B,MORD)
c
c***************************************************************
c     SINGLE PERIOD FOURIER FIT TO THE HYDRODYNAMIC TIME SERIES
c***************************************************************
c
c
c---------------------------------------------------------------
c
c     M(bol, Sun) = 4.75 ---> Lang: "Astrophysical Data: ..."
c
c     BC = BC (Kurucz'92) + 0.120 
c                        ---> Procyon, Sun, B-W  M_b <--> M_v 
c                             (Liu & Janes 1990, ApJ 354, 273)
c
c---------------------------------------------------------------
c
C     FOURIER DECOMPOSITION HAS THE FOLLOWING FORM:
C
C     A0 + A1*SIN(2*PI*F1*(T-T0)+FI1) + ...
C
C
C     INPUT  :  TEMP.DAT --- T(I), X(I), I=1,2,...,N
C
C     OUTPUT = COEFF.DAT --- FOURIER COEFFICIENTS FOR PLOTTING
C            =  FDEC.DAT --- STAR NAME, DISPERSION, AMPLITUDES, E.T.C.
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      character *10 name

      character *4 e8
      DOUBLE PRECISION X(5000),T(5000),G(200,201),W(200),F(200)
      DOUBLE PRECISION A(200),B(200),AF(200), pha(10000)
C
      double precision ar(200),pr(200),acc(5000)
      double precision svel(5000),slum(5000),stef(5000),srad(5000)
      double precision t1(2000),vel1(2000),rl(2000),rphot(2000)
     c ,tphot(2000),flm(2000)
C
C     T0    = EPOCH
C
C     MORD  = ORDER OF THE FIT (NUMBER OF THE FREQUENCY COMPONENTS)
C
C     INTER = 1 -- INTERPOLATE KURUCZ TABLES
C           = 0 -- USE ANALYTIC FIT TO THE KURUCZ TABLES
C
c
c**open file to output only the coeff
c
      open(55, file = 'mycoeff.dat', status = 'unknown')

      T0=0.0D0
C     MORD=6
      INTER=1
C
C*************************
      p7   = 8.0d0*datan(1.0d0)
      grav = 6.67259d-8
      sunt = 5780.0d0
      sunm = 1.9891d33
      sunr = 6.9599d10
      sunl = 3.8268d33
      sunb = 4.75d0
      sunz = 0.020d0
      bcc  = 0.120d0
C*************************
C
      ian = 3
c     write(*,*) ' Analyse  V_r, V_mag, M_bol ( 1, 2, 3 ):'
c     read(*,*) ian
C
C     READ T(I),X(I),  I=1,2,...,N
C
c** open file
c
c1      OPEN(1,FILE= 'test.dat')
c      READ(1,*) smass,slu,strad,steff,xc,zc
c      READ(1,*) P0
c1     N=0
c
c** reads the name, period and the number of obs
c
c1     read(1, *) name, P0, Npts
c
 400  CONTINUE
      N=N+1
c      READ(1,*,END=401) t(n),svel(n),slum(n),stef(n),srad(n)
c
c** read the time, phase, and the velocity
c
c1     read(1, *,end=401) t1(n),d0,d1,rl(n),e,e1,e2,e3,e4,e5,e6
c1    c,e7,e8,e9,e10
c1     write(28, *) t1(n),d0,d1,rl(n),e
c     read(1, *, end=401) T(n), pha(n), svel(n)
c
c1     t1(n) = t1(n) - t1(1)
c1     t1(n) = t1(n)/P0
c1     t1(n) = t1(n) - int(t1(n))
       Do 54 n = 1, Npts
       T(n) = t1(n)
        x(n) = rl(n)
        write(75,*)t1(n),x(n)
  54   continue
c     if(ian.eq.3) x(n) = rl(n)/100.0
c     if(ian.eq.3) x(n) = -2.5d0*dlog10(rl(n)) + sunb
C
c1     GO TO 400
 401  CONTINUE
c1     CLOSE(1)
c1     N=N-1
       N = Npts
c
c
c1     if(inter.eq.0) write(*,*) 'Working !'
c1     if(ian.ne.2) go to 1488
c1     if(inter.eq.1) write(*,*) 'Working ! ... it takes a bit long'
 1488 continue
c
c1     if(ian.ne.2) go to 410
c
c*********************************************
c     Compute bolometric correction
c
c*********************************************
c     Initialize interpolation routine
c
c      if(inter.eq.1) call KURINT(cfeh,clogg,cteff,0,value)
cc
c      feh   = dlog10(zc/sunz)
c      fac   = 1.0d5/(3600.0d0*24.0d0)
c      gmin  = 1.0d30
c      gmax  = -gmin
c      delma = -1.0d30
c      ibad  = 0
c
c      fac   = 1.0d5/(3600.0d0*24.0d0)
c
c     Compute acceleration
c
c      i1=1
c      i2=2
c      acc(1) = fac*(svel(i2)-svel(i1))/(t(i2)-t(i1))
c      i1=n-1
c      i2=n
c      acc(n) = fac*(svel(i2)-svel(i1))/(t(i2)-t(i1))
c      n1=n-1
c      do 1515 i=2,n1
c      i1    = i-1
c      i2    = i+1
c      acc(i)= fac*(svel(i2)-svel(i1))/(t(i2)-t(i1))
c 1515 continue
c
c     Check deviating points
c
c      sfac=0.02d0
c      is=idint(dfloat(n)*sfac)
c      if(is.lt.1) go to 1519
c      ie=n-is
c      is1=is+1
c      s=0.0d0
c
c      do 1517 i=is1,ie
c      accm=0.0d0
c      do 1520 k=1,is
c      accm=accm+acc(i-k)+acc(i+k)
c 1520 continue
c      accm=accm/dfloat(2*is)
c      d=acc(i)-accm
c      s=s+d*d
c 1517 continue
c      s=dsqrt(s/dfloat(n))
c      s3=3.0d0*s
c
c      do 1521 i=is1,ie
c      accm=0.0d0
c      do 1522 k=1,is
c      accm=accm+acc(i-k)+acc(i+k)
c 1522 continue
c      accm=accm/dfloat(2*is)
c      d=dabs(acc(i)-accm)
c      if(d.lt.s3) go to 1521
c      acc(i)=accm 
c 1521 continue
c 1519 continue
c
c     Compute  g_eff
c
c      ic=0
c      icm=50
c
c      do 490 i=1,n
c
c      ic=ic+1
c      if(ic.lt.icm) go to 2222
c      ic=0
c      write(*,*) i,n
c 2222 continue
c
c      radi  = sunr*srad(i)
c      gef   = grav*smass*sunm/(radi*radi) + acc(i)  
c
c>>>>>>>>>>>>>>>>>>>>>>  gef < 10.0  guard
c
c      if(gef.gt.10.0d0) ggg=gef
c      if(gef.lt.10.0d0) ibad=ibad+1
c      if(gef.lt.10.0d0) gef=ggg
c
c      gef = dlog10(gef) 
c
c      gmin=dmin1(gmin,gef)
c      gmax=dmax1(gmax,gef)
c
c**********************************************************
c     Analytic formula for the bolometric correction
c
c      bc   = -5.8602 + 1.5393*dlog10(stef(i))
c      bc   = bc + 0.0672*feh + 0.0110*feh*feh - 0.0362*gef
c
c     Calculate interpolated bolometric correction
c      
c      if(inter.eq.1) call KURINT(feh,gef,stef(i),1,value)
c
c      delta=dabs(value-bc)
c
c      delma=dmax1(delma,delta)
c
c      if(inter.eq.1) bc = value
c
c      bc   = bc + bcc
c      x(i) = -2.5d0*dlog10(slum(i)) + sunb - bc
c
c 490  continue
c
c      write(*,500) n,ibad
c 500  format('Number of data:',i5,/,'Number of small g_eff :',i5)
c      if(inter.eq.1) write(*,*) delma
c
c      write(*,488) gmin
c      write(*,489) gmax
c 488  format(1x,'g_min =',f7.2)
c 489  format(1x,'g_max =',f7.2,/)
C
C------------------------------
C     FIT TIME SERIES
C------------------------------
C
 410  continue
C
      M=MORD
      W(1) = 1.d0
c      W(1)=1.0D0/P0
      DO 887 I=1,M
      W(I)=W(1)*DFLOAT(I)
 887  CONTINUE
C
C     PUNCH PARAMETERS
C
c**** dummy parameters
c
      model=1
      index=1
      zz=1.0d0
      pfrac=1.0d0
c
      OPEN(7,FILE='coeff.dat')
      WRITE(7,*) MODEL
      WRITE(7,*) INDEX
      WRITE(7,*) ZZ
      WRITE(7,*) PFRAC
      WRITE(7,*) M
      DO 702 J=1,M
      WRITE(7,*) W(J)
 702  CONTINUE
C
       DO 30 I=1,N
       T(I)=T(I)-T0
 30    CONTINUE
       DO 300 J=1,M
       F(J)=W(J)
 300   CONTINUE
C
       M2=2*M
       M3=M2+1
       M4=M3+1
       CALL NORM(N,M,T,X,F,G)
       CALL MATINV(M,G)
       A0=0.0D+00
       DO 31 J=1,M3
  31   A0=A0+G(1,J)*G(J,M4)
       DO 32 J=2,M3
       S=0.0D+00
       DO 33 K=1,M3
  33   S=S+G(J,K)*G(K,M4)
       A(J)=S
       AF(J)=S
  32   CONTINUE
C
C************************************************
C     CALCULATION OF THE AMPLITUDES AND PHASES
C
       K=0
       DO 37 J=2,M3,2
       K=K+1
       J1=J+1
       AMP=DSQRT(A(J)*A(J)+A(J1)*A(J1))
       PH=A(J1)/AMP
       PH=DACOS(PH)
       IF(A(J).LT.0.0D+00) PH=P7-PH
       B(K)=PH
       A(K)=AMP
  37   continue
C************************************************
C
C     PUNCH THE AMPLITUDES AND PHASES
C
      WRITE(7,*) T0
      WRITE(7,*) A0
      WRITE(55,*) A0
      DO 704 J=1,M
         WRITE(55,*) A(J),B(J)
         WRITE(7,*) A(J),B(J)
 704  CONTINUE
C
      CLOSE(7)
C
C     CALCULATION OF THE DISPERSION
C       
       DISP=0.0D+00
       DO 2000 I=1,N
       S=A0
       TIME=T(I)
       DO 2001 K=1,M
       PHASE=F(K)*TIME
       S=S+A(K)*DSIN(P7*PHASE+B(K))
 2001  CONTINUE
       S=X(I)-S
       DISP=DISP+S*S
 2000  CONTINUE
C
C      UNBIASED ESTIMATION OF THE DISPERSION
C
       DISP=DSQRT(DISP/DFLOAT(N-2*M-1))
C
C      PRINT OUT FOURIER DECOMPOSITION IN A NICE FORMAT
C
      OPEN(3,FILE='fdec.dat')
C
      WRITE(3,823) smass,slu,strad,steff,xc,zc
 823  format(1x,f5.2,2x,f5.1,2x,f5.2,2x,f6.0,2x,f5.2,2x,f6.4)
      WRITE(3,826) P0,T0,A0,N,DISP
 826  FORMAT(1X,F13.10,3X,F7.4,1X,F7.3,1X,I7,1X,F7.4,1X,F5.2)
C
      kk=m+1
      do 1137 k=kk,50
      a(k)=0.0d0
      b(k)=0.0d0
 1137 continue
C
      WRITE(3,814) A(1),B(1),A(2),B(2),A(3),B(3),A(4),B(4),A(5),B(5)
      WRITE(3,814) A(6),B(6),A(7),B(7),A(8),B(8),A(9),B(9),A(10),B(10)
      WRITE(3,814) A(11),B(11),A(12),B(12),A(13),B(13),A(14),B(14),
     &A(15),B(15)
 814  FORMAT(10(1X,F7.4))
      sum = 0.0d0
      do 58 i = 1, Npts
       sum = A0
       do 59 j = 1, MORD
         t5 = 2.0d0*3.141d0*t1(i)*dble(j) + B(j)
         sum = sum + (A(j)*sin(t5))
  59   continue
  58   continue
c 58   write(66,*)t1(i),rl(i),sum
         
c
c     Compute phase differences and amplitude ratios
c
      i=1
      write(*,914) i,a(1),p0
 914  format(1x,i3,3x,f5.3,4x,f5.3)
      do 910 i=2,5
      ar(i)=a(i)/a(1)
      pr(i)=b(i)-i*b(1)
      if(ian.ne.1) go to 912
      if(i.eq.2) pr(i)=pr(i)-3.141593
      if(i.eq.4) pr(i)=pr(i)-3.141593
 912  continue
      if(pr(i).gt.0.0) go to 911 
      pr(i)=pr(i)+6.283185
      go to 912
 911  continue
      write(*,913) i,ar(i),pr(i)
 913  format(1x,i3,2(3x,f5.2))
 910  continue
c
      CLOSE(3)
C
      RETURN
      END
C
       SUBROUTINE MATINV(M,G)
C      MATRIX INVERTATION
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      INPUT: G(M3,M3), OUTPUT: G(M3,M3)
       DOUBLE PRECISION G(200,201)
       M3=2*M+1
       DO 11 I=1,M3
       W=G(I,I)
       Q=1.0D+00/W
       DO 12 K=1,M3
       G(I,K)=Q*G(I,K)
       IF(K.NE.I) GO TO 12
       G(I,K)=Q
  12   CONTINUE
       DO 13 J=1,M3
       IF(J.EQ.I) GO TO 13
       TT=G(J,I)
       Q1=-TT/W
       DO 14 K=1,M3
       G(J,K)=G(J,K)-TT*G(I,K)
       IF(K.NE.I) GO TO 14
       G(J,K)=Q1
  14   CONTINUE
  13   CONTINUE
  11   CONTINUE
       RETURN
       END
C
       SUBROUTINE NORM(N,M,T,X,F,G)
C      CALCULATION OF THE ELEMENTS OF THE NORMAL MATRIX AND R.H.S.
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DOUBLE PRECISION G(200,201),T(5000),X(5000),F(200)
       P7=8.0D0*DATAN(1.0D0)
       M3=2*M+1
       M4=M3+1
       L=0
       G(1,1)=FLOAT(N)
       DO 1 J=2,M3,2
       S=0.0D+00
       C=0.0D+00
       L=L+1
       J1=J+1
       FL=F(L)
       DO 2 I=1,N
       PH=T(I)*FL
       PH=P7*PH
       S=S+DSIN(PH)
  2    C=C+DCOS(PH)
       G(1,J)=C
  1    G(1,J1)=S
       L1=0
       DO 3 J=2,M3,2
       L1=L1+1
       J1=J+1
       L2=L1-1
       FL1=F(L1)
       DO 4 K=J,M3,2
       K1=K+1
       L2=L2+1
       FL2=F(L2)
       S1=0.0D+00
       S2=0.0D+00
       C1=0.0D+00
       C2=0.0D+00
       DO 5 I=1,N
       PH1=T(I)*FL1
       PH1=P7*PH1
       PH2=T(I)*FL2
       PH2=P7*PH2
       SPH1=DSIN(PH1)
       CPH1=DCOS(PH1)
       SPH2=DSIN(PH2)
       CPH2=DCOS(PH2)
       C1=C1+CPH1*CPH2
       S1=S1+CPH2*SPH1
       C2=C2+SPH2*CPH1
  5    S2=S2+SPH1*SPH2
       G(J,K)=C1
       G(J,K1)=C2
       G(J1,K)=S1
  4    G(J1,K1)=S2
  3    CONTINUE
       S=0.0D+00
       DO 6 I=1,N
  6    S=S+X(I)
       G(1,M4)=S
       L=0
       DO 7 J=2,M3,2
       S=0.0D+00
       C=0.0D+00
       L=L+1
       J1=J+1
       FL=F(L)
       DO 8 I=1,N
       PH=T(I)*FL
       PH=P7*PH
       S=S+X(I)*DSIN(PH)
  8    C=C+X(I)*DCOS(PH)
       G(J,M4)=C
  7    G(J1,M4)=S
       DO 9 J=1,M3
       DO 10 K=J,M3
  10   G(K,J)=G(J,K)
   9   CONTINUE
       RETURN
       END
C
      SUBROUTINE KURINT(cfeh,clogg,cteff,index,value)
c
c     THIS ROUTINE INTERPOLATES SOME TABULATED QUANTITY IN THE
c     ATMOSPHERE MODEL (E.G. THAT OF KURUCZ) 
c
c     Input:     cfeh  --- [Fe/H]
c                clogg --- Log g
c                cteff --- T_eff
c                index ----------- 0 ... just read in the atm. model
c                      ----------- 1 ... Bolometric Correction
c                      ----------- 2 ... B-V
c                      ----------- 3 ... U-B
c
c     Output:    value
c
c     Possible  ranges of [Fe/H], Log g and T_eff are given 
c     with the parameters (femin,femax), (gmin,gmax), (tmin,tmax)
c
c..................................................................
c
      implicit real*8 (a-h,o-z)
c
      character*70 aline
      dimension feh(5000),gl(5000),tef(5000)
      dimension ub(5000),bv(5000),bc(5000)
      dimension vr(5000),vi(5000),vk(5000)
      dimension x(2000),y(2000),z(2000)
      dimension grav(50),fehc(20),val(100),vf(100)
      dimension fehd(19)
c
      data fehd /-5.0,-4.5,-4.0,-3.5,-3.0,-2.5,-2.0,-1.5,
     &-1.0,-0.5,-0.3,-0.2,-0.1,0.0,0.1,0.2,0.3,0.5,1.0/
c
c     Some constants (number of grid points, step size, ... )
c
      ng=11
      g1=0.0d0
      dg=0.5d0
      nf=19
      t2=35000.0d0
c
c     [Fe/H], Log g and T_eff ranges
c
      tmin=4000.0d0
      tmax=9000.0d0
      gmin=0.0d0
      gmax=4.5d0
      femin=-3.0d0
      femax=0.0d0
c
      itab = 1
c
c     if itab=1 :  read in  'kurucz92.col'  
c                           from /users/kovacs/kurucz.dir/
c
c     if itab=2 :  read in  'cd13.tab'
c                           from /d1/jurcsik/kurucz/
c
c     !!!!!!!!!!!  WITH  itab=2  IT DOES NOT WORK AT THE MOMENT  !!!!!
c   
c
      if(index.ne.0) go to 1000
c
      if(itab.eq.2) go to 2000
c
c*********************************
c     Read in Kurucz'92 data     *
c*********************************
c
c.... Read in  'kurucz92.col'  from /users/kovacs/kurucz.dir/
c   
      open(1,file='/scratch/kovacs/kurucz92.col')
c 
      n=0
      read(1,503) aline
      read(1,503) aline
 503  format(a50)
c
      do 510 j=1,nf
      feoh=fehd(j)
      read(1,503) aline
c
cc      write(*,553) aline,feoh
cc 553  format(a25,3x,f5.2)
c
 514  continue
      do 511 k=1,ng
      grav(k)=g1+dg*dfloat(k-1)
c
      read(1,*) teff,bcv,umb,bmv
c
      if(feoh.lt.femin) go to 511
      if(feoh.gt.femax) go to 511
      if(teff.lt.tmin) go to 511
      if(teff.gt.tmax) go to 511
      if(grav(k).lt.gmin) go to 511
      if(grav(k).gt.gmax) go to 511
      n=n+1
      feh(n)=feoh
      tef(n)=teff
      gl(n)=grav(k)
      bv(n)=bmv
      ub(n)=umb
      bc(n)=bcv
c
c      write(*,513) n,feh(n),tef(n),gl(n),bv(n),ub(n)
 513  format(2x,i5,3x,f5.2,2x,f6.0,2x,f4.2,2(2x,f5.3))
c
 511  continue
c
      if(dabs(teff-t2).gt.0.1d0) go to 514
c
 510  continue
c
      close(1)
      nm=n
c
      write(*,*) nm
c
      return
c
c***********************************************************
c
 2000 continue
c
c
c*********************************
c     Read in Kurucz'92 data     *
c*********************************
c
c.... Read in  'cd13.tab'  from /d1/jurcsik/kurucz/
c   
      open(1,file='/d1/jurcsik/kurucz/cd13.tab')
c 
      n=0
      read(1,2001) aline
 2001 format(a70)
 2002 continue
      n=n+1
      read(1,*,end=2010) tef(n),gl(n),feh(n),bc(n),vv,
     &ub(n),bv(n),vr(n),vi(n),vk(n)
c
cc      write(*,513) n,feh(n),tef(n),gl(n),bv(n),ub(n)
c
      go to 2002
 2010 continue
      close(1)
      nm=n-1
      write(*,*) nm
      return
c
c***********************************************************

c
 1000 continue
c
c**************************
c     Interpolation       *
c**************************
c
c     Check ranges
c
      igo=1
      if(cfeh.lt.femin) igo=0
      if(cfeh.gt.femax) igo=0
      if(clogg.lt.gmin) igo=0
      if(clogg.gt.gmax) igo=0
      if(cteff.lt.tmin) igo=0
      if(cteff.gt.tmax) igo=0
      if(igo.eq.0) write(*,*) 'Out of range error!'
c
      if(igo.eq.1) go to 2100
c
      write(*,*) cfeh,clogg,cteff
      read(*,*) qq
c
cccccc      if(igo.eq.0) stop
c
 2100 continue
c
c     Choose the variable to be interpolated
c     
      do 80 i=1,nm
      if(index.eq.1) z(i)=bc(i)
      if(index.eq.2) z(i)=bv(i)
      if(index.eq.3) z(i)=ub(i)
 80   continue 
c
c     Start interpolation (quadratic in each variable)
c
c.... [Fe/H] loop
c
      jf=0
c
      do 100 nk=1,nf
c
      feoh=fehd(nk)
      if(feoh.lt.femin) go to 100
      if(feoh.gt.femax) go to 100
      jf=jf+1
      fehc(jf)=feoh
c
c**** Interpolate in the  (T_eff,log g)  plane  ---  [Fe/H] is fixed
c
      jg=0
c
      do 101 ni=1,ng
c
      gravi=g1+dg*dfloat(ni-1)
      if(gravi.lt.gmin) go to 101
      if(gravi.gt.gmax) go to 101
      jg=jg+1
      grav(jg)=gravi
c
c**** Interpolate in T_eff
c
      n=0
      eps=1.0d-5
c
      do 102 i=1,nm
c
      if(dabs(feoh-feh(i)).gt.eps) go to 102
      if(dabs(gravi-gl(i)).gt.eps) go to 102
c
      n=n+1
      x(n)=tef(i)
      y(n)=z(i)
c
 102  continue
c
      if(n.lt.3) go to 101
      CALL INTIM(X,Y,n,cteff,val(jg))
c
 101  continue
c
c**** Interpolate in log g
c
      n=jg
      do 200 i=1,n
      x(i)=grav(i)
      y(i)=val(i)
 200  continue
c
      if(n.lt.3) go to 100
      CALL INTIM(X,Y,n,clogg,vf(jf))
c
 100  continue
c
c**** Interpolate in [Fe/H]
c
      n=jf
      do 300 i=1,n
      x(i)=fehc(i)
      y(i)=vf(i)
 300  continue
c
      if(n.lt.3) stop
      CALL INTIM(X,Y,n,cfeh,value)
c
      return
      end
C
      SUBROUTINE INTIM(X,Y,NN,XX,YY)
C
C     This routine performs quadratic interpolation in the the sequence 
C     (x(i),y(i)), i=1,2,...,nn
C
C     The interpolated value is computed at    XX
C     The interpolated value is                YY 
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION Y(2000),X(2000)
C
      irev=0
      if(x(1).gt.x(nn)) irev=1
      if(irev.eq.1) go to 10
C
      if(x(1).gt.xx) write(*,*) 'Out of range error!'
cccccc      if(x(1).gt.xx) stop
      if(x(nn).lt.xx) write(*,*) 'Out of range error!'
cccccc      if(x(nn).lt.xx) stop
      go to 11
 10   continue
      if(x(1).lt.xx) write(*,*) 'Out of range error!'
      if(x(1).lt.xx) stop
      if(x(nn).gt.xx) write(*,*) 'Out of range error!'
      if(x(nn).gt.xx) stop
 11   continue
C
C     Compute the array indices
C
      d=1.0d30
      do 1 i=1,nn
      dd=dabs(xx-x(i))
      if(dd.gt.d) go to 1
      d=dd
      j1=i-1
      j2=i
      j3=i+1
      if(j1.lt.0) go to 2
      if(j3.gt.nn) go to 3
      go to 1
 2    continue
      j1=1
      j2=2
      j3=3
      go to 1
 3    continue
      j1=nn-2
      j2=nn-1
      j3=nn
 1    continue
C      
      AA=(XX-X(J2))*(XX-X(J3))/(X(J1)-X(J2))/(X(J1)-X(J3))
      BB=(XX-X(J1))*(XX-X(J3))/(X(J2)-X(J1))/(X(J2)-X(J3))
      CC=(XX-X(J1))*(XX-X(J2))/(X(J3)-X(J1))/(X(J3)-X(J2))
C
      YY=Y(J1)*AA+Y(J2)*BB+Y(J3)*CC
C
      RETURN
      END