******************************************************************
  
  subroutine pda_rfftf(n,r,wsave)
  
  ******************************************************************
  
  subroutine pda_rfftf computes the fourier coefficients of a real
  perodic sequence (fourier analysis). the transform is defined
  below at output parameter r.
  
  input parameters
  
  n       the length of the array r to be transformed.  the method
          is most efficient when n is a product of small primes.
          n may change so long as different work arrays are provided
  
  r       a real array of length n which contains the sequence
          to be transformed
  
  wsave   a work array which must be dimensioned at least 2*n+15.
          in the program that calls pda_rfftf. the wsave array must be
          initialized by calling subroutine pda_rffti(n,wsave) and a
          different wsave array must be used for each different
          value of n. this initialization does not have to be
          repeated so long as n remains unchanged thus subsequent
          transforms can be obtained faster than the first.
          the same wsave array can be used by pda_rfftf and pda_rfftb.
  
  
  output parameters
  
  r       r(1) = the sum from i=1 to i=n of r(i)
  
          if n is even set l =n/2   , if n is odd set l = (n+1)/2
  
            then for k = 2,...,l
  
               r(2*k-2) = the sum from i = 1 to i = n of
  
                    r(i)*cos((k-1)*(i-1)*2*pi/n)
  
               r(2*k-1) = the sum from i = 1 to i = n of
  
                   -r(i)*sin((k-1)*(i-1)*2*pi/n)
  
          if n is even
  
               r(n) = the sum from i = 1 to i = n of
  
                    (-1)**(i-1)*r(i)
  
   *****  note
               this transform is unnormalized since a call of pda_rfftf
               followed by a call of pda_rfftb will multiply the input
               sequence by n.
  
  wsave   contains results which must not be destroyed between
          calls of pda_rfftf or pda_rfftb.