! This version is updated on Jul 10, 2012. Fixed a bug about when k1 and k2 are ! exactly antiparallel to each other, we rearrage them to be parallel before ! calling the ellispes !********************************************* !Common Variables !********************************************* MODULE Common Implicit None Double Precision::a1, b1, c1, a2, b2 ,c2 ! length of semiaxes Double Precision::l1(3), m1(3), n1(3), l2(3), m2(3), n2(3), d(3), p0(3), p0d(3) ! all vectors END MODULE Common !******************************************* !Subroutine to calculate the distance of closest approach of two ellipsoids !************************************** Subroutine ellipsoids(dist) Use Common Implicit None Double Precision,Parameter::tau=1.61803398 Double Precision::theta1,theta2,theta3,theta4,pi Double Precision::dist,f1,f2,f3,f4,f5,f6,err,h=1.D-05 Double Precision,External::norm pi=atan(1.d0)*4.d0 !normalize the vectors l1=l1/norm(l1) m1=m1/norm(m1) n1=n1/norm(n1) l2=l2/norm(l2) m2=m2/norm(m2) n2=n2/norm(n2) d=d/norm(d) !get p0 call cross(d,l1,p0) If(norm(p0)<1.0d-14) then call cross(d,m1,p0) End If p0=p0/norm(p0) call cross(p0,d,p0d) !initialize my choice of theta theta1=0.d0 theta2=pi theta3=theta2-1.d0/tau*(theta2-theta1) theta4=theta1+1.d0/tau*(theta2-theta1) call plane_ellp(theta1,f1) call plane_ellp(theta3,f3) call plane_ellp(theta4,f4) Do while(abs(theta2-theta1)>1.D-9) !******************case a If((f1<=f3 .and. f3<=f4) .or. (f3<=f1 .and. f1<=f4)) then theta1=theta3 f1=f3 theta3=theta4 f3=f4 theta4=theta1+1.d0/tau*(theta2-theta1) call plane_ellp(theta4,f4) !******************case b ElseIf((f1<=f4 .and. f4<=f3) .or. (f4<=f1 .and. f1<=f3)) then theta2=theta4 theta4=theta3 f4=f3 theta3=theta2-1.d0/tau*(theta2-theta1) call plane_ellp(theta3,f3) !******************case c ElseIf((f3<=f4 .and. f4<=f1) .or. (f4<=f3 .and. f3<=f1)) then call plane_ellp(theta1+h,f5) If(f5>f1) then theta2=theta3 theta3=theta2-1.d0/tau*(theta2-theta1) theta4=theta1+1.d0/tau*(theta2-theta1) call plane_ellp(theta3,f3) call plane_ellp(theta4,f4) Else call plane_ellp(theta1-h,f6) if(f6a22)then kpmp=1.d0/dsqrt(1.d0-e1*k1d**2)*b1/a1*k1d elseif(a11