/*@@
   @file      InitialData.F77
   @date      
   @author    Tom Goodale
   @desc 
              Initial data for the 3D Wave Equation
   @enddesc 
 @@*/

#include "cctk.h" 
#include "cctk_parameters.h"
#include "cctk_arguments.h"



 /*@@
   @routine    WaveToyF77_InitialData
   @date       
   @author     Tom Goodale
   @desc 
               Set up initial data for the wave equation
   @enddesc 
   @calls      
   @calledby   
   @history 
 
   @endhistory 

@@*/

      subroutine IDScalarWave_InitialData(CCTK_FARGUMENTS)

      implicit none

      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS

      INTEGER CCTK_Equals

      INTEGER    i,j,k
      CCTK_REAL  dt,omega, pi
      CCTK_REAL  dx,dy,dz

      pi = 4.0*atan(1.0)

c     Grid spacing shortcuts
c     ----------------------
      dx = cctk_delta_space(1)
      dy = cctk_delta_space(2)
      dz = cctk_delta_space(3)
      dt = cctk_delta_time

      omega = sqrt(kx**2+ky**2+kz**2)

      if (CCTK_Equals(initial_data,"plane").eq.1) then 

         do k=1,cctk_lsh(3)
            do j=1,cctk_lsh(2)
               do i=1,cctk_lsh(1)

                  phi(i,j,k)     = amplitude*cos(kx*x(i,j,k)+ky*y(i,j,k)
     &                             +kz*z(i,j,k)+omega*cctk_time)
                  phi_old(i,j,k) = amplitude*cos(kx*x(i,j,k)+ky*y(i,j,k)
     &                             +kz*z(i,j,k)+omega*(cctk_time-dt))
                  
               end do
            end do
         end do

      else if (CCTK_Equals(initial_data,"gaussian").eq.1) then 

         do k=1,cctk_lsh(3)
            do j=1,cctk_lsh(2)
               do i=1,cctk_lsh(1)

                  phi(i,j,k)     = amplitude*exp( -(sqrt(x(i,j,k)**2
     &                      +y(i,j,k)**2+z(i,j,k)**2)-radius)**2/sigma**2)
                  phi_old(i,j,k) = amplitude*exp( -(sqrt(x(i,j,k)**2
     &                      +y(i,j,k)**2+z(i,j,k)**2)-radius-dt)**2/sigma**2)

               end do
            end do
         end do

      else if (CCTK_Equals(initial_data, "box").eq.1) then

c        Use kx,ky,kz as number of modes in each direction.

         do k=1,cctk_lsh(3)
            do j=1,cctk_lsh(2)
               do i=1,cctk_lsh(1)

                  phi(i,j,k) = amplitude*sin(kx*(x(i,j,k)-0.5)*pi)*
     $                 sin(ky*(y(i,j,k)-0.5)*pi)*
     $                 sin(kz*(z(i,j,k)-0.5)*pi)*
     $                 cos(omega*cctk_time*pi)
                  
                  phi_old(i,j,k)= amplitude*sin(kx*(x(i,j,k)-0.5)*pi)*
     $                 sin(ky*(y(i,j,k)-0.5)*pi)*
     $                 sin(kz*(z(i,j,k)-0.5)*pi)*
     $                 cos(omega*(cctk_time-dt)*pi)


               end do
            end do
         end do

      end if

c     Apply symmetry boundary conditions
c     ----------------------------------
      call ApplySymmetry(cctkGH,"wavetoy::scalarevolve")

c     Synchronise
c     -----------
      call CCTK_SyncGroup(cctkGH,"wavetoy::scalarevolve")

      return
      end