      program triplot

      type surface
      integer nsurface
      real*4,dimension(:),pointer::r,s,x,y,z,xn,yn,zn
      end type surface

      integer window(4),red(230),green(230),blue(230)
      integer xchoice
      external xchoice
      real*4 xr(8),yr(8),zr(8),tr(8),psi(8)
      real*4 ur(8),vr(8),wr(8)
      real*4 xc(7),yc(7),zc(7),tc(7)
      real*4 uc(7),vc(7),wc(7)
      real*4,dimension(:,:),allocatable::lsfc,lsfr

      character fnme*256
      real*4,dimension(:),allocatable::x,y,z,u,v,w,t,p,s
      real*4,dimension(:,:),allocatable::lsf
      integer,dimension(:,:),allocatable::icon,iface
      integer,dimension(:),allocatable::octree,leaf,node,ftr,rtf,face
      integer,dimension(:),allocatable::kfixx,kfixy,kfixz,kfixt
      logical,dimension(:),allocatable::influid
      type (surface),dimension(:),allocatable::surf

      data xr/0.,1.,0.,1.,0.,1.,0.,1./
      data yr/0.,0.,1.,1.,0.,0.,1.,1./
      data zr/0.,0.,0.,0.,1.,1.,1.,1./

      window(1)=100
      window(2)=100
      window(3)=700
      window(4)=700
      call xopen (window,char(0))
      red=0
      green=0
      blue=0
      red(1)=180
      green(1)=120
      green(3)=255
      ncol=20
        do i=1,ncol
        red(100+i)=int(255.*float(i)/ncol)
        blue(100+i)=255-int(255.*float(i)/ncol)
        enddo
      call xcmap (red,green,blue)
      call xscale (0,window(3),0,window(4),-0.05,1.05,1.05,-0.05,1)
!      call xscale (0,window(3),0,window(4),-0.01,0.01,0.01,-0.01,1)
!      call xscale (0,window(3),0,window(4),-0.05,.25,.85,.55,1)
!      call xscale (0,window(3),0,window(4),.75,1.05,.85,.55,1)
      call xpen (1)

      open (7,file='triplot.txt',status='old')
      read (7,*) xn,yn,zn,d
      read (7,'(a)') fnme
      nfnme=len_trim(fnme)
      close (7)

      xyzn=sqrt(xn**2+yn**2+zn**2)
      xn=xn/xyzn
      yn=yn/xyzn
      zn=zn/xyzn
      d=d/xyzn

      open (7,file=fnme(1:nfnme),status='old')
      read (7,*) noctree,nnode,nleaves,nface,nlsf
      allocate (x(nnode),y(nnode),z(nnode),u(nnode),v(nnode),w(nnode))
      allocate (t(nnode),s(nnode),p(nnode),lsf(nnode,nlsf),lsfr(8,nlsf),lsfc(7,nlsf))
      allocate (kfixt(nnode),kfixx(nnode),kfixy(nnode),kfixz(nnode))
      allocate (surf(nlsf))
      do k=1,nnode
      read (7,*) x(k),y(k),z(k),u(k),v(k),w(k), &
                 (lsf(k,j),j=1,nlsf),t(k),s(k),p(k),kfixx(k), &
                 kfixy(k),kfixz(k),kfixt(k)
      enddo
      allocate (icon(8,nleaves))
      read (7,*) ((icon(k,i),k=1,8),i=1,nleaves)
      allocate (octree(noctree))
      read (7,*) (octree(k),k=1,noctree)
      allocate (iface(9,nface))
      do i=1,nface
      read (7,*) (iface(k,i),k=1,9)
      enddo
      allocate (node(nnode),leaf(nnode),ftr(nnode))
      allocate (rtf(nnode),influid(nnode))
      do k=1,nnode
      read (7,*) node(k),leaf(k),ftr(k),rtf(k),influid(k)
      enddo
      allocate (face(nface))
      do k=1,nface
      read (7,*) face(k)
      enddo
        do k=1,nlsf
        read (7,*) nsurface
        surf(k)%nsurface=nsurface
        allocate (surf(k)%x(nsurface),surf(k)%y(nsurface),surf(k)%z(nsurface), &
                  surf(k)%xn(nsurface),surf(k)%yn(nsurface),surf(k)%zn(nsurface), &
                  surf(k)%r(nsurface),surf(k)%s(nsurface))
          do i=1,nsurface
          read (7,*) rs,ss,xs,ys,zs,xns,yns,zns
          surf(k)%x(i)=xs
          surf(k)%y(i)=ys
          surf(k)%z(i)=zs
          surf(k)%xn(i)=xns
          surf(k)%yn(i)=yns
          surf(k)%zn(i)=zns
          surf(k)%r(i)=rs
          surf(k)%s(i)=ss
          enddo
        enddo
      close (7)

      do k=1,nlsf
      print*,'Surface ',k,surf(k)%nsurface
      print*,minval(surf(k)%x(:)),maxval(surf(k)%x(:))
      print*,minval(surf(k)%y(:)),maxval(surf(k)%y(:))
      print*,minval(surf(k)%z(:)),maxval(surf(k)%z(:))
      print*,minval(surf(k)%xn(:)),maxval(surf(k)%xn(:))
      print*,minval(surf(k)%yn(:)),maxval(surf(k)%yn(:))
      print*,minval(surf(k)%zn(:)),maxval(surf(k)%zn(:))
      print*,minval(surf(k)%r(:)),maxval(surf(k)%r(:))
      print*,minval(surf(k)%s(:)),maxval(surf(k)%s(:))
      enddo

      fac=maxval(sqrt(u**2+v**2+w**2))*20.

      do k=1,nnode
!      if (abs(z(k)-.8125).lt.1.e-6) print*,x(k),y(k),z(k),lsf(k,1),influid(k),t(k)
      enddo

      do ileave=1,nleaves

        do k=1,8
        xr(k)=x(icon(k,ileave))
        yr(k)=y(icon(k,ileave))
        zr(k)=z(icon(k,ileave))
        tr(k)=t(icon(k,ileave))
        ur(k)=u(icon(k,ileave))
        vr(k)=v(icon(k,ileave))
        wr(k)=w(icon(k,ileave))
          do j=1,nlsf
          lsfr(k,j)=lsf(icon(k,ileave),j)
          enddo
        psi(k)=xn*xr(k)+yn*yr(k)+zn*zr(k)-d
        enddo

      ncross=0
      do k=1,2
      do j=1,2
      do i=1,2
      ijk=i+(j-1)*2+(k-1)*4
        if (i.ne.2) then
        ijkp=(i+1)+(j-1)*2+(k-1)*4
          if (psi(ijk)*psi(ijkp).le.0) then
          ncross=ncross+1
          rat=abs(psi(ijk))/(abs(psi(ijk))+abs(psi(ijkp)))
          xc(ncross)=xr(ijk)+(xr(ijkp)-xr(ijk))*rat
          yc(ncross)=yr(ijk)+(yr(ijkp)-yr(ijk))*rat
          zc(ncross)=zr(ijk)+(zr(ijkp)-zr(ijk))*rat
          tc(ncross)=tr(ijk)+(tr(ijkp)-tr(ijk))*rat
          uc(ncross)=ur(ijk)+(ur(ijkp)-ur(ijk))*rat
          vc(ncross)=vr(ijk)+(vr(ijkp)-vr(ijk))*rat
          wc(ncross)=wr(ijk)+(wr(ijkp)-wr(ijk))*rat
            do l=1,nlsf
            lsfc(ncross,l)=lsfr(ijk,l)+(lsfr(ijkp,l)-lsfr(ijk,l))*rat
            enddo
          endif
        endif
        if (j.ne.2) then
        ijkp=i+j*2+(k-1)*4
          if (psi(ijk)*psi(ijkp).le.0) then
          ncross=ncross+1
          xc(ncross)=xr(ijk)+(xr(ijkp)-xr(ijk))*rat
          yc(ncross)=yr(ijk)+(yr(ijkp)-yr(ijk))*rat
          zc(ncross)=zr(ijk)+(zr(ijkp)-zr(ijk))*rat
          tc(ncross)=tr(ijk)+(tr(ijkp)-tr(ijk))*rat
          uc(ncross)=ur(ijk)+(ur(ijkp)-ur(ijk))*rat
          vc(ncross)=vr(ijk)+(vr(ijkp)-vr(ijk))*rat
          wc(ncross)=wr(ijk)+(wr(ijkp)-wr(ijk))*rat
            do l=1,nlsf
            lsfc(ncross,l)=lsfr(ijk,l)+(lsfr(ijkp,l)-lsfr(ijk,l))*rat
            enddo
          endif
        endif
        if (k.ne.2) then
        ijkp=i+(j-1)*2+k*4
          if (psi(ijk)*psi(ijkp).le.0) then
          ncross=ncross+1
          xc(ncross)=xr(ijk)+(xr(ijkp)-xr(ijk))*rat
          yc(ncross)=yr(ijk)+(yr(ijkp)-yr(ijk))*rat
          zc(ncross)=zr(ijk)+(zr(ijkp)-zr(ijk))*rat
          tc(ncross)=tr(ijk)+(tr(ijkp)-tr(ijk))*rat
          uc(ncross)=ur(ijk)+(ur(ijkp)-ur(ijk))*rat
          vc(ncross)=vr(ijk)+(vr(ijkp)-vr(ijk))*rat
          wc(ncross)=wr(ijk)+(wr(ijkp)-wr(ijk))*rat
            do l=1,nlsf
            lsfc(ncross,l)=lsfr(ijk,l)+(lsfr(ijkp,l)-lsfr(ijk,l))*rat
            enddo
          endif
        endif
      enddo
      enddo
      enddo

        if (ncross.gt.0) then

        call order (xc,yc,zc,tc,uc,vc,wc,lsfc,nlsf,ncross,xn,yn,zn)

        xc(ncross+1)=xc(1)
        yc(ncross+1)=yc(1)
        zc(ncross+1)=zc(1)
        tc(ncross+1)=tc(1)
        uc(ncross+1)=uc(1)
        vc(ncross+1)=vc(1)
        wc(ncross+1)=wc(1)
          do k=1,nlsf
          lsfc(ncross+1,k)=lsfc(1,k)
          enddo

          do k=1,ncol
          val1=float(k-1)/float(ncol)
          if (k.eq.1) val1=val1+1.e-5
          val2=float(k)/float(ncol)
          if (k.eq.ncol) val2=val2+1.e-5
          icol=100+k
!          call cont (xc,yc,zc,tc,ncross,val1,val2,icol,xn,yn,zn)
!          call contl (xc,yc,zc,tc,ncross,val1,1,xn,yn,zn)
          enddo

        call xpen (3)
        call project (xc(ncross),yc(ncross),zc(ncross),xn,yn,zn,xcp,ycp)
        call xplot (xcp,ycp,3)
          do k=1,ncross
          call project (xc(k),yc(k),zc(k),xn,yn,zn,xcp,ycp)
          call xplot (xcp,ycp,2)
          enddo

        call xpen (4)
        rad=1./4000.
        rad=1./400.
          do k=1,8
            if (kfixt(icon(k,ileave)).eq.1) then
            call project (xr(k),yr(k),zr(k),xn,yn,zn,xrp,yrp)
!            call xcircle (xrp,yrp,rad)
            endif
          enddo

        call xpen (5)
        xm=sum(xc(1:ncross))/ncross
        ym=sum(yc(1:ncross))/ncross
        zm=sum(zc(1:ncross))/ncross
        um=sum(uc(1:ncross))/ncross
        vm=sum(vc(1:ncross))/ncross
        wm=sum(wc(1:ncross))/ncross
        call project (xm,ym,zm,xn,yn,zn,xmp,ymp)
        call project (xm+um/fac,ym+vm/fac,zm+wm/fac,xn,yn,zn,xmpp,ympp)
        call arrow (xmp,ymp,xmpp,ympp)

          do k=1,nlsf
          call contl (xc,yc,zc,lsfc(1,k),ncross,0.0001,2,xn,yn,zn)
          enddo

        endif

      enddo

      call xpen (6)
        do k=1,nlsf
          do i=1,surf(k)%nsurface
          call project (surf(k)%x(i),surf(k)%y(i),surf(k)%z(i),xn,yn,zn,xmp,ymp)
!          if (k.eq.2) call xcirclef (xmp,ymp,rad)
!          if (k.eq.2) call xcirclef (surf(k)%r(i),surf(k)%s(i),rad)
          enddo
        enddo

      call xbox (xn,yn,zn)

      call xsave ()
      ichoice=xchoice ('Continue            ',1)
      call xclose ()

      end

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

      subroutine order (x,y,z,t,u,v,w,lsf,nlsf,n,xn,yn,zn)

      real*4,dimension(:),allocatable::xp,yp,zp,tp,up,vp,wp,an
      real*4,dimension(:,:),allocatable::lsfp
      integer,dimension(:),allocatable::indx
      real*4 x(7),y(7),z(7),t(7),u(7),v(7),w(7),lsf(7,nlsf)

      allocate (xp(n),yp(n),zp(n),tp(n),up(n),vp(n),wp(n), &
                lsfp(n,nlsf),an(n),indx(n))

      xp=x(1:n)
      yp=y(1:n)
      zp=z(1:n)
      tp=t(1:n)
      up=u(1:n)
      vp=v(1:n)
      wp=w(1:n)
      lsfp=lsf(1:n,:)

      xc=sum(x(1:n))/n
      yc=sum(y(1:n))/n
      zc=sum(z(1:n))/n
        do i=1,n
        an(i)=atan2(((y(i)-yc)*(z(1)-zc)-(z(i)-zc)*(y(1)-yc))*xn+ &
                    ((z(i)-zc)*(x(1)-xc)-(x(i)-xc)*(z(1)-zc))*yn+ &
                    ((x(i)-xc)*(y(1)-yc)-(y(i)-yc)*(x(1)-xc))*zn, &
                     (x(i)-xc)*(x(1)-xc)+ &
                     (y(i)-yc)*(y(1)-yc)+ &
                     (z(i)-zc)*(z(1)-zc))
        enddo

      call indexx (n,an,indx)

        do i=1,n
        x(i)=xp(indx(i))
        y(i)=yp(indx(i))
        z(i)=zp(indx(i))
        t(i)=tp(indx(i))
        u(i)=up(indx(i))
        v(i)=vp(indx(i))
        w(i)=wp(indx(i))
          do k=1,nlsf
          lsf(i,k)=lsfp(indx(i),k)
          enddo
        enddo

      deallocate (xp,yp,zp,tp,up,vp,wp,lsfp,an,indx)

      return
      end

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

      SUBROUTINE INDEXX(N,ARRIN,INDX)
      DIMENSION ARRIN(N),INDX(N)
      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END

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

      subroutine cont (x,y,z,t,n,v1,v2,ic,xn,yn,zn)

      real x(n+1),y(n+1),z(n+1),t(n+1)
      real xpp(10),ypp(10),zpp(10),xppp(10),yppp(10)

      nplot=0
        do k=1,n
          if ((t(k)-v1)*(t(k)-v2).le.0.) then
          nplot=nplot+1
          xpp(nplot)=x(k)
          ypp(nplot)=y(k)
          zpp(nplot)=z(k)
          endif
        ratio1=(v1-t(k))/(t(k+1)-t(k))
        ratio2=(v2-t(k))/(t(k+1)-t(k))
        ratio=min(ratio1,ratio2)
          if (ratio*(ratio-1.).le.0.) then
          nplot=nplot+1
          xpp(nplot)=x(k)+ratio*(x(k+1)-x(k))
          ypp(nplot)=y(k)+ratio*(y(k+1)-y(k))
          zpp(nplot)=z(k)+ratio*(z(k+1)-z(k))
          endif
        ratio=max(ratio1,ratio2)
          if (ratio*(ratio-1.).le.0.) then
          nplot=nplot+1
          xpp(nplot)=x(k)+ratio*(x(k+1)-x(k))
          ypp(nplot)=y(k)+ratio*(y(k+1)-y(k))
          zpp(nplot)=z(k)+ratio*(z(k+1)-z(k))
          endif
        enddo
        if (nplot.ne.0) then
        call xpen(ic)
          do k=1,nplot
          call project (xpp(k),ypp(k),zpp(k),xn,yn,zn,xppp(k),yppp(k))
          enddo
        call xfill (xppp,yppp,nplot)
        endif

      return
      end

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

      subroutine contl (x,y,z,t,n,v,ic,xn,yn,zn)

      real x(n+1),y(n+1),z(n+1),t(n+1)

      ipen=1
      call xpen (ic)
        do k=1,n
          if ((t(k)-v)*(t(k+1)-v).le.0.) then
          ipen=-ipen
          rat=.5
          if (t(k+1).ne.t(k)) rat=(v-t(k))/(t(k+1)-t(k))
          xp=x(k)+rat*(x(k+1)-x(k))
          yp=y(k)+rat*(y(k+1)-y(k))
          zp=z(k)+rat*(z(k+1)-z(k))
          call project (xp,yp,zp,xn,yn,zn,xpp,ypp)
          if (ipen.gt.0) call xplot (xpp,ypp,2)
          if (ipen.lt.0) call xplot (xpp,ypp,3)
          endif
        enddo

      return
      end

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

      subroutine xbox (xn,yn,zn)

      call xpen (1)

      call project (0.,0.,0.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,3)
      call project (1.,0.,0.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,2)
      call project (1.,1.,0.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,2)
      call project (0.,1.,0.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,2)
      call project (0.,0.,0.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,2)
      call project (0.,0.,1.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,2)
      call project (1.,0.,1.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,2)
      call project (1.,1.,1.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,2)
      call project (0.,1.,1.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,2)
      call project (0.,0.,1.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,2)
      call project (1.,0.,0.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,3)
      call project (1.,0.,1.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,2)
      call project (1.,1.,0.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,3)
      call project (1.,1.,1.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,2)
      call project (0.,1.,0.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,3)
      call project (0.,1.,1.,xn,yn,zn,xp,yp)
      call xplot (xp,yp,2)

      return
      end

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

      subroutine project (x,y,z,xn,yn,zn,xp,yp)

      if (abs(xn).ge.abs(yn) .and. abs(xn).ge.abs(zn)) then
      xp=-yn*x+xn*y
      yp=-zn*x+xn*z
      return

      elseif (abs(yn).ge.abs(zn)) then
      xp=yn*x-xn*y
      yp=-zn*y+yn*z
      return

      else
      xp=zn*x-xn*z
      yp=zn*y-yn*z
      return

      endif

      end

!-----
      subroutine arrow (x1,y1,x2,y2)

      call xplot (x1,y1,3)
      call xplot (x2,y2,2)

      if (y2-y1.eq.0. .and. x2-x1.eq.0.) then
      ang=0.
      else
      ang=atan2(y2-y1,x2-x1)
      endif
      ang1=ang+3.141592654+.5
      ang2=ang+3.141592654-.5
      dl=sqrt((x2-x1)**2+(y2-y1)**2)/3.

      call xplot (x2+dl*cos(ang1),y2+dl*sin(ang1),3)
      call xplot (x2,y2,2)
      call xplot (x2+dl*cos(ang2),y2+dl*sin(ang2),2)

      return
      end

