Thanks for the reply gummibaer/Norbert. I know no one has need of the code, but I meant in case anyone was willing to take a look at it. Also, thanks for your recommendation on the book. I'll try to get it ASAP.
As for the code:
subroutine cylinderdist_full(bc1,bc2,r1,r2,mu_a,mu_b,dist,ID)
real,intent(in)::bc1(2,3),bc2(2,3),r1,r2
real,intent(out)::mu_a(1,3),mu_b(1,3),dist
integer,intent(out)::ID
real::A(1,3),B(1,3),C(1,3),D(1,3),M1(3,2),M2(3,2)
real:

(1,2),q(1,2),r(1,2),s(1,2),n1(1,2),n2(1,2)
real::d1,d2,u,v
real::left,right
logical judge1, judge2
interface mean
function mean1(v,d)
real,intent(in)::v

)
integer,intent(in)::d
real::mean1
endfunction mean1
function mean2(v,d)
real,intent(in)::v

,

integer,intent(in)::d
real,target::mean21(1,size(v,2)),mean22(size(v,1),1)
real,pointer::mean2

,

endfunction mean2
function mean3(v,d)
real,intent(in)::v

,:,

integer,intent(in)::d
real,target::mean31(1,size(v,2),size(v,3)),mean32(size(v,1),1,size(v,3))&
&,mean33(size(v,1),size(v,2),1)
real,pointer::mean3

,:,

endfunction mean3
function mean4(v,d)
real,intent(in)::v

,:,:,

integer,intent(in)::d
real,target::mean41(1,size(v,2),size(v,3),size(v,4)),mean42(size(v,1),1,size(v,3),size(v,4))&
&,mean43(size(v,1),size(v,2),1,size(v,4)),mean44(size(v,1),size(v,2),size(v,3),1)
real,pointer::mean4

,:,:,

endfunction mean4
endinterface mean
A=reshape((/bc1(1,

/),(/1,3/))
B=reshape((/bc1(2,

/),(/1,3/))
C=reshape((/bc2(1,

/),(/1,3/))
D=reshape((/bc2(2,

/),(/1,3/))
M1=null_3b3(A-B);
r=matmul((C-A),M1);
s=matmul((D-A),M1);
n1=matmul((s-r),reshape((/0,1,-1,0/),(/2,2/)))
n1=n1/norm_bo(n1)
v=dot(n1,s)
d1=norm_bo(v*n1-r)/norm_bo(s-r)*cos_bo((v*n1-r),(s-r))
j=0
if(d1>1.0 .or. d1==1.0) then
d1=1.0
j=j+1
elseif (d1<0.0 .or. d1==0.0) then
d1=0.0
j=j+1
endif
mu_b=C+d1*(D-C)
M2=null_3b3(D-C)
p=matmul((A-D),M2)
q=matmul((B-D),M2)
n2=matmul((q-p),reshape((/0,1,-1,0/),(/2,2/)))
n2=n2/norm_bo(n2)
u=dot(n2,q);
d2=norm_bo(u*n2-p)/norm_bo(q-p)*cos_bo((u*n2-p),(q-p))
if (d2>1.0.or.d2==1.0)then
d2=1.0
j=j+2
elseif (d2<0.0 .or. d2==0.0)then
d2=0.0
j=j+2
endif
mu_a=A+d2*(B-A)
select case(j)
case(1)
d2=dot((mu_b-A),(B-A))/(norm_bo(B-A))**2
if (d2>1.0 .or.d2==1.0)then
d2=1.0
j=j+2
elseif (d2<0.0 .or. d2==0)then
d2=0.0
j=j+2
endif
mu_a=A+d2*(B-A)
case(2)
d1=dot((mu_a-C),(D-C))/(norm_bo(D-C))**2;
if (d1>1.0.or.d1==1)then
d1=1.0
j=j+1
elseif (d1<0.0.or.d1==0)then
d1=0.0
j=j+1
endif
mu_b=C+d1*(D-C)
case(3)
d2=dot((mu_b-A),(B-A))/(norm_bo(B-A))**2
d1=dot((mu_a-C),(D-C))/(norm_bo(D-C))**2
if (d2<1)then
if(d2>0)then
j=j-2
mu_a=A+d2*(B-A)
endif
endif
if (d1<1)then
if(d1>0)then
j=j-1
mu_b=C+d1*(D-C)
endif
endif
endselect
ID=0
if (norm_bo(mu_a-mu_b)<(r1+r2))then
select case(j)
case(0)
ID=1
case(1)
if ((norm_bo(mu_b-mu_a)-r1)<r2*sin_bo((mu_b-mu_a),(C-D)))then
ID=1
endif
case(2)
if ((norm_bo(mu_b-mu_a)-r2)<r1*sin_bo((mu_b-mu_a),(A-B)))then
ID=1
endif
case(3)
if(abs(dot((mu_b-mu_a),(C-D)))/norm_bo(C-D)<r1)then
if(abs(dot((mu_b-mu_a),(A-B)))/norm_bo(A-B)<r2)then
if(norm_bo(mu_a-mu_b)<sqrt(r1**2+r2**2+2*r1*r2*cos_bo((mu_a-mean(bc1,1)),(mu_bmean(
bc2,1)))))then
ID=1
endif
endif
endif
endselect
endif
dist=norm_bo(mu_a-mu_b)
endsubroutine cylinderdist_full
As you can see, I have no idea how I'm supposed to get around to optimizing this. I understand some of the equations, but the programming aspect is not something I have preparation in. I went ahead and downloaded Silverfrost FTN95 so I can try and see how this works.
From here on, I count on the guidance of the wise ones on the forums...