!---------------------------------------------------------------------
!
!   calculate the fourth point of the parallelism
!
!   already known the three points coordinates: (x1,y1,z1),(x2,y2,z2)
!                                               (x3,y3,z3)
!
!----------------------------------------------------------------------


          subroutine parallelpt(x1,y1,z1,x2,y2,z2,x3,y3,z3
     &                       ,xtmp,ytmp,ztmp)

          implicit none

           real,intent(inout)::x1,y1,z1,x2,y2,z2,x3,y3,z3
           real,intent(out)::xtmp,ytmp,ztmp
           real::c1,c2,c3,d1,d2,a,b,f,g,h,i,j,k,xtmp1,xtmp2
     &           ,ytmp1,ytmp2,ztmp1,ztmp2,distance

           d1=1-0.5*((x2-x3)**2+(y2-y3)**2+(z2-z3)**2)
           d2=1-0.5*((x2-x1)**2+(y2-y1)**2+(z2-z1)**2)

           a=(d1*z3-d2*z1)/(y1*z3-y3*z1)
           b=(x1*z3-x3*z1)/(y1*z3-y3*z1)
           f=x1*x1+z1*z1
           g=y1*y1+z1*z1
           h=2*x1*y1
           i=-2*x1*d1
           j=-2*y1*d1
           k=d1*d1-z1*z1

           c1=f+g*b*b-h*b
           c2=-2*g*a*b+h*a+i-j*b
           c3=g*a*a+j*a+k

           xtmp1=(-c2+sqrt(c2*c2-4*c1*c3))*0.5/c1
           xtmp2=(-c2-sqrt(c2*c2-4*c1*c3))*0.5/c1

           ytmp1=a-b*xtmp1
           ytmp2=a-b*xtmp2

           ztmp1=-sqrt(1-xtmp1*xtmp1-ytmp1*ytmp1)
           ztmp2=-sqrt(1-xtmp2*xtmp2-ytmp2*ytmp2)
           
            if((xtmp1-x2)**2+(ytmp1-y2)**2+(ztmp1-z2)**2
     &        .gt.(xtmp2-x2)**2+(ytmp2-y2)**2+(ztmp2-z2)**2
     &         )then
               xtmp=xtmp1
               ytmp=ytmp1
               ztmp=ztmp1
            else
               xtmp=xtmp2
               ytmp=ytmp2
               ztmp=ztmp2
            endif

          end subroutine parallelpt