program MAIN
*             ====
*     Efficient 2-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=100,MANY=187)
      integer kon(MANY,3),nok(LIMIT,20)
      real x(LIMIT),y(LIMIT)
      real xe(3),ye(3)
      character cpx*5,cpy*5
      logical done(MANY),next
*
* ... Read coordinates:
      open(1,file='coor2d.dat')
      rewind(1)
      do 10 k=1,LIMIT
      nok(k,1)=0
10    read(1,*,end=11) x(k),y(k)
11    close(1)
*
* ... Read topology:
      open(1,file='topol2d.dat')
      rewind(1)
      do 20 k=1,MANY
      done(k)=.false.
20    read(1,*,end=21) (kon(k,m),m=1,3)
21    close(1)
*
* ... Dual connectivity:
      do 30 k=1,MANY
      do 35 m=1,3
      i=kon(k,m)+1
      nok(i,1)=nok(i,1)+1
      n=nok(i,1)
35    nok(i,n+1)=k
30    continue
* ... Tested: OK.
*
* ... Coordinates of probe:
      call getarg(1,cpx)
      call getarg(2,cpy)
      read(cpx,*,err=99) xp
      read(cpy,*,err=99) yp
*
* ... Visualization here.
*
      nel=1
* ... Main loop:
100   continue
      if(next) then
        write(*,*) 'Previous search failed'
* ...   Take next unused triangle:
        do 45 k=1,MANY
        if(done(k).eq..true.) goto 45
	nel=k
	goto 15
45      continue
        stop 'Point outside mesh'
      endif
15    continue
* ... Local coordinates:
      do 40 m=1,3
      xe(m)=x(kon(nel,m)+1)
      ye(m)=y(kon(nel,m)+1)
40    continue
      x21=xe(2)-xe(1)
      x31=xe(3)-xe(1)
      y21=ye(2)-ye(1)
      y31=ye(3)-ye(1)
      xp1=xp-xe(1)
      yp1=yp-ye(1)
* ... Isoparametric mapping:
      det=x21*y31-x31*y21
      if(det.eq.0.) stop 'Triangle area zero'
      p2=(xp1*y31-x31*yp1)/det
      p3=(x21*yp1-xp1*y21)/det
      p1=1.-p2-p3
*
* ... Inside or outside:
      fun=min(p1,p2,p3)
      if(fun.ge.0.) stop 'OK: triangle found'
      if(p1.eq.fun) memo=1
      if(p2.eq.fun) memo=2
      if(p3.eq.fun) memo=3
      k1=mod(memo,3)+1
      k2=mod(k1,3)+1
      i1=kon(nel,k1)+1
      i2=kon(nel,k2)+1
*
* ... Find neighbour:
      m1=nok(i1,1)
      m2=nok(i2,1)
      k1=2
      k2=2
      next=.false.
      do 50 maal=1,MANY
      lef1=nok(i1,k1)
      lef2=nok(i2,k2)
      if(lef1.eq.nel.or.lef1.lt.lef2) then
         k1=k1+1
* ...    New start-triangle?
         if(k1.gt.m1+1) next=.true.
         goto 50
      endif
      if(lef2.eq.nel.or.lef2.lt.lef1) then
         k2=k2+1
* ...    New start-triangle?
         if(k2.gt.m2+1) next=.true.
         goto 50
      endif
* ... Next triangle found:
      if(lef1.eq.lef2) goto 51 
50    if(next) goto 51
51    write(*,*) lef1
      done(nel)=.true.
      nel=lef1
      if(done(nel)) next=.true.
      goto 100
*
99    stop 'Wrong parameters'
      end