#include <stdlib.h>
#include <math.h>

void fourier

 /********************************************************************\
 * Returns the Fourier coefficients u_hat(k) of u(x) for problems     *
 * in which u(x) has Dirichlet boundary conditions at x=0 and         *
 * x=l.                                                               *
 *                                                                    *
 * In particular, the returned coefficients are such that             *
 *   u(j) = u(0) (1 - x/l) + u(j_max) x/l                             *
 *          + {sum from k=1 to k=j_max-1) u_hat(k) sin(k pi x(j)/l)   *
 *                                                                    *
 * In other words, the u_k are the Fourier sine series coefficients.  *
 * They are equal to the twice the imaginary parts of the complex     *
 * Fourier series coefficients (the real parts of the complex         *
 * coefficients are zero.)                                            *
 \********************************************************************/

(
 int j_max,         // input: maximum value of the mesh index j
 float u[],         // input: one-dimensional array of u-values
 float u_hat[]      // output: values of the Fourier coefficients
)

{
  /************************ local variables *************************/

  // mesh point index
  int j;

  // wave number
  int k;

  // u value minus the linear part
  float du;

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

  // temporary
  float tmp;

  // pi  (the correct value 3.14... of pi will be set later)
  static float pi=0.;

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

  //initialize all coefficients to zero

  for (k=0; k<=j_max; k++) u_hat[k]=0.;

  //add the contributions of each point to the coefficients

  for (j=1; j<j_max; j++)
    {
      // u value minus the linear part
      du=u[j]-u[0]-(u[j_max]-u[0])*j/j_max;

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

      // add the contributions of point j to all coefficients u_hat(k)
      for (k=1; k<j_max; k++)
	{
	  // 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;

	  // add the contribution
	  u_hat[k]=u_hat[k]+du*sink*2./j_max;

	}

    }

  return;
}