subroutine fourier(j_dim,j_max,u,u_hat)

c     Returns the Fourier coefficients u_hat(k) of u(x) for problems
c     in which u(x) has Dirichlet boundary conditions at x=0 and
c     x=l.

c     In particular, the returned coefficients are such that
c       u(j) = u(0) (1 - x/l) + u(j_max) x/l
c              + {sum from k=1 to k=j_max-1) u_hat(k) sin(k pi x(j)/l)

c     In other words, the u_k are the Fourier sine series coefficients.
c     They are equal to the twice the imaginary parts of the complex
c     Fourier series coefficients (the real parts of the complex
c     coefficients are zero.)

      implicit none

c input

c     declared maximum value of the mesh index j
      integer j_dim

c     actual maximum value of the mesh index j
      integer j_max

c     one-dimensional array of u-values u(j)
      real u(0:j_dim)

c output

c     values of the Fourier coefficients u_hat(k)
      real u_hat(0:j_dim)

c local

c     mesh point index
      integer j

c     wave number
      integer k

c     u value minus the linear part
      real du

c     the values of cos(pi x(j)/l) and sin(pi x(j)/l)
      real cos1,sin1
c     the values of cos(k pi x(j)/l) and sin(k pi x(j)/l)
      real cosk,sink

c     temporary
      real tmp

c     pi  (the correct value 3.14... of pi will be set later)
      real pi
      save pi
      data pi/0./

c     compute pi if it has not been done yet
      if(pi.eq.0.)pi=4.*atan(1.)

c     initialize all coefficients to zero

      do 100 k=0,j_max
         u_hat(k)=0.
 100  continue

c     add the contributions of each point to the coefficients

      do 200 j=1,j_max-1

c        u value minus the linear part
         du=u(j)-u(0)-(u(j_max)-u(0))*j/j_max

c        initialize the needed cosines and sines
         cos1=cos(pi*j/j_max)
         sin1=sin(pi*j/j_max)
         cosk=1.
         sink=0.

c        add the contributions of point j to all coefficients u_hat(k)
         do 150 k=1,j_max-1

c           compute cos(k pi x(j)/l) and sin(k pi x(j)/l) using trig
            tmp=cosk
            cosk=cosk*cos1-sink*sin1
            sink=sink*cos1+tmp*sin1

c           add the contribution
            u_hat(k)=u_hat(k)+du*sink*2./j_max

 150     continue

 200  continue

      return
      end