program plane 
  use sep
  implicit none

  real   ,pointer,dimension(:,:,:,:)  :: dat
  integer:: ii,jj,n1,n2,t,n3,np,ip
  real :: d2,d3,o2,o3,x,y,dp,op
  real:: d1,o1,mag,a,b,c,d
  logical:: add

  call sep_init()
  call from_param('add',add,.false.)
  call from_param('a',a,1.)
  call from_param('op',op,0.0006)
  call from_param('c',c,0.0006)
  call from_param('d',d,1.5)

  call from_param('d1',d1,0.004)
  call from_param("n1",n1,1024)
  call from_param('o1',o1,0.00)
  call from_param("mag",mag,1.)

  call from_param("n2",n2,128)
  call from_param("d2",d2,10.)
  call from_param("o2",o2,0.)

  call from_param("n3",n3,64)
  call from_param("d3",d3,10.)
  call from_param("o3",o3,0.)

  call from_param("np",np,1)
  call from_Param("dp",dp,0.00003)

  call to_history('esize',4)
  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",np);call to_history("d4",dp);call to_history("o4",op)
  call sep_close()

  allocate(dat(n1,n2,n3,np)); dat=0.

  if (add) call sep_read(dat)

  do ip=1,np
     b=(ip-1.)*dp+op

     do jj=1,n3
        y = (jj-1)*d3+o3

        if (b .le. 0. ) then

           do ii=1,n2
              x = (ii-1)*d2 + o2
              t = 1.5 +(b*x+c*y+d)/(d1*a)

              if (t .le. n1 .and. t .ge. 1) then 
                 dat(t,n2-ii+1,jj,ip)=dat(t,ii,jj,ip)+mag
              end if
           end do

        else

           do ii=1,n2
              x = (ii-1.)*d2 + o2
              t = 1.5 + (b*(x-(n2-1.)*d2+o2 )+c*y+d  )/(d1*a) 
              if (t .le. n1 .and. t .ge. 1) then 
                 dat(t,n2-ii+1,jj,ip)=dat(t,ii,jj,ip)+mag
              end if
           end do

        end if

     end do
     write(0,*)'WORKING ON ',ip,' OF ',np, b

  end do
  call sep_write(dat)
end program plane
  
