c dwr_ab(A,B,R)
c Gaussian Integral between two segments A,B
c the same format as ~/works/dwr/I1.m

       double precision function dwr_ab(A,B,R)
        implicit double precision(a-h,o-z)
        double precision A(3),B(3),R(3),R2(3)
	double precision x1,aa,bb,ab,rr,ar,ar,g1,g2,g3
	complex*8  y1,y2,I
        double precision p1,p2,q1,r1,u1,v1,u2,v2,D,tmp,tmp2
	
	
c R2=cross(A,B)
	R2(1)=A(2)*B(3)-A(3)*B(2)
	R2(2)=A(3)*B(1)-A(1)*B(3)
	R2(3)=A(1)*B(2)-A(2)*B(1)
	x1=dot(R2,R)

	aa=dot(A,A); bb=dot(B,B); ab=dot(A,B); rr=dot(R,R)
	ar=dot(A,R); br=dot(B,R)
	g1=4.0*(aa*bb-ab*ab);g2=bb;g3=-4.0*ab

	p1=4*(aa-ar)/g3;p2=-4*ar/g3
	q1=8*(br*aa-ar*ab)/g1
	r1=4*(aa*rr-ar*ar)/g1
	u1=2*(br-ab)/g2;v1=(rr+aa-2*ar)/g2
	u2=2*br/g2;v2=rr/g2
	D=g3/g1/sqrt(g2)

        tmp=1.0;tmp2=0.0
	y1=I(tmp,p1,q1,r1,u1,v1)-I(tmp2,p1,q1,r1,u1,v1)
	y2=I(tmp,p2,q1,r1,u2,v2)-I(tmp2,p2,q1,r1,u2,v2)

	dwr_ab=x1*D*(y1-y2)

	return
	end