FOURIER COEFFICIENTS

c     coeficientes de fourier
c     ncoef numero de coeficientes  determinar
      dimension coefcos(0:20), coefsin(0:20)
      ncoef =6
      a=1.
      period=1.
      call fourier(a,ncoef,period,coefcos, coefsin)
      print*,'a0,a1,a2,b0,b1,b2=',coefcos(0),coefcos(1),coefcos(2),
     $coefsin(0),coefsin(1),coefsin(2)
      stop
      end


      subroutine fourier(a,ncoef,period,coefcos, coefsin)
      dimension coefcos(0:20), coefsin(0:20)
      force(t)=A*t/period
      pi=2.*asin(1.)

      t1=-period/2.
      t2=period/2.
      w0=2.*pi/period
      nstep=100
      dt=(t2-t1)/float(nstep)
      do 10 i=0,ncoef
      coefcos(i)=0.
      coefsin(i)=0.
      do 20 j=1,nstep
      t=t1+dt*float(j)
      coefcos(i)= coefcos(i)+force(t-dt)*cos(float(i)*w0*(t-dt))+
     $force(t)*cos(float(i)*w0*t)
      coefsin(i)= coefsin(i)+force(t-dt)*sin(float(i)*w0*(t-dt))+
     $  force(t)*sin(float(i)*w0*t)
20    continue
      coefcos(i)= (2./period)*(dt/2.)*coefcos(i)
      coefsin(i)=(2./period)*(dt/2.)*coefsin(i)
10    continue
      return
      end
      
      Run
       a0,a1,a2,b0,b1,b2=  8.94069618E-009  1.46031374E-008
       -2.5331973E-008   0.  0.318205029 -0.158945516