program MAIN
*             ====
*     Efficient 3-D point probe
*     -------------------------
* ... This code was developed & is CopyLefted by:
*
* Han de Bruijn; Production&Services "A little bit of Physics would be  (===)
* TUD Computing Centre; P.O. Box 354  NO Idleness in Mathematics"(HdB) @-O^O-@
* 2600 AJ  Delft; The Netherlands -- http://pubwww.tudelft.nl/~rcpshdb  #/_\#
* E-mail: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 278 3787. Tel: 2751.   ###
*
      parameter (LIMIT=25,MANY=75)
      integer kon(MANY,4),nok(LIMIT,25)
      real x(LIMIT),y(LIMIT),z(LIMIT)
      real xe(4),ye(4),ze(4),fp(4)
      character cpx*5,cpy*5,cpz*5
      integer merge(MANY),node(3),more(3)
      logical done(MANY)
*
* ... Read coordinates:
      open(1,file='coor3d.dat')
      rewind(1)
      do 10 k=1,LIMIT
      nok(k,1)=0
10    read(1,*,end=11) x(k),y(k),z(k)
11    close(1)
*
* ... Read topology:
      open(1,file='topol3d.dat')
      rewind(1)
      kount=0
      do 20 k=1,MANY
      done(k)=.false.
      merge(k)=0
      kount=kount+1
20    read(1,*,end=21) (kon(k,m),m=1,4)
21    close(1)
      kount=kount-1
*
* ... Dual connectivity:
      do 30 k=1,kount
      do 35 m=1,4
      i=kon(k,m)+1
      nok(i,1)=nok(i,1)+1
      n=nok(i,1)
35    nok(i,n+1)=k
30    continue
*
* ... Coordinates of probe:
      call getarg(1,cpx)
      call getarg(2,cpy)
      call getarg(3,cpz)
      read(cpx,*,err=99) px
      read(cpy,*,err=99) py
      read(cpz,*,err=99) pz
*
      nel=1
* ... Main loop:
100   continue
      write(*,*) nel
* ... Local coordinates:
      do 40 m=1,4
      xe(m)=x(kon(nel,m)+1)
      ye(m)=y(kon(nel,m)+1)
      ze(m)=z(kon(nel,m)+1)
40    continue
* ... X coordinates:
      x1=xe(2)-xe(1)
      x2=xe(3)-xe(1)
      x3=xe(4)-xe(1)
* ... Y coordinates:
      y1=ye(2)-ye(1)
      y2=ye(3)-ye(1)
      y3=ye(4)-ye(1)
* ... Z coordinates:
      z1=ze(2)-ze(1)
      z2=ze(3)-ze(1)
      z3=ze(4)-ze(1)
* ... Probe:
      xp=px-xe(1)
      yp=py-ye(1)
      zp=pz-ze(1)
*
* ... Isoparametric mapping. Thanks to MAPLE (-:
      det=x1*y2*z3-x1*y3*z2-y1*x2*z3+y1*x3*z2+z1*x2*y3-z1*x3*y2
      if(det.eq.0.) stop 'Tetrahedron volume zero'
* ... Function  fp  belongs to node (1); replace in 'det' ?1 by ?p :
      fp(2)=(xp*y2*z3-xp*y3*z2-yp*x2*z3+yp*x3*z2+zp*x2*y3-zp*x3*y2)/det
* ... Function  fp  belongs to node (2); replace in 'det' ?2 by ?p :
      fp(3)=(x1*yp*z3-x1*y3*zp-y1*xp*z3+y1*x3*zp+z1*xp*y3-z1*x3*yp)/det
* ... Function  fp  belongs to node (3); replace in 'det' ?3 by ?p :
      fp(4)=(x1*y2*zp-x1*yp*z2-y1*x2*zp+y1*xp*z2+z1*x2*yp-z1*xp*y2)/det
* ... Last but not least:
      fp(1)=1.-fp(2)-fp(3)-fp(4)
*
* ... Inside or outside:
      fun=min(fp(1),fp(2),fp(3),fp(4))
      if(fun.ge.0.) stop 'OK: tetrahedron found'
      do 60 k=1,4
60    if(fp(k).eq.fun) k0=k
      k1=mod(k0,4)+1
      node(1)=kon(nel,k1)+1
      k2=mod(k1,4)+1
      node(2)=kon(nel,k2)+1
      k3=mod(k2,4)+1
      node(3)=kon(nel,k3)+1
*
* ... Add tetrahedrons to dual:
      do 65 k=1,3
65    more(k)=nok(node(k),1)
      do 55 k=1,3
      knoop=node(k)
      do 50 m=1,more(k)
      n=nok(knoop,m+1)
50    merge(n)=merge(n)+1
55    continue
* ... Find neighbour now:
      know=0
      do 90 m=1,more(1)
      n=nok(node(1),m+1)
90    if(merge(n).eq.3 .and. n.ne.nel) know=n
* ... Substract tetrahedrons in dual:
      do 85 k=1,3
85    more(k)=nok(node(k),1)
      do 75 k=1,3
      knoop=node(k)
      do 70 m=1,more(k)
      n=nok(knoop,m+1)
70    merge(n)=merge(n)-1
75    continue
*
* ... Next tetrahedron found:
      done(nel)=.true.
      nel=know
      if(done(nel)) nel=0
      if(nel.eq.0) then
        write(0,*) 'Previous search failed'
* ...   Take next unused tetrahedron:
        do 45 k=1,kount
        if(done(k)) goto 45
	nel=k
	goto 15
45      continue
        stop 'Point outside mesh'
      endif
15    continue
      goto 100
*
99    stop 'Wrong parameters'
      end