--- /dev/null
+ subroutine odlodc(r1,r2,a,b,uu,vv,aa,bb,dd)
+ implicit real*8 (a-h,o-z)
+ dimension r1(3),r2(3),a(3),b(3),x(3),y(3)
+ odl(u,v) = (r1(1)-r2(1))**2+(r1(2)-r2(2))**2+(r1(3)-r2(3))**2
+ & + 2*ar*u - 2*br*v - 2*ab*u*v + aa*u**2 + bb*v**2
+c print *,"r1",(r1(i),i=1,3)
+c print *,"r2",(r2(i),i=1,3)
+c print *,"a",(a(i),i=1,3)
+c print *,"b",(b(i),i=1,3)
+ aa = a(1)**2+a(2)**2+a(3)**2
+ bb = b(1)**2+b(2)**2+b(3)**2
+ ab = a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
+ ar = a(1)*(r1(1)-r2(1))+a(2)*(r1(2)-r2(2))+a(3)*(r1(3)-r2(3))
+ br = b(1)*(r1(1)-r2(1))+b(2)*(r1(2)-r2(2))+b(3)*(r1(3)-r2(3))
+ det = aa*bb-ab**2
+c print *,'aa',aa,' bb',bb,' ab',ab,' ar',ar,' br',br,' det',det
+ uu = (-ar*bb+br*ab)/det
+ vv = (br*aa-ar*ab)/det
+c print *,u,v
+ uu=dmin1(uu,1.0d0)
+ uu=dmax1(uu,0.0d0)
+ vv=dmin1(vv,1.0d0)
+ vv=dmax1(vv,0.0d0)
+ dd1 = odl(uu,vv)
+ dd2 = odl(0.0d0,0.0d0)
+ dd3 = odl(0.0d0,1.0d0)
+ dd4 = odl(1.0d0,0.0d0)
+ dd5 = odl(1.0d0,1.0d0)
+ dd = dsqrt(dmin1(dd1,dd2,dd3,dd4,dd5))
+ if (dd.eq.dd2) then
+ uu=0.0d0
+ vv=0.0d0
+ else if (dd.eq.dd3) then
+ uu=0.0d0
+ vv=1.0d0
+ else if (dd.eq.dd4) then
+ uu=1.0d0
+ vv=0.0d0
+ else if (dd.eq.dd5) then
+ uu=1.0d0
+ vv=1.0d0
+ endif
+c Control check
+c do i=1,3
+c x(i)=r1(i)+u*a(i)
+c y(i)=r2(i)+v*b(i)
+c enddo
+c dd1 = (x(1)-y(1))**2+(x(2)-y(2))**2+(x(3)-y(3))**2
+c dd1 = dsqrt(dd1)
+ aa = dsqrt(aa)
+ bb = dsqrt(bb)
+c write (8,*) uu,vv,dd,dd1
+c print *,dd,dd1
+ return
+ end