subroutine Cone( Nb, Xc, Yc, Zc, THETA, PHI, BONDLENGTH, Xn, Yn, Zn, Seed )

implicit none

! This subroutine randomly places Nb beads at an angle THETA(i) from a line
! formed by two beads with coordinates Xc(1:2), Yc(1:2), and Zc(1:2).  
! The second and later beads are placed an angle PHI(i) from the initial 
! bead. The angle PHI(i) is the angle of rotation about the axis formed by
! the two initial beads in Xc, Yc, Zc.  The bead is placed a distance 
! BONDLENGTH(i) from the position of the second existing bead.  The 
! coordinates of the new beads are Xn(i), Yn(i), and Zn(i).

integer, intent(in)								:: Nb
real, dimension(2), intent(in)					:: Xc, Yc, Zc
real, dimension(Nb), intent(in)					:: THETA, PHI
real, dimension(Nb), intent(in)					:: BONDLENGTH
real, dimension(Nb),intent(out)					:: Xn, Yn, Zn
integer, intent(inout)							:: Seed
real, external									:: ran2

! Local Variables

integer											:: i
real											:: h
real											:: a, b, c, d, s
real, dimension(3)								:: T, U, V, W

W(1) = Xc(1) - Xc(2)
W(2) = Yc(1) - Yc(2)
W(3) = Zc(1) - Zc(2)

W = W / sqrt( dot_product( W, W ) )

V(1) = 2.0 * ran2(Seed)	- 1.0
V(2) = 2.0 * ran2(Seed)	- 1.0

do while ( V(1) * V(1) + V(2) * V(2) > 1.0 )
	V(1) = 2.0 * ran2(Seed)	- 1.0
	V(2) = 2.0 * ran2(Seed)	- 1.0
end do

a = W(1) * W(1) + W(2) * W(2) - 1.0
b = - 2.0 * ( V(1) * W(1) * W(3) + V(2) * W(2) * W(3) )	
c = V(2) * V(2) * W(1) * W(1) + V(1) * V(1) * W(2) * W(2) + &
	V(1) * V(1) * W(3) * W(3) + V(2) * V(2)	* W(3) * W(3) - &
	2.0 * V(1) * V(2) * W(1) * W(2) - V(1) * V(1) - V(2) * V(2)

d = b * b - 4.0 * a * c

if( d < 0.0 .AND. d > -1.0e-7 ) then
	
	V(3) = -b / ( 2.0 * a )

else if( ran2(Seed) < 0.5 ) then
	
	V(3) = ( -b + sqrt( d ) ) / ( 2.0 * a	)

else

	V(3) = ( -b - sqrt( d ) ) / ( 2.0 * a	)

end if

V = V / sqrt( dot_product( V, V ) )

do i = 1, Nb

	h = 1.0 / tan( THETA(i) )

	c = cos( PHI(i) )
	s = sin( PHI(i) )

	U(1) = c * V(1) + s * W(3) * V(2) - s * W(2) * V(3)
	U(2) = c * V(2) + s * W(1) * V(3) - s * W(3) * V(1)		
	U(3) = c * V(3) + s * W(2) * V(1) - s * W(1) * V(2)

	T(1) = h * W(1) +  U(1)
	T(2) = h * W(2) +  U(2)
	T(3) = h * W(3) +  U(3)

	T = T / sqrt( dot_product( T, T ) )

	Xn(i) = Xc(2) + BONDLENGTH(i) * T(1)
	Yn(i) = Yc(2) + BONDLENGTH(i) * T(2)
	Zn(i) = Zc(2) + BONDLENGTH(i) * T(3)

end do

return 

end	subroutine Cone