#!/usr/bin/tclsh
#
# This scripts displays the rolling sphere used in the tesselation procedure
# to compute the molecular surface via VMD's graphics interface.
# It is used during debugging and all its parameters are thrown out by an
# error message in CandidateForTesselation.cpp.


proc animate_sphere { molid atomid1 atomid2 atomid3 rad x1 y1 z1 x2 y2 z2 {step 25} } {
	global SphereRadius
	global vector1
	global vector2
	global CircleCenter
	global angle
	global radius
	global M_PI
	set SphereRadius $rad

	# avoid some 'animate dup' bug
	set molid [molinfo $molid get id]
	# check molid
	if {[molinfo $molid get numframes] != 1} {
		error "Can't find $molid or more than 1 frame."
	}
	# step must be integer and greater 2
	if {$step != int($step)} {
		error "step must be integer."
	}
	if {$step < 2} {
		error "step should be greater than 2."
	}
	# make step-1 copies of current frame
	for {set i 1} {$i <$step} {incr i} {
		animate dup frame 0 $molid
	}

	# select both atoms
	incr atomid1 -1
	incr atomid2 -1
	incr atomid3 -1
	set selid1 [atomselect $molid "index $atomid1" frame 0]
	set selid2 [atomselect $molid "index $atomid2" frame 0]
	set selid3 [atomselect $molid "index $atomid3" frame 0]

	# calculate origin of rotation
	set Atom1 [$selid1 get {x y z}]
	puts "Atom1 is $Atom1."
	set Atom2 [$selid2 get {x y z}]
	puts "Atom2 is $Atom2."
	set Atom3 [$selid3 get {x y z}]
	puts "Atom3 is $Atom3."
	set CircleCenter [vecadd [vecscale 0.5 [expr $Atom1]] [vecscale 0.5 [expr $Atom2]]]
	set CircleCenter [vecadd [vecscale 0.5 [expr $Atom1]] [vecscale 0.5 [expr $Atom2]]]
	puts "CircleCenter is at $CircleCenter."
	
	# calculate normal
	set Normaltmp [vecsub [expr $Atom1] [expr $Atom2]]
	set CirclePlaneNormal [vecnorm $Normaltmp]
	puts "CirclePlaneNormal is at $CirclePlaneNormal."

	# calculate radius
	set SphereCenter [vecsub [list $x1 $y1 $z1] $CircleCenter]
	puts "SphereCenter is at $SphereCenter."
	set OtherSphereCenter [vecsub [list $x2 $y2 $z2] $CircleCenter]
	puts "OtherSphereCenter is at $OtherSphereCenter."
	set radius [veclength $SphereCenter]
	puts "Radius of circle is $radius."
	set vector1 [vecnorm $SphereCenter]
	puts "vector1 is $vector1."

	# calculate other normal vector on plane, pointing out of triangle
	set vector2 [veccross $CirclePlaneNormal $vector1]
	set helper [vecsub $CircleCenter [expr $Atom3]]
	if {[vecdot $helper $vector2] < -1e-10} {
		puts "vector2 points inwards, rescaling."
		set vector2 [vecscale -1. $vector2]
	}
	puts "vector2 is $vector2."

	# get angle
	set angle [expr acos([vecdot $SphereCenter $OtherSphereCenter]/[veclength $SphereCenter]/[veclength $OtherSphereCenter])]

	# make angle in right direction (rotating outward in the sense of vector2
	set SKP [vecdot $OtherSphereCenter $vector2]
	puts "SKP with SearchDirection is $SKP."
	if {$SKP < -1e-10} {
		set angle [expr 2.*$M_PI - $angle]
	}
	puts "angle is $angle."

	global vmd_frame
	trace variable vmd_frame([molinfo top]) w update_current_sphere
	display update
	#enumerate_atoms 0
	return
}	

proc update_current_sphere {name index op} { 
	global SphereRadius
	global radius
	global angle
	global vector1
	global vector2
	global CircleCenter
	# draw sphere
	set CurrentAngle [expr ($angle/[expr [molinfo $index get numframes]-1])*[molinfo $index get frame]]
	set center [vecadd [vecscale $radius [vecadd [vecscale [expr cos($CurrentAngle)] $vector1] [vecscale [expr sin($CurrentAngle)] $vector2]]] $CircleCenter]
	#puts "center is now at $center with CurrentAngle is $CurrentAngle and radius $SphereRadius."
	draw material Transparent
	draw color silver
	if {[draw exists 2] != 0} {
		draw replace 2
	}
	draw sphere $center radius [expr $SphereRadius] resolution 30
	return
}

proc animate_sphere_off {} {
	global vmd_frame
	trace vdelete vmd_frame([molinfo top]) w update_current_sphere
	draw delete all
	animate delete beg 1 end 25 skip 0 [molinfo top]
	return
}

#enumerate_atoms 0
#animate_sphere 0 15 6 3. 4.00145 1.4968 1.28282 1.97354 -2.53854 -1.81601