program main

  real(4), parameter :: PI = 3.141592

  integer :: imax
  integer :: jmax
  integer :: kmax

  real(4) :: X0, DelX
  real(4) :: Y0, DelY
  real(4) :: Z0, DelZ
  real(4), allocatable :: x_X(:)
  real(4), allocatable :: y_Y(:)
  real(4), allocatable :: z_Z(:)

  integer       :: FU
  character(64) :: FileName = 'out/test1.vtk'

  integer :: i
  integer :: j
  integer :: k

  imax = 21
  jmax = 21
  kmax = 21

  allocate( x_X( imax ) )
  allocate( y_Y( jmax ) )
  allocate( z_Z( kmax ) )

  X0 = 0.0; DelX = 2.0*PI/real(imax)
  Y0 = 0.0; DelY = 2.0*PI/real(jmax)
  Z0 = 0.0; DelZ = 2.0*PI/real(kmax)
  do i = 1, imax
    x_X(i) = X0 + DelX * (i-1)
  end do
  do j = 1, jmax
    y_Y(j) = Y0 + DelY * (j-1)
  end do
  do k = 1, kmax
    z_Z(k) = Z0 + DelZ * (k-1)
  end do

  FU = 50
  open( FU, file = FileName, status='unknown' )

  write( FU, '(a)' ) '# vtk DataFile Version 1.0'
  write( FU, '(a)' ) trim( FileName )
  write( FU, '(a)' ) 'ASCII'
  write( FU, '(a)' ) 'DATASET STRUCTURED_POINTS'
  write( FU, '(a,3(a,i3))' ) 'DIMENSIONS', ' ', imax, ' ', jmax, ' ', kmax

  write( FU, '(a,3(a,f4.1))' ) 'ORIGIN', ' ', X0, ' ', Y0, ' ',  Z0
  write( FU, '(a,3(a,f4.1))' ) 'SPACING', ' ', DelX, ' ', DelY, ' ', DelZ
!!$
!!$  write( FU, '(a,i3)' ) 'POINTS ', imax*jmax*kmax
!!$  do k = 1, kmax
!!$    do j = 1, jmax
!!$      do i = 1, imax
!!$        write( FU, '(f5.2)' ) cos(x_X(i))*cos(y_Y(j))*z_Z(k)
!!$      end do
!!$    end do
!!$  end do

  write( FU, '(a,i8)' ) 'POINT_DATA ', imax*jmax*kmax

  write( FU, '(a)' ) 'SCALARS intensity float'
  write( FU, '(a)' ) 'LOOKUP_TABLE default'

  do k = 1, kmax
    do j = 1, jmax
      do i = 1, imax
        write( FU, '(f5.2)' ) cos(x_X(i))*cos(y_Y(j))*z_Z(k)
      end do
    end do
  end do

  close( FU )

end program main
