program relax1
C
C  illustrate solution of boundary-value problem by relaxation
C  improved sweep
C
C  example taken from Koonin & Meredith
C
C   d^2 phi/dx^2 = -12 x^2
C
C  version 1: Dirichlet boundary conditions
C
C  CWJ SDSU May 2005
C

      implicit none

      integer size
      parameter(size = 40)
   
      real phi(0:size)
      real source(0:size)

      real dx,x
      real w	! relaxtion parameter

      real E
      integer nsweep,n

      dx = 1.0/float(size)

      call fillsource(size,dx,source)


      print*,' enter relaxation parameter w '
      read*,w

      print*,' Enter # of sweeps '
      read*,nsweep

C............. make initial guess.............

      do n = 0,size
         phi(n) = 0.0
      enddo

      do n = 1,nsweep
         call sweep(size,dx,w,source,phi)
         call energy(size,dx,source,phi,E)
         print*,' Sweep # ',n,' Energy = ',E

      enddo
      do n = 0,size
        write(17,*)n*dx,phi(n)
        write(18,*)n*dx,n*dx*(1-(n*dx)**3)
      enddo
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine sweep(size,dx,w,source,phi)
      implicit none
      integer size
      real phi(0:size),source(0:size)
 
      real dx
      integer i
      real w

      phi(0) = 0.
      phi(size) = 0.

      do i = 1,size-1
        phi(i) = (1.-w)*phi(i) + 
     & 0.5*w*( phi(i+1)+ phi(i-1)+dx**2*source(i))

      enddo

      return
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine energy(size,dx,source,phi,E)
C
C  computes the energy of phi 
C
      implicit none

      integer size
      real phi(0:size),source(0:size)
 
      real dx
      integer i
      real e

      E = 0.0
      do i = 1,size
        e = e + ( ( phi(i)-phi(i-1) )/dx)**2/2.
        e = e-source(i)*phi(i)
      
      enddo
      e = e*dx

      return
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine fillsource(size,dx,source)

      implicit none
      integer size

      real source(0:size)

      real dx
      real x
      integer i

      do i = 0,size
         x = i*dx
         source(i) = 12*x*x
      enddo
      return
      end
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC