c      
c   Copyright (c) 1992  by  Dwight Meglan
c
      subroutine PlotHiddenLine( BxPt,NmBx,MAX_BOX,VwPt,VwDir,VwWin,error )
      implicit none

      integer*4 NmBx,MAX_BOX
      real*4    BxPt(3,8,MAX_BOX),VwPt(3),VwDir(3),VwWin(4)
      logical*4 error

      real*4    vecA(3),vecB(3),vecC(3),Vec_Dprod,node(3,1600),NrmVw(3),proj
      integer*4 FcPt(4,6),FcNm,num,NmTri,axis,x,y,z,NmPt,i,TriNode(3,600)
      logical*4 vnrm
      
      data FcPt / 1,2,6,5, 2,3,7,6, 3,4,8,7, 4,1,5,8, 4,3,2,1, 7,8,5,6 /
                ! f1       f2       f3       f4       f5       f6
      data x,y,z / 1,2,3 /

c
c     First perform hidden face removal to set up an array containing triangles
c     of the faces of each segment which are visible and an array containing
c     all lines which are part of these faces
c
      error= .false.
      if ( .not. vnrm( VwDir,NrmVw ) ) goto 1000  ! error...
      NmTri= 0
      NmPt= 0
      do num= 1, NmBx
        do i= 1, 8
          do axis= x, z
            node(axis,NmPt+i)= BxPt(axis,i,num)
          end do ! axis
        end do ! i
        do FcNm= 1, 6
          call Vec_Sub( BxPt(x,FcPt(2,FcNm),num),
     &                  BxPt(x,FcPt(1,FcNm),num),
     &                  vecA        )
          call Vec_Sub( BxPt(x,FcPt(3,FcNm),num),
     &                  BxPt(x,FcPt(2,FcNm),num),
     &                  vecB        ) 
          call Vec_Xprod( vecB,vecA,vecC )
          if ( .not. vnrm(vecC,vecC) ) goto 1010 ! error...
          proj= Vec_Dprod( vecC,NrmVw )
c
c         The abs value test is used to check for faces which are almost
c         visible but are calculated as such due to numerical precision
c         error.  A better way to eliminate this is to calculate the
c         projs of all the faces and then select the three most negative
c         projs
c
          if ( (proj .lt. 0.0) .and. (abs(proj) .gt. 0.001) ) then ! face visible
            TriNode(1,NmTri+1)= FcPt(1,FcNm) + NmPt
            TriNode(2,NmTri+1)= FcPt(2,FcNm) + NmPt
            TriNode(3,NmTri+1)= FcPt(3,FcNm) + NmPt
            TriNode(1,NmTri+2)= FcPt(3,FcNm) + NmPt
            TriNode(2,NmTri+2)= FcPt(4,FcNm) + NmPt
            TriNode(3,NmTri+2)= FcPt(1,FcNm) + NmPt
            NmTri= NmTri + 2
          end if ! face visible
        end do ! FcNm
        NmPt= NmPt+8
      end do ! num
c
c     Now call the hidden line removal and drawing routine
c
      call RemoveHiddenLine( NmTri,TriNode,node,VwPt,VwDir,VwWin,error )
      return

1000  write(6,1005)
1005  format(' Error- First vector normalization failed in PlotHiddenLine')
      error= .true.
      return

1010  write(6,1015)
1015  format(' Error- Second vector normalization failed in PlotHiddenLine')
      error= .true.
      return
      end                                                            
c                                            
c   Copyright (c) 1992  by  Dwight Meglan                        
c
      subroutine RemoveHiddenLine( NmTri,TriNode,node,VwPt,VwDir,VwWin,error )
      implicit none

      include 'vogle.fi'

      integer*4 NmTri,TriNode(3,600)
      real*4    node(3,1600),VwPt(3),VwDir(3),VwWin(4)
      logical*4 error

      integer*4 NmLn,TriLn(2,3),TriNm,TriSide,LnNm,pnt,axis,udiv,
     &          vdiv,uwmin,uwmax,vwmin,vwmax,x,y,z,u,v,row,col,TriCnt(40,40),
     &          TriRct(40,40,100),col1,col2,row1,row2,NmTriChk,num,
     &          TriChk(100),NmSeg,SegNm,NodeCon(10,1600),NmCon(1600),
     &          LnNode(2,1600),node1,node2,i,tpnt,lpnt,PtSame
      real*4    uvtri(2,3,600),uvline(2,2,1600),phi,theta,
     &          st,ct,sp,cp,NrmVw(3),utmin(600),utmax(600),vtmin(600),
     &          vecA(3),vecB(3),vecC(3),TriNrm(3,600),ulmin(1600),ulmax(1600),
     &          vlmin(1600),vlmax(1600),udel,vdel,urmin,urange,vrange,
     &          urmax,vrmin,vrmax,SegLst(2,100),s1,s2,s(3),temp,SegMdPt(3),
     &          lambda,Vec_Dprod,xtemp,ytemp,ztemp,vtmax(600),TmpPt(3),
     &          projA,projB
      logical*4 AddTri,SkipLn,cross(3),SkipTri,AddLn,IntsctFnd,vnrm

      data TriLn / 1,2, 2,3, 3,1 /
      data x,y,z,u,v / 1,2,3,1,2 /
      data uwmin,uwmax,vwmin,vwmax / 1,2,3,4 /
c                           
c     Decompose the triangles into lines and store this info.  The arrays
c     NodeCon and NmCon are for bookeeping.  The first records the
c     id number of the nodes which a node is connected to, while the second
c     contains the number of nodes the a node connects to.  This data is
c     used to prevent overlapping lines from being placed in the line array
c     and thereby reducing the number of lines sent through the hidden line
c     removal algorithm.  Due to the special way in which the faces of the
c     boxes which represent the body segments were broken into 2 triangles
c     each, the third line in each triangle crosses the face of the box
c     and therefore should not be displayed.  Thus, the counter for TriSide
c     goes from 1 to 2 instead of 1 to 3.  This eliminates one third of
c     the lines but makes this subroutine not completely general.  For use
c     on a general object change this counter back to 1 to 3.
c
      error= .false.
      do pnt= 1, 3*NmTri
        NmCon(pnt)= 0
      end do
      LnNm= 0
      do TriNm= 1, NmTri
        do TriSide= 1, 2        ! due to the way in which the triangles
          AddLn= .true.       ! are made from the four sided faces of
          node1= TriNode(TriLn(1,TriSide),TriNm) ! the boxes the last
          node2= TriNode(TriLn(2,TriSide),TriNm) ! line is not displayed
          i= 1                                         
          do while( (i .le. NmCon(node1)) .and. (AddLn) )
            if ( node2 .eq. NodeCon(i,node1) ) AddLn= .false.
            i= i + 1
          end do
          if ( AddLn ) then
            NmCon(node1)= NmCon(node1) + 1
            NmCon(node2)= NmCon(node2) + 1
            NodeCon(NmCon(node1),node1)= node2
            NodeCon(NmCon(node2),node2)= node1
            LnNm= LnNm + 1
            LnNode(1,LnNm)= node1
            LnNode(2,LnNm)= node2
          end if ! AddLn
        end do ! TriSide
      end do ! TriNm
      NmLn= LnNm
c                                                             
c     First compute the view angles based upon the view direction.  This
c     is a vector in the direction of viewing and at present is assumed
c     to pass through the origin (0.,0.,0.) of the Global Coordinate System.
c     Once this is done compute the location of the vertices of the triangles
c     and the endpoints of the lines in the picture plane (UV plane).
c
      udiv= 20
      vdiv= 20
      if ( .not. vnrm(VwDir, NrmVw) ) goto 1000 ! error...
      if ( abs(NrmVw(x)) .gt. 0.0001 ) then
        if ( abs(NrmVw(y)) .lt. 0.0001 ) then
          theta= 0.0
        else
          theta= atan2( -NrmVw(y),-NrmVw(x) )
        end if
      else
        if ( -NrmVw(y) .gt. 0.0 ) then
          theta=  3.1415926/2.0
        else
          theta= -3.1415926/2.0
        end if
      endif
      st= sin(theta)
      ct= cos(theta)
      temp= sqrt(NrmVw(x)**2+NrmVw(y)**2)
      if ( temp .gt. 0.0001 ) then
        if ( abs(NrmVw(z)) .lt. 0.0001 ) then
          phi= 0.0
        else
          phi= atan2( -NrmVw(z), temp )
        end if
      else
        if ( -NrmVw(z) .gt. 0.0 ) then
          phi=  3.1415926/2.0
        else
          phi= -3.1415926/2.0
        end if
      end if
      sp= sin(phi)
      cp= cos(phi)
      do TriNm= 1, NmTri
        do pnt= 1, 3
          uvtri(u,pnt,TriNm)= -node(x,TriNode(pnt,TriNm))*st + 
     &                         node(y,TriNode(pnt,TriNm))*ct
          uvtri(v,pnt,TriNm)= -node(x,TriNode(pnt,TriNm))*ct*sp - 
     &                         node(y,TriNode(pnt,TriNm))*st*sp +
     &                         node(z,TriNode(pnt,TriNm))*cp
        end do ! pnt                                        
        utmin(TriNm)= min(uvtri(u,1,TriNm),uvtri(u,2,TriNm),uvtri(u,3,TriNm))
        utmax(TriNm)= max(uvtri(u,1,TriNm),uvtri(u,2,TriNm),uvtri(u,3,TriNm))
        vtmin(TriNm)= min(uvtri(v,1,TriNm),uvtri(v,2,TriNm),uvtri(v,3,TriNm))
        vtmax(TriNm)= max(uvtri(v,1,TriNm),uvtri(v,2,TriNm),uvtri(v,3,TriNm))
        call Vec_Sub( node(x,TriNode(2,TriNm)),node(x,TriNode(1,TriNm)),vecA )
        call Vec_Sub( node(x,TriNode(3,TriNm)),node(x,TriNode(1,TriNm)),vecB )
        call Vec_Xprod( vecB,vecA,vecC )
        if ( .not. vnrm(vecC,TriNrm(x,TriNm)) ) goto 1010 ! error...
      end do ! TriNm
      do LnNm= 1, NmLn
        do pnt= 1, 2
          uvline(u,pnt,LnNm)= -node(x,LnNode(pnt,LnNm))*st + 
     &                         node(y,LnNode(pnt,LnNm))*ct
          uvline(v,pnt,LnNm)= -node(x,LnNode(pnt,LnNm))*ct*sp -
     &                         node(y,LnNode(pnt,LnNm))*st*sp +
     &                         node(z,LnNode(pnt,LnNm))*cp
        end do ! pnt
        ulmin(LnNm)= min(uvline(u,1,LnNm),uvline(u,2,LnNm))
        ulmax(LnNm)= max(uvline(u,1,LnNm),uvline(u,2,LnNm))
        vlmin(LnNm)= min(uvline(v,1,LnNm),uvline(v,2,LnNm))
        vlmax(LnNm)= max(uvline(v,1,LnNm),uvline(v,2,LnNm))
      end do ! LnNm
c
c     Subdivide the picture plane into rectangles and scan triangles to
c     determine which rectangles they lie within
c
      urange= VwWin(uwmax)-VwWin(uwmin)
      vrange= VwWin(vwmax)-VwWin(vwmin)
      udel= urange/udiv  ! the num of divisions here may need to be modified
      vdel= vrange/vdiv  ! depending on the average size of the input triangles
      do row= 1, vdiv
        do col= 1, udiv
          urmin= VwWin(uwmin) + (col-1)*udel
          urmax= urmin + udel
          vrmin= VwWin(vwmin) + (row-1)*vdel
          vrmax= vrmin + vdel
          TriCnt(row,col)= 0
          do TriNm= 1, NmTri
            if ((utmin(TriNm).le.urmax) .and. (utmax(TriNm).ge.urmin) .and.
     &          (vtmin(TriNm).le.vrmax) .and. (vtmax(TriNm).ge.vrmin))then
              TriCnt(row,col)= TriCnt(row,col) + 1   
              TriRct(row,col,TriCnt(row,col))= TriNm
            end if
          end do
        end do ! col
      end do ! row
c
c     Now check each line by first finding which triangles MAY intersect the
c     line by using bounding rectangles.  Next, check for intersection with
c     each of the triangle's lines.  If intersection is found, then check
c     whether line is in front or behind triangle.  If behind clip the line
c     into segments.  Draw segment when finished and move to next line.
c     NOTE!!!!  This algorithm will not work properly if a triangle passes
c     through another.
c
      do LnNm= 1, NmLn
        col1= (ulmin(LnNm)-VwWin(uwmin))/udel - 1
        if ( col1 .lt. 1 ) col1= 1
        col2= (ulmax(LnNm)-VwWin(uwmin))/udel + 1
        if ( col2 .gt. udiv ) col2= udiv
        row1= (vlmin(LnNm)-VwWin(vwmin))/vdel - 1
        if ( row1 .lt. 1 ) row1= 1
        row2= (vlmax(LnNm)-VwWin(vwmin))/vdel + 1
        if ( row2 .gt. vdiv ) row2= vdiv
        NmTriChk= 0
        do col= col1, col2
          do row= row1, row2
            do TriNm= 1, TriCnt(row,col)
              num= 0
              AddTri= .true.
              do while( (num .lt. NmTriChk) .and. (AddTri) )
                num= num + 1
                if ( TriRct(row,col,TriNm) .eq. TriChk(num) ) AddTri=.false.
              end do                         
              if ( AddTri ) then
                NmTriChk= NmTriChk + 1
                TriChk(NmTriChk)= TriRct(row,col,TriNm)
              end if
            end do ! TriNm
          end do ! row
        end do ! col
c
c       NmTriChk now contains how many triangles must be checked against
c       the line for possible intersection and TriChk contains the actual
c       identification number of the triangles.  Now check for intersection
c       of line with each of these rectangles
c
        NmSeg= 1          ! line will start out as a single segment
        SegLst(1,1)= 0.0  ! line is described parametrically between
        SegLst(2,1)= 1.0  ! s= 0 to 1
        num= 0
        SkipLn= .false.
        do while ( (num .lt. NmTriChk) .and. (.not. SkipLn) )
          num= num + 1
          s1= 0.0
          s2= 1.0
          TriNm= TriChk(num)
c
c         If line is one of the sides of the triangle being checked then
c         skip this triangle since it will not affect the visiblity of the
c         line
c
          SkipTri= .false.
          do tpnt= 1, 3
            do lpnt= 1, 2
              if ( TriNode(tpnt,TriNm) .eq. LnNode(lpnt,LnNm) )
     &          SkipTri= .true.
            end do ! lpnt
          end do ! tpnt
          if ( .not. SkipTri ) then
            call InsctTri( uvline(u,1,LnNm),uvtri(u,1,TriNm),
     &                      uvtri(u,2,TriNm),s(1),cross(1)       )
            call InsctTri( uvline(u,1,LnNm),uvtri(u,2,TriNm),
     &                      uvtri(u,3,TriNm),s(2),cross(2)       )
            call InsctTri( uvline(u,1,LnNm),uvtri(u,3,TriNm),
     &                      uvtri(u,1,TriNm),s(3),cross(3)       )
            if ( cross(1) .and. cross(2) ) then ! crosses tri sides 1 and 2
              s1= s(1)
              s2= s(2)
            else if ( cross(1) .and. cross(3) ) then ! crosses tri sides 1 and 3
              s1= s(1)
              s2= s(3)
            else if ( cross(2) .and. cross(3) ) then ! crosses tri sides 2 and 3
              s1= s(2)
              s2= s(3)
            else ! line does not cross the triangle
              SkipTri= .true.
            end if
            if ( ((s1 .lt. 0.0) .and. (s2 .lt. 0.0)) .or. ! line not hidden
     &           ((s1 .gt. 1.0) .and. (s2 .gt. 1.0))      ) SkipTri= .true.  
c
c         Compute the segment of the line covered by the present triangle and
c         update the segment visible array.  First check whether the endpoints
c         of the line lie on either side of the plane defined by the triangle.
c         If so, then see if the line intersects the triangle, and if it
c         does then update the s parameters to reflect the fact that only
c         part of the line is behind the triangle.  If the line does not
c         intersect the triangle then use another method to check whether
c         the line is infront or behind the triangle.  If it is behind then
c         the triangle will hide part (or all) of the line.
c                                                                          
            if ( .not. SkipTri ) then
              if ( s1 .gt. s2 ) then
                temp= s1
                s1= s2
                s2= temp
              end if
              if ( s1 .lt. 0.0 ) s1= -0.01
              if ( s2 .gt. 1.0 ) s2=  1.01
              call Vec_Sub( node(x,LnNode(1,LnNm)),
     &                      node(x,TriNode(1,TriNm)),vecA )
              call Vec_Sub( node(x,LnNode(2,LnNm)),
     &                      node(x,TriNode(1,TriNm)),vecB )
              projA= Vec_Dprod( TriNrm(x,TriNm),vecA )
              projB= Vec_Dprod( TriNrm(x,TriNm),vecB )
              IntsctFnd= .false.
              if ( ( (projA .lt. 0) .and. (projB .gt. 0) ) .or.  ! endpoints on
     &             ( (projA .gt. 0) .and. (projB .lt. 0) )      )! opposite
     &          call LnXTri( node(x,TriNode(1,TriNm)),    ! sides of
     &                       node(x,TriNode(2,TriNm)),    ! triangle
     &                       node(x,TriNode(3,TriNm)),
     &                       node(x,LnNode(1,LnNm)),
     &                       node(x,LnNode(2,LnNm)),
     &                       projA,projB,TriNrm(x,TriNm),NrmVw,
     &                       IntsctFnd,s1,s2     )
              if ( .not. IntsctFnd ) then
                do axis= x, z
                  SegMdPt(axis)= node(axis,LnNode(1,LnNm))+0.5*(s1+s2)*
     &                         ( node(axis,LnNode(2,LnNm))-
     &                           node(axis,LnNode(1,LnNm)) )
                end do ! axis     
                call Vec_Sub( SegMdPt,VwPt,vecA )
                call Vec_Sub( node(x,TriNode(1,TriNm)),VwPt,vecB )
                lambda= Vec_Dprod( TriNrm(x,TriNm),vecB )/
     &                  Vec_Dprod( TriNrm(x,TriNm),vecA )
c                projA= projA*Vec_Dprod( TriNrm(x,TriNm),NrmVw )
c                projB= projB*Vec_Dprod( TriNrm(x,TriNm),NrmVw )
              end if ! check for line one side or other of triangle
c
c             Note that to allow for precision errors, lambda is tested
c             to be less than 1.01 rather 1.0
c
              if ( (lambda .lt. 1.0) .or. IntsctFnd )
c              if ( ((projA .le. 0.0) .and. (projB .le. 0.0)) .or. IntsctFnd )
     &          call UpdateSegLst(SegLst,NmSeg,s1,s2,SkipLn)
            end if ! line does cross present triangle
          end if ! line not a side of the present triangle
        end do ! num .and. .not. SkipLn
c
c       Now draw the segments of the line which have been determined to
c       be visible
c
        if ( .not. SkipLn ) then
          do SegNm= 1, NmSeg
            if ( SegLst(1,SegNm) .ne. SegLst(2,SegNm) ) then ! nonzero
              do axis= x, z
                TmpPt(axis)= node(axis,LnNode(1,LnNm)) + SegLst(1,SegNm)*    
     &                        ( node(axis,LnNode(2,LnNm)) -
     &                          node(axis,LnNode(1,LnNm))   )
              end do ! axis
              call move( TmpPt(x),TmpPt(y),TmpPt(z) )
              do axis= x, z
                TmpPt(axis)= node(axis,LnNode(1,LnNm)) + SegLst(2,SegNm)*
     &                          ( node(axis,LnNode(2,LnNm)) -
     &                            node(axis,LnNode(1,LnNm))   )
              end do ! axis
              call draw( TmpPt(x),TmpPt(y),TmpPt(z) )
            end if
          end do ! SegNm                                
        end if ! .not. SkipLn
      end do ! LnNm
      return

1000  write(6,1005)
1005  format(' Error- First vector normalization failed in RemoveHiddenLine')
      error= .true.
      return

1010  write(6,1015)
1015  format(' Error- Second vector normalization failed in RemoveHiddenLine')
      error= .true.
      return
      end   
c   
c   Copyright (c) 1992  by  Dwight Meglan
c               
      subroutine InsctTri( line,vrtx1,vrtx2,s,cross )
      implicit none

      real*4    line(2,2),vrtx1(2),vrtx2(2),s
      logical*4 cross

      integer*4 u,v
      real*4    denom,numer,t,InsctPt(2),udel,vdel
      data u,v / 1,2 /

      cross= .false.
      s= 0.0   ! parametric intersection point along line
      denom= ( line(u,2)-line(u,1) )*( vrtx1(v)-vrtx2(v) ) -
     &       ( line(v,2)-line(v,1) )*( vrtx1(u)-vrtx2(u) )
      if ( abs(denom) .lt. 0.0001 ) return  ! line and side do not cross
      
      numer= ( line(u,2)-line(u,1) )*( vrtx1(v)-line(v,1) ) -
     &       ( line(v,2)-line(v,1) )*( vrtx1(u)-line(u,1) )
      t= numer/denom  ! parametric intersection point along triangle side
      if ( (t .lt. 0.0) .or. (t .gt. 1.0) ) return ! no intersection
      
      cross= .true.
      InsctPt(u)= vrtx1(u) + t*( vrtx2(u)-vrtx1(u) )
      InsctPt(v)= vrtx1(v) + t*( vrtx2(v)-vrtx1(v) )
      udel= line(u,2) - line(u,1)
      vdel= line(v,2) - line(v,1)
      if ( abs(udel) .gt. abs(vdel) ) then ! line more horizontal and vertical
        if ( abs(udel) .gt. 0.0001 )  s= ( InsctPt(u)-line(u,1) )/udel
      else                                 ! line more vertical than horizontal
        if ( abs(vdel) .gt. 0.0001 )  s= ( InsctPt(v)-line(v,1) )/vdel
      end if
      return
      end
c
c   Copyright (c) 1992  by  Dwight Meglan
c                    
      subroutine LnXTri( pln1,pln2,pln3,linA,linB,projA,projB,PlnNrm,
     &                   VwDir,IntsctFnd,s1,s2 )
      implicit none

      real*4    pln1(3),pln2(3),pln3(3),linA(3),linB(3),projA,projB,PlnNrm(3),
     &          VwDir(3),s1,s2
      logical*4 IntsctFnd

      real*4    matrix(3,3),AB(3,4),params(3),proj,Vec_Dprod,det,temp
      integer*4 i,j
      logical*4 error

      do i= 1, 3
        matrix(i,1)=  pln2(i)
        matrix(i,2)=  pln3(i)
        matrix(i,3)= -linB(i)
      end do
      det= matrix(1,1)*( matrix(2,2)*matrix(3,3)-matrix(2,3)*matrix(3,2) )
     &   - matrix(1,2)*( matrix(2,1)*matrix(3,3)-matrix(2,3)*matrix(3,1) )
     &   + matrix(1,3)*( matrix(2,1)*matrix(3,2)-matrix(2,2)*matrix(3,1) )
      if ( abs(det) .lt. 0.0001 ) then ! line and plane are parallel
        IntsctFnd= .false.
        return
      else        
        do i= 1, 3
          do j= 1, 3
            AB(i,j)= matrix(i,j)
          end do                                                     
          AB(i,4)= linA(i)-pln1(i)
        end do
        call GaussElim( 3,AB,params,error )
        if ( error ) then
          IntsctFnd= .false.
        else
          IntsctFnd= .true.
          do i= 1, 3
            if ( (params(i) .lt. 0.0) .or. (params(i) .gt. 1.0) ) 
     &        IntsctFnd=.false.
          end do     
        end if
        if ( .not. IntsctFnd ) return ! intersection outside triangle
        if ( (params(3) .gt. s1) .and. (params(3) .lt. s2) ) then
          proj= projA*Vec_Dprod( PlnNrm,VwDir ) ! or not on line (s=0to1)
          if ( proj .lt. 0 ) then ! point A of line is infront of triangle
            s1= params(3) ! line hidden from plane intsct to s2
          else
            s2= params(3) ! line hidden from s1 to plane intsct
          end if
          if ( s1 .gt. s2 ) then
            temp= s2
            s2= s1
            s1= temp
          end if
        else
          IntsctFnd= .false.
        end if
      end if
      return
      end
c         
c
c   
        subroutine GaussElim( size,AB,out,error )
        implicit none  

        integer*4 size
        real*4    AB(size,size+1),out(size)
        logical*4 error
                             
        integer*4 i,j,k,row,col,pivot,vnum
        real*4    temp,ratio,value

        error= .false.
        do i= 1, size-1
          pivot= i
          do j= pivot+1, size
            if ( abs(AB(pivot,i)) .lt. abs(AB(j,i)) ) pivot= j
          end do
          if ( abs(AB(pivot,i)) .lt. 0.00001 ) goto 999
          if ( pivot .ne. i ) then
            do col= i, size
              temp= AB(i,col)
              AB(i,col)= AB(pivot,col)
              AB(pivot,col)= temp
            end do
          end if
          do row= i+1,size
            if ( AB(row,i) .ne. 0.0 ) then
              ratio= AB(row,i)/AB(i,i)
              do col= i+1,size
                AB(row,col)= AB(row,col) - ratio*AB(i,col)
              end do
            end if
          end do ! row
        end do ! i
        if ( abs(AB(size,size)) .lt. 0.00001 ) goto 999
        AB(size,size+1)= AB(size,size+1)/AB(size,size)
        do j= 2, size
          vnum= size+1-j
          value= AB(vnum,size+1)
          do k= vnum+1, size
            value= value - AB(vnum,k)*AB(k,size+1)
          end do ! k
          AB(vnum,size+1)= value/AB(vnum,vnum)
        end do
        do i= 1, 3
          out(i)= AB(i,size+1)
        end do
        return
999     error= .true.
        return
        end
c
c   Copyright (c) 1992  by  Dwight Meglan
c
      subroutine UpdateSegLst( SegLst,NmSeg,HideStart,HideEnd,SkipLn )
      implicit none

      integer*4 NmSeg
      real*4    SegLst(2,100),HideStart,HideEnd
      logical*4 SkipLn

      integer*4 start,end,FirstSeg,LastSeg,SegNm,i,j,k
      real*4    temp,StartDif,EndDif
      logical*4 SplitSeg

      integer*4 MIN_LEN

      parameter (MIN_LEN=0.00001)

      data start,end / 1,2 /
c
c     NOTE!!! Segment begin and end points are stored in parametric form
c     in SegLst in ascending order such that SegLst(1-2,1) is the first
c     segment and SegLst(1-2,NmSeg) is the last segment.
c
c     First consider several cases in which the present segments are not modified
c     because the new hidden segment defined by s1 to s2 has already been
c     removed or all the remaining segments are covered by the new hidden
c     segment. NOTE that the difference tests are present to account for
c     numerical precision error.
c
      StartDif= abs(HideStart-SegLst(end,NmSeg))
      EndDif=   abs(HideEnd-SegLst(start,1))
      SkipLn= .false.
      if ( (HideStart .ge. SegLst(end,NmSeg)) .or. (StartDif .lt. MIN_LEN) ) return ! no effect
      if ( (HideEnd .le. SegLst(start,1)) .or. (EndDif .lt. MIN_LEN) ) return ! no effect
      if ( ( (HideStart .le. SegLst(start,1)) .or. (StartDif .lt. MIN_LEN) ) .and.
     &     ( (HideEnd .ge. SegLst(end,NmSeg)) .or. (EndDif .lt. MIN_LEN) )  )  then 
        NmSeg= 0                   ! all visible segments
        SkipLn= .true.            ! are now hidden
        return
      end if
c
c     If new hidden segment has survived this long then each visible line
c     segment must be compared to the new hidden segment.  First Check special
c     case of a single segment.  Check for the beginning being hidden, the
c     end being hidden, or the segment being split into 2 pieces.
c
      if ( NmSeg .eq. 1 ) then
        if ( (HideStart .lt. SegLst(start,1)) .or. 
     &       (abs(HideStart-SegLst(start,1)) .lt. MIN_LEN) ) then
          SegLst(start,1)= HideEnd
        else if ( (HideEnd .gt. SegLst(end,1)) .or.
     &            (abs(HideEnd-SegLst(end,1)) .lt. MIN_LEN) ) then
          SegLst(end,1)= HideStart
        else
          SegLst(end,2)= SegLst(end,1)
          SegLst(end,1)= HideStart
          SegLst(start,2)= HideEnd
          NmSeg= 2
        end if
        return
      end if
c
c     Now consider the case of more than one segment.  First see if the
c     new hidden section starts before or ends after the visible segments
c     and if so then delete the covered segments and clip the new first
c     or last segment if it is partially covered by the hidden section.
c     Lastly, look at the most general case of the hidden section being
c     in the middle of the visible segments.
c
c     First hidden section starting before visible segments
c
      if ( HideStart .le. SegLst(start,1) ) then  ! segments at beginning
        FirstSeg= 1                                 ! will be eliminated
        do while( (SegLst(end,FirstSeg) .le. HideEnd) .and.
     &            (FirstSeg .lt. NmSeg)                      )
          FirstSeg= FirstSeg + 1
        end do
        if ( SegLst(end,FirstSeg)-HideEnd .lt. MIN_LEN ) ! check for near
     &     FirstSeg=FirstSeg+1                           ! miss
        if ( FirstSeg .gt. NmSeg ) then ! near miss occured on last segment
          NmSeg= 0                       ! now all visible segments are
          SkipLn= .true.                ! hidden
          return
        else if ( FirstSeg .eq. NmSeg ) then ! only the last segment is left
          SegLst(end,1)= SegLst(end,NmSeg)
          if ( HideEnd .gt. SegLst(start,NmSeg) ) then
            SegLst(start,1)= HideEnd
          else
            SegLst(start,1)= SegLst(start,NmSeg)
          end if
          NmSeg= 1
        else if ( FirstSeg .eq. 1 ) then ! only first segment clipped
          SegLst(start,1)= HideEnd
        else ! segments previous to FirstSeg must be deleted
          if ( SegLst(start,FirstSeg) .lt. HideEnd ) ! first remaining segment
     &       SegLst(start,FirstSeg)= HideEnd         ! will be clipped
          do SegNm= FirstSeg, NmSeg
            SegLst(start,SegNm-FirstSeg+1)= SegLst(start,SegNm)
            SegLst(end,SegNm-FirstSeg+1)=   SegLst(end,SegNm)
          end do
          NmSeg= NmSeg - FirstSeg + 1
        end if
c
c     Hidden section ending after visible segments
c
      else if ( HideEnd .gt. SegLst(end,NmSeg) ) then   ! segments at end
        LastSeg= NmSeg                                    ! of list will be
        do while( (SegLst(start,LastSeg) .ge. HideStart) .and. ! eliminated
     &            (LastSeg .gt. 1)                                )            
          LastSeg= LastSeg - 1
        end do
        if ( HideStart-SegLst(start,LastSeg) .lt. MIN_LEN ) ! check for near
     &     LastSeg= LastSeg - 1                             ! miss
        if ( LastSeg .lt. 1 ) then ! near miss occured on the first segment
          NmSeg= 0                ! so now all segments are hidden
          SkipLn= .true.
          return
        end if
        if ( SegLst(end,LastSeg) .gt. HideStart ) ! last segment will be
     &     SegLst(end,LastSeg)= HideStart         ! clipped
        NmSeg= LastSeg
c
c     Hidden section lies in the midst of visible segments
c
      else
        SegNm= 0 ! first check for hidden section lying completely within
        SplitSeg= .false.  ! a visible segment
        do while ( (SegNm .lt . NmSeg) .and. (.not. SplitSeg) )
          SegNm= SegNm + 1
          if ( ( (HideStart .ge. SegLst(start,SegNm)) .or.
     &           (abs(HideStart-SegLst(start,SegNm)) .lt. MIN_LEN) ) .and.
     &         ( (HideEnd .le. SegLst(end,SegNm)) .or.
     &           (abs(HideEnd-SegLst(end,SegNm)) .lt. MIN_LEN) )           ) 
     &         SplitSeg=.true.
        end do
        if ( SplitSeg ) then ! break segment into 2 pieces
          StartDif= abs(HideStart - SegLst(start,SegNm))
          EndDif=   abs(SegLst(end,SegNm) - HideEnd)
          if ( (StartDif .lt. MIN_LEN) .and. (EndDif .lt. MIN_LEN) ) then
            do i= SegNm, NmSeg-1 ! delete segment SegNm
              SegLst(start,i)= SegLst(start,i+1)
              SegLst(end,i)=   SegLst(end,i+1)
            end do
            NmSeg= NmSeg - 1
          else if ( StartDif .lt. MIN_LEN ) then ! delete front part of SegNm
            SegLst(start,SegNm)= HideEnd
          else if ( EndDif .lt. MIN_LEN ) then ! delete back part of SegNm
            SegLst(end,SegNm)= HideStart
          else ! split segment
            do i= NmSeg, SegNm, -1 ! make room for a new segment
              SegLst(start,i+1)= SegLst(start,i)
              SegLst(end,i+1)=   SegLst(end,i)
            end do
            SegLst(start,SegNm+1)= HideEnd
            SegLst(end,SegNm)=     HideStart
            NmSeg= NmSeg + 1
          end if
          return
        end if ! split segment into 2 pieces
c
c       Now the most general case... 
c       Find first and last segments either touched by or within the hidden
c       section
c
        FirstSeg= 1                                        ! find first segment
        do while( (SegLst(end,FirstSeg) .le. HideStart) .and.
     &            (FirstSeg .lt. NmSeg)                ) ! hidden by the new
          FirstSeg= FirstSeg + 1                          ! hidden section
        end do
        if ( SegLst(start,FirstSeg) .lt. HideStart ) then
          SegLst(end,FirstSeg)= HideStart  ! first segment end clipped
          FirstSeg= FirstSeg + 1             ! set FirstSeg to next segment
        end if
        LastSeg= NmSeg                                   ! find last segment
        do while( (SegLst(start,LastSeg) .ge. HideEnd) .and.
     &            (LastSeg .gt. 1)                      )  ! hidden by the new
          LastSeg= LastSeg - 1                            ! hidden section
        end do
        if ( SegLst(end,LastSeg) .gt. HideEnd ) then
          SegLst(end,LastSeg)= HideEnd     ! last segment start clipped
          LastSeg= LastSeg - 1               ! set LastSeg to next segment
        end if
        if ( LastSeg .ge. FirstSeg ) then ! delete segments
          do i= 1, LastSeg-FirstSeg+1
            SegLst(start,FirstSeg+i-1)= SegLst(start,LastSeg+i)
            SegLst(end,FirstSeg+i-1)=   SegLst(end,LastSeg+i)
          end do
          NmSeg= NmSeg - (LastSeg-FirstSeg+1)
        end if ! segments need to be deleted
      end if  ! new hidden segment lies at beginning, end or middle existing
c               visible segments
c     Look for any segments which have their beginnings and ends at the
c     same point.  This is an artifact of how the segments are split up
c     Eliminate these zero length segments
c
c      k= NmSeg
c      do i= k, 1, -1
c        if ( abs(SegLst(start,i)-SegLst(end,i)) .lt. 0.00001 ) then
c          do j= i, NmSeg
c            SegLst(start,j)= SegLst(start,j+1)
c            SegLst(end,j)=   SegLst(end,j+1)
c          end do ! j
c          NmSeg= NmSeg - 1
c        end if
c      end do ! k
c      if ( NmSeg .le. 0 ) SkipLn= .true.
      return
      end
