program Splash2_
  use sep

  implicit none
  real,    pointer        :: map(:,:), pck(:,:)
  real                    :: xx,zz,rr,rad
  integer                 :: ix,iz
  real                    ::  x, z

  integer                 :: np,nc,ip
  integer                 :: nz,nx
  real                    :: oz,ox
  real                    :: dz,dx

  integer                 :: kx,kz
  integer                 :: mx,mz

  call sep_init()

  call from_history("n1",np,1)
  call from_history("n2",nc); if(nc/=3) call seperr("nc /=3")

  call from_param("rad",rad,1.)
  call to_history("esize",4)

  call from_param("nz",nz,1 ); call to_history("n1",nz)
  call from_param("oz",oz,0.); call to_history("o1",oz)
  call from_param("dz",dz,1.); call to_history("d1",dz)

  call from_param("nx",nx,1 ); call to_history("n2",nx)
  call from_param("ox",ox,0.); call to_history("o2",ox)
  call from_param("dx",dx,1.); call to_history("d2",dx)

  call sep_close()

  kx=1+nint(rad/dx)
  kz=1+nint(rad/dz)

  allocate(map(nz,nx))
  allocate(pck(np,nc))

  call sep_read(pck)

  do ip=1,np
     zz=pck(ip,1)
     xx=pck(ip,2)
     rr=pck(ip,3)

     mx=1+nint((xx-ox)/dx)
     mz=1+nint((zz-oz)/dz)

     if(mx<kx .or. mx>nx-kx) cycle
     if(mz<kz .or. mz>nz-kz) cycle

     do iz=mz-kz,mz+kz;    z=oz+(iz-1)*dz
        do ix=mx-kx,mx+kx; x=ox+(ix-1)*dx
           if( sqrt( (xx-x)**2 + (zz-z)**2) < rad) then
              map(iz,ix) = rr
           end if
        end do
     end do
  end do

  call sep_write(map)

  call exit (0)
end program SPLASH2_




