program Splash_
  use sep
  use util_

  implicit none
  type(myaxis)            :: ax,az
  real,    pointer        :: map(:,:)
  real                    :: xx,zz,mag,rad
  integer                 :: ix,iz
  real                    ::  x, z

  call sep_init()

  call hinaxis(az,1)
  call hinaxis(ax,2)

  call from_param("zz",zz,0.)
  call from_param("xx",xx,0.)
  call from_param("rad",rad,1.)
  call from_param("mag",mag,1.)

  call to_history("esize",4)

  call sep_close()

  allocate(map(az%n,ax%n))

  call sep_read(map)

  do iz=1,az%n;    z=az%o+(iz-1)*az%d
     do ix=1,ax%n; x=ax%o+(ix-1)*ax%d
        if( sqrt( (xx-x)**2 + (zz-z)**2) < rad) then
           map(iz,ix) = mag
        end if
     end do
  end do

  call sep_write(map)

  call exit (0)
end program SPLASH_




