program calc_coher
use puck_mod
use sep3d_struct_mod
use sep
implicit none
type(sep3d) :: image
integer :: ndim,n(3)
real,allocatable,dimension(:,:,:) :: in,dip,coher
real,allocatable,dimension(:,:) :: i2d,d2d,c2d
integer ::nbox(2),i2,i3

call sep_init(SOURCE)
call init_sep3d("in",image,"INPUT")
ndim=sep3d_ndims(image)
nbox=2
call from_param("nbox",nbox,nbox)

if(ndim/=2 .and. ndim/=3) &
  call seperr("can only hand 2 and 3-D data")


n=1
n(1:ndim)=image%n(1:ndim)
allocate(in(n(1),n(2),n(3)))
allocate(coher(n(1),n(2),n(3)))
call sreed("in",in,size(in)*4)

allocate(i2d(n(1),n(2)))
allocate(d2d(n(1),n(2)))
allocate(c2d(n(1),n(2)))
do i3=1,n(3)
  i2d=in(:,:,i3)
  if(.not. puck_2d(i2d,c2d,d2d,nbox))&
    call seperr("trouble running puck2d")
  coher(:,:,i3)=c2d**2
end do
deallocate(i2d,d2d,c2d)
allocate(i2d(n(1),n(3)))
allocate(d2d(n(1),n(3)))
allocate(c2d(n(1),n(3)))
if(n(3)>1) then
    do i2=1,n(2)
    i2d=in(:,i2,:)
    if(.not. puck_2d(i2d,c2d,d2d,nbox))&
      call seperr("trouble running puck2d")
    coher(:,i2,:)=coher(:,i2,:)+c2d**2
  end do
end if
deallocate(i2d,d2d,c2d)
coher=sqrt(coher)
call srite("out",coher,size(coher)*4)

end program
