subroutine extend1(p1,p2,pp)

       implicit none
       real,dimension(3)::p1,p2,pp,p3
       real::p1p2,dpp,theta

       p1p2=sqrt((p1(1)-p2(1))**2+(p1(2)-p2(2))**2+(p1(3)-p2(3))**2)
       p1p2=0.5*p1p2
       theta=2.*asin(p1p2)
       p3=p2*cos(theta)
       pp=2*p3-p1
       dpp=sqrt(pp(1)**2+pp(2)**2+pp(3)**2)
       pp=pp/dpp

       end subroutine extend1