program plane
  use sep

  implicit none

  integer           :: n1,n2
  integer           :: nhor,ndown,ihor,idown,index
  real              :: dip,ztop,ddown,sx1,sy1,sx2,sy2,pert
  real              :: dx,dy,dz,tg,sxmin,symin,pi
  real, allocatable :: arrayout(:,:)
  
  call sep_init ()
  call from_par ("dip",dip,0.)
  call from_par ("ztop",ztop,0.)
  call from_par ("sx1",sx1,0.)
  call from_par ("sy1",sy1,0.)
  call from_par ("sx2",sx2,2.)
  call from_par ("sy2",sy2,2.)
  call from_par ("nhor",nhor,2)
  call from_par ("ndown",ndown,2)
  call from_par ("ddown",ddown,1.)
  call from_par ("pert",pert,1.)
  n1=4
  n2=nhor*ndown

  call to_history ("n1",n1)
  call to_history ("n2",n2)
  call sep_close ()

  if(nhor < 2) call seperr("nhor must be at least 2")
  if(ndown < 2) call seperr("ndown must be at least 2")

  !-----------
  pi=3.141592
  allocate(arrayout(n1,n2))

  dx=(sx2-sx1)/(nhor-1)
  dy=(sy2-sy1)/(nhor-1)
  dz=ddown*sin(dip*pi/180.)
  tg=ddown*cos(dip*pi/180.)/sqrt(dx*dx+dy*dy)

  !write(0,*) dx,dy,dz,tg,dip

  do ihor=1,nhor
     do idown=1,ndown
        
        index=(idown-1)*nhor+ihor
        arrayout(1,index)= ztop + (idown-1)*dz                    ! z
        arrayout(2,index)= sx1+(ihor-1)*dx + (idown-1)*tg*abs(dy) ! x
        arrayout(3,index)= sy1+(ihor-1)*dy + (idown-1)*tg*abs(dx) ! y
        arrayout(4,index)= pert                        ! perturbation

     end do
  end do

  !-----------

  call sep_write(arrayout)

  deallocate(arrayout)
  call exit(0)

end program plane

















