program OW_mkcig
! 3D-linear Radon transform

  use sep

  implicit none

  real, allocatable, dimension (:,:,:,:,:)            :: out
  integer					:: stat,n1,n2,n3,n4,n5,k1,k2,k3,k4,k5,iy,ihx,ix,idz,izd,ihy
  real                                            :: px,phx,py,phy,o1,o2,o3,o4,o5,d1,d2,d3,d4,d5,gz,fz,rdz,zd,x,hx,y,hy

  call sep_init()
  call sep_begin_prog()

  call from_param("n1",n1,101);call from_param("d1",d1,20.);call from_param("o1",o1,-1000.)  
  call from_param("n2",n2,21);call from_param("d2",d2,20.);call from_param("o2",o2,-200.)
  call from_param("n3",n3,101);call from_param("d3",d3,10.);call from_param("o3",o3,0.)  
  call from_param("n4",n4,51);call from_param("d4",d4,20.);call from_param("o4",o4,-500.)
  call from_param("n5",n5,1);call from_param("d5",d5,20.);call from_param("o5",o5,0.)
  call from_param("px",px,15.);call from_param("phx",phx,0.);
  call from_param("py",py,15.);call from_param("phy",phy,0.);
  call from_param("k1",k1,n1/2);call from_param("k2",k2,n2/2);call from_param("k3",k3,n3/2)  
  call from_param("k4",k4,n4/2);call from_param("k5",k5,n5/2)
  
  call to_history("n1",n1);call to_history("d1",d1);call to_history("o1",o1)  
  call to_history("n2",n2);call to_history("d2",d2);call to_history("o2",o2)
  call to_history("n3",n3);call to_history("d3",d3);call to_history("o3",o3)  
  call to_history("n4",n4);call to_history("d4",d4);call to_history("o4",o4)
  call to_history("n5",n5);call to_history("d5",d5);call to_history("o5",o5)
  call to_history("label1","x")
  call to_history("label2","hx")
  call to_history("label3","z")
  call to_history("label4","y")
  call to_history("label5","hy")
  call to_history("esize",4)
  allocate(out(n1,n2,n3,n4,n5));out=0

  do iy=1,n4
     y=o4+(iy-1)*d4
     do ihx=1,n2
        hx=o2+(ihx-1)*d2
        do ihy=1,n5
           hy=o5+(ihy-1)*d5
           do ix=1,n1
              x=o1+(ix-1)*d1
              zd  = (k3-1)*d3 - tand(phx)*hx - tand(phy)*hy - tand(px)*x - tand(py)*y
              rdz = (zd - o3) / d3
              idz = rdz
              izd = idz + 1
              fz = rdz - idz
              gz = 1.0 - fz
              if ( izd > n3 .or. izd < 1) cycle
              out(ix,ihx,idz,iy,ihy)=out(ix,ihx,idz,iy,ihy)+gz
              out(ix,ihx,izd,iy,ihy)=out(ix,ihx,izd,iy,ihy)+fz
           end do
        end do
     end do
  end do
  call srite("out",out,4*n1*n2*n3*n4*n5)
  deallocate(out)
  call sep_end_prog()
end program 
