parameter (N=64, NN=2*N+1)
       implicit real*8(a-h,o-z)
       real x(0:NN), uapp(0:NN), exact(0:NN)
       complex uhat(0:NN)
       complex w, A, y, ss, ss1, ss2, A1

       pi=4.0*atan(1.0)
       w=cmplx(0, 1)
       pa=1.9

       do k=0, NN-1
           k1=k-N
           p1=(k1+0.5)*pa
           p2=(k1-0.5)*pa
           s1=(-cos(p1)/(k1+0.5)+cos(p2)/(k1-0.5))/(2*pi)
           s2=(sin(p1)/(k1+0.5)-sin(p2)/(k1-0.5))/(2*pi)
           uhat(k)=cmplx(s1,s2)
       enddo

       ss=(N-1)*uhat(N+N-1)/(N*uhat(N+N))
       y=-w*log(ss)
       y=real(y)
       err1=abs(real(y)-pa)
       print *, err1, '   error in position'

       ss=w*(N-1)*y
       A=w*(N-1)*uhat(N+N-1)*exp(ss)
       A=real(A)
       err2=abs(-sin(0.95)/pi-real(A))
       print *, err2, '   error in strength'

      do j=0, 2*N
         x(j)=2*pi*j/(2*N+1)
          if (x(j) .le. real(y)) then
             F=-real(A)*x(j)
             exact(j)=sin(0.5*x(j))
          else
             F=real(A)*(2*pi-x(j))
             exact(j)=-sin(0.5*x(j))
          endif
         
          sum=0.0
          do k=1, N
             tau=exp(-15*log(10.0)*(1.0*k/N)**8)
             ss1=-w*k*y
             ss2=w*k*x(j)
             ss=(uhat(k+N)+w*A*exp(ss1)/k)*exp(ss2)
             sum=sum+tau*ss
             k1=-k
             ss1=-w*k1*y
             ss2=w*k1*x(j)
             ss=(uhat(k1+N)+w*A*exp(ss1)/k1)*exp(ss2)
             sum=sum+tau*ss
          enddo

          uapp(j)=real(uhat(N)-A*(pi-y)+sum+F)
c          print *, uapp(j), exact(j), abs(uapp(j)-exact(j))
       enddo

       write (17, *) 'x=['
       do j=0, NN-1
          write (17, *) x(j)
       enddo
       write (17, *) '];'

        write (18, *) 'y=['
       do j=0, NN-1
          write (18, *) real(uapp(j))
       enddo
       write (18, *) '];'

       write (19, *) 'y=['
       do j=0, NN-1
          write (19, *) abs(uapp(j)-exact(j))
       enddo
       write (19, *) '];'

c          write (19, *) abs(uapp(j)-exact(j))

          
       stop
       end