!-----------------------------------------------------------------------
!   Laplace's equation: Gauss-Seidel method 伊理・藤野(1985) p.131
!
!   coded by Shigeo Yoden on Sept. 13, 2000
!-----------------------------------------------------------------------
program laplace_gs

  use dcl        ! DCL Fortran90 (dcl-5.1/dcl-f90)

  integer, parameter :: imax=20
  real,    parameter :: pi=3.141593, dx=pi/imax
  real, dimension(0:imax,0:imax) :: vold, vnew, van
  character :: ct*30, cl*30, csgi*1

!-- 境界条件 & 初期推定値 ----
    do i=0, imax
      vold(i,imax) = sin(i*dx)
    end do
    do j=0, imax-1
      vold(:,j) = j*dx*vold(:,imax)
    end do
    vold(imax,:) = 0

    vnew = vold
   
!!!    write(*,*) 'epsilon=?'
!!!    read (*,*)  epsilon
    epsilon  = 0.0001
    test = 2*epsilon
    lmn  = 0

! execute G-S iteration and calculate TEST where
! TEST is the maximum component of (vnew - vold)
    do while(test > epsilon)      ! convergence check
       test = 0.
       do j=1, imax-1
         do i=1, imax-1
           vnew(i,j) = (vold(i-1,j) + vold(i,j-1)      &
                      + vold(i+1,j) + vold(i,j+1)) /4
           diff = abs(vnew(i,j)-vold(i,j))
           if(diff > test) test = diff
         end do
       end do
       vold = vnew
       lmn = lmn + 1
    end do

    write(*,*) 'iteration # = ', lmn, '  DIFmax = ', test

!-- 解析解 ----
    do j=0, imax
      do i=0, imax
           van(i,j) = sinh(j*dx)*sin(i*dx) /sinh(pi)
      end do
    end do
    van(imax,:) = 0

!-- グラフ化 ----
    call DclOpenGraphics()
    call DclNewFrame

    call DclSetWindow  (0.0, 1.0, 0.0, 1.0)
    call DclSetViewPort(0.2, 0.8, 0.2, 0.8)
    call DclSetTransFunction

    call DclSetTitle('x', 'y', csgi(194)//' '//csgi(167),  &
                               csgi(194)//' '//csgi(167))
    call DclDrawScaledAxis

    call DclSetParm('ENABLE_CONTOUR_MESSAGE', .false.)
    call DclSetParm('ENABLE_NEGATIVE_CONTOUR', .false.)
    call DclDrawContour(van)           ! 解析解(実線)
    call DclSetParm('POSITIVE_CONTOUR_TYPE', 3)
    call DclDrawContour(vnew)          ! 数値解(点線)

      ct = 'iteration # =?????'
      write(ct(14:18),'(i5)') lmn
    call DclDrawTitle('t', ct, position=-1.)
      cl = '=??.??E???s'
      write(cl(2:10),'(e9.2)') epsilon
      ct = csgi(156) // cl  
    call DclDrawTitle('t', ct, position=-0.8)
    call DclDrawTitle('t', ' Gauss-Seidel method', position=-1.)
    call DclDrawTitle('t', 'Laplace'//csgi(178)//'s Equation:', &
                                                   position=-1.)
    call DclCloseGraphics
end program laplace_gs
!-----------------------------------------------------------------------
!   Laplace's equation: Successive overrelaxation(SOR) method
!         Johnson and Riess, 1982: "Numerical Analysis(2nd edition)",
!                     Addison-Wesley Publishing Company, p.472
!
!   coded by Shigeo Yoden on Sept. 13, 2000
!-----------------------------------------------------------------------
program laplace_sor

  use dcl        ! DCL Fortran90 (dcl-5.1/dcl-f90)

  integer, parameter :: imax=20
  real,    parameter :: pi=3.141593, dx=pi/imax
  real, dimension(0:imax,0:imax) :: vold, vnew, van
  character :: ct*30, cl*30, csgi*1

!-- 境界条件 & 初期推定値 ----
    do i=0, imax
      vold(i,imax) = sin(i*dx)
    end do
    do j=0, imax-1
      vold(:,j) = j*dx*vold(:,imax)
    end do
    vold(imax,:) = 0

    vnew = vold
   
!!!    write(*,*) 'epsilon=?'
!!!    read (*,*)  epsilon
    epsilon  = 0.001
    test = 2*epsilon
    lmn  = 0

!-- SOR ----
    omega = 2/(1 + sin(dx))

! execute SOR iteration and calculate TEST where
! TEST is the maximum component of (vnew - vold)
    do while(test > epsilon)      ! convergence check
       test = 0.
       do j=1, imax-1
         do i=1, imax-1
           vnew(i,j) = vold(i,j) + omega*(vnew(i-1,j) + vnew(i,j-1)  &
                          + vold(i+1,j) + vold(i,j+1) -4*vold(i,j))/4
           diff = abs(vnew(i,j)-vold(i,j))
           if(diff > test) test = diff
         end do
       end do
       vold = vnew
       lmn = lmn + 1
    end do

    write(*,*) 'iteration # = ', lmn, '  DIFmax = ', test

!-- 解析解 ----
    do j=0, imax
      do i=0, imax
           van(i,j) = sinh(j*dx)*sin(i*dx) /sinh(pi)
      end do
    end do
    van(imax,:) = 0

!-- グラフ化 ----
    call DclOpenGraphics()
    call DclNewFrame

    call DclSetWindow  (0.0, 1.0, 0.0, 1.0)
    call DclSetViewPort(0.2, 0.8, 0.2, 0.8)
    call DclSetTransFunction

    call DclSetTitle('x', 'y', csgi(194)//' '//csgi(167),  &
                               csgi(194)//' '//csgi(167))
    call DclDrawScaledAxis

    call DclSetParm('ENABLE_CONTOUR_MESSAGE', .false.)
    call DclSetParm('ENABLE_NEGATIVE_CONTOUR', .false.)
    call DclDrawContour(van)           ! 解析解(実線)
    call DclSetParm('POSITIVE_CONTOUR_TYPE', 3)
    call DclDrawContour(vnew)          ! 数値解(点線)

      ct = 'iteration # =?????'
      write(ct(14:18),'(i5)') lmn
    call DclDrawTitle('t', ct, position=-1.)
      cl = '=??.??E???s'
      write(cl(2:10),'(e9.2)') epsilon
      ct = csgi(156) // cl  
    call DclDrawTitle('t', ct, position=-0.8)
    call DclDrawTitle('t', ' SOR method', position=-1.)
    call DclDrawTitle('t', 'Laplace'//csgi(178)//'s Equation:',  &
                                          position=-1.)
    call DclCloseGraphics
end program laplace_sor
!-----------------------------------------------------------------------
!   Poisson's equation: Gauss-Seidel method 伊理・藤野(1985) p.131
!
!   coded by Shigeo Yoden on Sept. 13, 2000
!-----------------------------------------------------------------------
program poisson_gs

  use dcl        ! DCL Fortran90 (dcl-5.1/dcl-f90)

  integer, parameter :: imax=60
  real,    parameter :: pi=3.141593, dx=pi/imax
  real, dimension(0:imax,0:imax) :: vold, vnew, rho
  character :: ct*30, cl*30, csgi*1

!-- 境界条件 & 初期推定値 ----
    vold = 0
    vnew = 0
    rho  = 0
    rho(  imax/3,  imax/3) =   1
    rho(2*imax/3,2*imax/3) = - 2
!!!    rho(imax/2,imax/2) = - 1

!!!    write(*,*) 'epsilon=?'
!!!    read (*,*)  epsilon
    epsilon  = 0.000001
    test = 2*epsilon
    lmn  = 0

! execute G-S iteration and calculate TEST where
! TEST is the maximum component of (vnew - vold)
    do while(test > epsilon)      ! convergence check
       test = 0.
       do j=1, imax-1
         do i=1, imax-1
           vnew(i,j) = (   vnew(i-1,j) + vnew(i,j-1)      &
                         + vold(i+1,j) + vold(i,j+1)      &
                                       - rho(i,j)*dx**2 ) /4
           diff = abs(vnew(i,j)-vold(i,j))
           if(diff > test) test = diff
         end do
       end do
       vold = vnew
       lmn = lmn + 1
    end do

    write(*,*) 'iteration # = ', lmn, '  DIFmax = ', test

!-- グラフ化 ----
    call DclOpenGraphics()
    call DclNewFrame

    call DclSetWindow  (0.0, 1.0, 0.0, 1.0)
    call DclSetViewPort(0.2, 0.8, 0.2, 0.8)
    call DclSetTransFunction

    call DclSetTitle('x', 'y', csgi(194)//' '//csgi(167),  &
                               csgi(194)//' '//csgi(167))
    call DclDrawScaledAxis

    call DclDrawContour(vnew)

      ct = 'iteration # =?????'
      write(ct(14:18),'(i5)') lmn
    call DclDrawTitle('t', ct, position=-1.)
      cl = '=??.??E???s'
      write(cl(2:10),'(e9.2)') epsilon
      ct = csgi(156) // cl  
    call DclDrawTitle('t', ct, position=-0.8)
    call DclDrawTitle('t', ' Gauss-Seidel method', position=-1.)
    call DclDrawTitle('t', 'Poisson'//csgi(178)//'s Equation:', position=-1.)
 
    call DclCloseGraphics
end program poisson_gs