program main

! routine to find trace of laser pointers on generally oriented planar surfaces
!  assumes separate surface specified for each trace (may be identical)

  implicit none
  real*8 :: s, ap1, ap2, x, y, z
  real*8, dimension(3) :: p1, p2 ! 2 points on defining line
  real*8, dimension(5) :: d   !  surface distance to origin (for normal form)
  real*8, dimension(3,5) :: xnorm, r1 ! surface data - normal & point-on-surface
  integer :: i, j, ntrace, nstep, nother
  !
  ! first read in surfaces in point-direct form and convert to normal form
  !  xnorm.(r-r1) = 0.  ! point-direction form
  !  xnorm.r + d  = 0.  ! normal form, d is normal distance from orignin to plane
  !
  read(*,*) ntrace, nstep, nother ! nother are other points that just get passe thru
  do j=1,ntrace
  read(*,*) xnorm(1:3,j), r1(1:3,j) ! unit normal & point on surface
  d(j) = -dot_product(r1(1:3,j), xnorm(1:3,j))
  enddo
  !
  do i=1,nstep
    do j=1,ntrace
      read(*,*) p1, p2
      !
      ap1 = dot_product(xnorm(:,j),p1)
      ap2 = dot_product(xnorm(:,j),p2)
      s=(-d(j)-ap1)/(ap2-ap1)
      x=p1(1) + s*(p2(1)-p1(1))
      y=p1(2) + s*(p2(2)-p1(2))
      z=p1(3) + s*(p2(3)-p1(3))
      write(*,*) x, y, z
    !
    enddo
    !
    do j=1,nother
      read(*,*) p1
      write(*,*) p1
    enddo
    !
  enddo
  !
  stop
end