source/wham/src/odlodc.f

   1       subroutine odlodc(r1,r2,a,b,uu,vv,aa,bb,dd)
   2       implicit real*8 (a-h,o-z)
   3       dimension r1(3),r2(3),a(3),b(3),x(3),y(3)
   4       odl(u,v) = (r1(1)-r2(1))**2+(r1(2)-r2(2))**2+(r1(3)-r2(3))**2
   5      & + 2*ar*u - 2*br*v - 2*ab*u*v + aa*u**2 + bb*v**2
   6 c      print *,"r1",(r1(i),i=1,3)
   7 c      print *,"r2",(r2(i),i=1,3)
   8 c      print *,"a",(a(i),i=1,3)
   9 c      print *,"b",(b(i),i=1,3)
  10       aa = a(1)**2+a(2)**2+a(3)**2
  11       bb = b(1)**2+b(2)**2+b(3)**2
  12       ab = a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
  13       ar = a(1)*(r1(1)-r2(1))+a(2)*(r1(2)-r2(2))+a(3)*(r1(3)-r2(3))
  14       br = b(1)*(r1(1)-r2(1))+b(2)*(r1(2)-r2(2))+b(3)*(r1(3)-r2(3))
  15       det = aa*bb-ab**2
  16 c      print *,'aa',aa,' bb',bb,' ab',ab,' ar',ar,' br',br,' det',det
  17       uu = (-ar*bb+br*ab)/det
  18       vv = (br*aa-ar*ab)/det
  19 c      print *,u,v
  20       uu=dmin1(uu,1.0d0)
  21       uu=dmax1(uu,0.0d0)
  22       vv=dmin1(vv,1.0d0)
  23       vv=dmax1(vv,0.0d0)
  24       dd1 = odl(uu,vv)
  25       dd2 = odl(0.0d0,0.0d0)
  26       dd3 = odl(0.0d0,1.0d0)
  27       dd4 = odl(1.0d0,0.0d0)
  28       dd5 = odl(1.0d0,1.0d0)
  29       dd = dsqrt(dmin1(dd1,dd2,dd3,dd4,dd5))
  30       if (dd.eq.dd2) then
  31         uu=0.0d0
  32         vv=0.0d0
  33       else if (dd.eq.dd3) then
  34         uu=0.0d0
  35         vv=1.0d0
  36       else if (dd.eq.dd4) then
  37         uu=1.0d0
  38         vv=0.0d0
  39       else if (dd.eq.dd5) then
  40         uu=1.0d0
  41         vv=1.0d0
  42       endif
  43 c Control check
  44 c      do i=1,3
  45 c        x(i)=r1(i)+u*a(i)
  46 c        y(i)=r2(i)+v*b(i)
  47 c      enddo
  48 c      dd1 = (x(1)-y(1))**2+(x(2)-y(2))**2+(x(3)-y(3))**2
  49 c      dd1 = dsqrt(dd1)
  50       aa = dsqrt(aa)
  51       bb = dsqrt(bb)
  52 c      write (8,*) uu,vv,dd,dd1
  53 c      print *,dd,dd1
  54       return
  55       end