!-----------------------------------------------------------------------
!
!  saite.f90 :   Loesung der 1D-Wellengleichung
!
!                              2
!                u  (x,t)  =  c * u  (x,t)
!                 tt               xx
!
!                auf dem Intervall [0,L=1] mit einem
!                expliziten Verfahren 2. Ordnung.
!
!                Anfangsbedingungen:
!
!                u(x,0)  = f(x)
!                u (x,0) = g(x)
!                 t
!
!                Randbedingungen (feste Enden):
!
!                u(0,t) = u(L,t) = 0
!
!  u(x,t) ...... Auslenkung
!  c ........... Phasengeschwindigkeit
!  N + 1 ....... Anzahl der Gitterpunkte im Ortsgitter x(0),...,x(N)
!  dx .......... Schrittweite im Ortsgitter
!  dt .......... Zeitschritt
!  tmax ........ Integrationsintervall [0,tmax]
!  iout ........ Ausgabe alle "iout" Zeitschritte
!  saite.dat ... Ausgabefile
!
!-----------------------------------------------------------------------

module constants_m
   implicit none
   integer, parameter, public ::                               &
            single = kind(1.0),                                &
            double = selected_real_kind(2*precision(1.0_single))
end module constants_m

!-----------------------------------------------------------------------

module global_m
   use constants_m
   implicit none
   private
   real(kind=double), public :: c
end module global_m

!-----------------------------------------------------------------------

module saite_sub_m

   use constants_m
   use global_m
   implicit none
   private
   public :: f, g, g1

contains

   !
   !  Anfangsbedingung  u(x,0) = f(x)
   !
   function f(x) result(f_result)
      real(kind=double), intent(in) :: x
      real(kind=double)             :: f_result
      f_result = exp( -100.0*(x-0.5)*(x-0.5) )
      return
   end function f

   !
   !  Anfangsbedingung  u (x,0) = g(x)
   !                     t
   !
   function g(x) result(g_result)
      real(kind=double), intent(in) :: x
      real(kind=double)             :: g_result
      g_result = 0.0 + (x-x)
      return
   end function g

   function g1(x) result(g1_result)
      real(kind=double), intent(in) :: x
      real(kind=double)             :: g1_result
      g1_result = -200.0*c*(x-0.5)*exp(-100.0*(x-0.5)*(x-0.5))
      return
   end function g1

end module saite_sub_m

!-----------------------------------------------------------------------

program saite

   use constants_m
   use global_m
   use saite_sub_m

   implicit none

   integer, parameter                :: N = 100
   real(kind=double), dimension(0:N) :: x, U, Uold, Unew
   real(kind=double)                 :: t, tmax, dt, dx, alpha, alpha2
   integer                           :: j, it, iout

   write( unit=*, fmt=* )
   write( unit=*, fmt="(a)", advance="no" ) " tmax [s]  (0.006): "
   read ( unit=*, fmt=* )  tmax
   write( unit=*, fmt="(a)", advance="no" ) " iout          (3): "
   read ( unit=*, fmt=* )  iout

   c  = 100.0_double
   dx = 1.0_double/N
   dt = dx/c

   alpha  = c*dt/dx
   alpha2 = alpha*alpha

   !  Ausgabefile :

   open( unit=1, file="saite.dat", status="replace", action="write" )

   !  Raeumliche Gitterpunkte :

   do j = 0, N
      x(j) = j * dx
   end do

   !  Anfangs- und Randbedingungen :

   t  = 0.0
   it = 0

   do j = 1, (N-1)
      Uold(j) = f(x(j))
   end do
   Uold(0) = 0.0
   Uold(N) = 0.0

   !  Ausgabe:
   do j = 0, N
      write( unit=1, fmt="(2(f8.5,tr1), es13.5)" )  x(j), t, Uold(j)
   end do
   write( unit=1, fmt="(/)" )

   !   (1)
   !  u      (erster Zeitschritt)
   !   j
   ! 

   t  = dt
   it =  1

   do j = 1, (N-1)
      U(j) = (1.0 - alpha2) * Uold(j)   +                            &
              0.5 * alpha2 * (Uold(j+1) + Uold(j-1))  +  dt * g1(x(j))
   end do
   U(0) = 0.0
   U(N) = 0.0


   !   (n+1)
   !  u
   !   j
   !

   do
      if ( modulo( it, iout ) == 0 ) then
         !  Ausgabe:
         do j = 0, N
            write( unit=1, fmt="(2(f8.5,tr1), es13.5)" )  x(j), t, U(j)
         end do
         write( unit=1, fmt="(/)" )
      end if

      do j = 1, (N-1)
         Unew(j) = 2.0 * (1.0 - alpha2) * U(j) +        &
                   alpha2 * (U(j+1) + U(j-1))  -  Uold(j)
      end do
      Unew(0) = 0.0
      Unew(N) = 0.0

      do j = 0, N
         Uold(j) = U(j)
         U(j)    = Unew(j)
      end do

      t  = t + dt
      it = it + 1

      if ( t > tmax ) then
         exit
      end if

   end do

   close( unit=1 )

   stop

end program saite

!-----------------------------------------------------------------------