C   prog2a.f -- Solve viscous Burgers Equation
C   Written Apr 6, 2003 by Eric Olson for Mathematics 484
C
C   For information about the Discrete Fourier Transform please see
C   Chapter 12 in Numerical Recipies by Press, Flannery, Teukolsky and
C   Vetterling.  Information about the heat and Burgers equations may
C   be found in chapters 13 and 14 in A First Course in the Numerical
C   Analysis of Differential Equations by Iserles.

      program main
      implicit none

      include "param2a.h"
C   param2a.h contains the following lines uncommented
C   Note N must be a power of 2
C     integer N
C     parameter(N=8)

      complex*16 u
      real*8 t
      integer i
      common /soln/u(N)
      write(*,10)
10    format(1X,'#prog2a.f -- Compute viscous Burger''s Equation')
      call doinit(u)
      call vectwr(u,0.0D0)
      t=0.0
      do i=1,8
        call euler(u,0.125D0,1024)
        call vectwr(u,i*0.125D0)
      enddo
      stop
      end

      subroutine doinit(u)
      implicit none
      include "param2a.h"
      complex*16 u
      real*8 x,xx,m_pi
      integer i
      dimension u(N)
      data m_pi/3.14159265358979323D0/
      do i=1,N
        x=2*M_PI*(i-1)/N
        xx=M_PI-x 
        xx=xx*xx
        u(i)=8*exp(-xx)
      enddo
      call fft(u,N,1)
      end

      subroutine force(f,u)
      implicit none
      include "param2a.h"
      complex*16 uu
      common /temp/uu(N)
      complex*16 f,u,ik
      real*8 a,b
      dimension f(N),u(N)
      integer i,k
      do i=1,N
        uu(i)=u(i)
      enddo
      call fft(uu,N,-1)
      do i=1,N
         uu(i)=DCMPLX(DBLE(uu(i))*DBLE(uu(i)),0.0D0)
      enddo
      call fft(uu,N,1)
      do i=1,N
        if(2*(i-1).gt.N) then
            k=i-1-N
        else
            k=i-1
        endif
        ik=DCMPLX(0.0D0,DBLE(k))
        f(i)=ik*(ik*u(i)-uu(i))
        if(2*(i-1).eq.N) then
            f(i)=DCMPLX(DBLE(f(i)),0.0D0)
        endif
      enddo
      end

      subroutine euler(u,tn,m)
      implicit none
      include "param2a.h"
      complex*16 u,f
      real*8 tn,dt
      integer m,i,j
      dimension u(N)
      common /fvect/f(N)
      dt=tn/m
      do j=1,m
        call force(f,u)
        do i=1,N
            u(i)=u(i)+dt*f(i)
        enddo
      enddo
      end

      subroutine vectwr(u,t)
      implicit none
      include "param2a.h"
      complex*16 u,uu
      common /temp/uu(N)
      real*8 t,m_pi
      integer i,k
      dimension u(N)
      data m_pi/3.14159265358979323D0/
      do i=1,N
        uu(i)=u(i) 
      enddo
      call fft(uu,N,-1)
      write(*,20) t
20    format(1X,//,1X,'#At time ',F8.5,' ...')
      do i=1,N
        if(2*(i-1).gt.N) then
            k=i-1-N
        else
            k=i-1
        endif
        write(*,30) 2*M_PI*(i-1)/N,uu(i),k,u(i)
30      format(1X,3E14.5,I6,2E14.5)
      enddo
      end
