# Bob Hanson 1:12 PM 7/2/2008  PISEMA dQa test script

/*


the script will 

1) load ala20.pdb
2) for EACH frame on the random helix, generate the Z-rotated frame for theta 0 to 355, every 5 degrees
3) for EACH of those quaternions qb, calculate dqab = qb / qa0
4) draw a point for dqab.

*/

function doit(ndeg,DRAW_ARCS)

draw delete
draw quaternions

var qRot = quaternion({0 0 1},nDeg) * quaternion({1 0 0},tau)
A = qRot%{0 0 1};
draw aaa scale @axisLength {0 0 0} @A
var n = random(1,360)
qRef = quaternion(A,n) * qRot * q0;
cutoff = 0.05

	var thetaPhi = data("thetaphi").split()
	var n = thetaPhi.size

	p = array()
	aminoresno = array();
	aminoname = array();
	aminoN = array();
	aminonorm = array();
	aminoOK = array();
	var pt = 1
	for (var i = 1; i <= n; i = i + 1)
		var S = thetaPhi[i].split(",")
		var theta = 0.0 + S[1]
		var phi = 0.0 + S[2]
		if ("" + theta == "NaN")
			continue
		endif
		# undo phi, then undo theta
		var b0 = quaternion({0 0 1},phi) * quaternion({0 1 0},theta);
		p[pt] = !b0
		aminoname[pt] = S[3]
		var rn = 0 + aminoname[pt].split("]")[2]
		var atomN = {*.N and resno=rn}
		aminoN[pt] = atomN
		var x = cross((quaternion({0 0 1},1) * p[pt]) %4, (quaternion({0 0 1},10) * p[pt]) %4)
		aminoNorm[pt] = x/x
		
		var x = "draw vA" + pt + " scale 2.1 vector " + atomN.xyz + " " + qRef%4 + " color "
		var c = qRef%4 * aminoNorm[pt]
		c = c * c
		if (c < cutoff)
			x = x + "red #" + aminoname[pt] + " " + c
			aminoOK[pt] = false
			if (SHOW_DETAIL)
				print x
			endif
		else
			x = x + "pink"
			aminoOK[pt] = true
		endif
		var x = script(x)

		if (SHOW_DETAIL)
			print " " + pt + " " + atomN.label() + " " + rn + " " + aminoname[pt] + " " + c
		endif
		for (j = 0; j < 360; j = j + 30)
			draw ID @{"ptbq" + pt+"_" + j} scale 1.5 vector @atomN @{(quaternion({0 0 1},j) * p[pt]) %4} color white
		end for

		draw ID @{"ptbx" + pt} vector @atomN  @{(p[pt]) %-3} color gray;
		draw ID @{"ptby" + pt} vector @atomN  @{(p[pt]) %-4} color orange;
		draw ID @{"ptbz" + pt} vector @atomN  @{(p[pt]) %-5} color purple;

		pt = pt + 1
	end for

# define all dq normals


var Rz = quaternion({0 0 1},10)

norms = array();
var nnorms = 0;
var normj = array()
var nset = 2
var n = p.size

var ni = 0;
var aalast = "x"
for (var j = 1; j <= n; j = j + 1)
	if (aminoname[j] != aalast)
		ni = ni + 1
		aalast = aminoname[j]
	endif
	var qb0 = p[j];
	var dqb0 = qb0 / qRef
	var dqb1 = Rz * dqb0
	var nb = cross((dqb0%-1),(dqb1%-1))
	var nb = nb/nb
	nnorms = nnorms + 1
	norms[nnorms] = nb
	normj[nnorms] = j
	if (SHOW_DETAIL)
		print "norm " + aalast + " " + nb
	endif
end for


if (DRAW_ARCS)
	offset = {0 0 0}
	var Rz = quaternion({0 0 1},1)
	for (var i = 1; i <= nnorms; i = i + 1)
		var qb0 = p[i];
		var p1 = 0
		for (var j = 0; j <= 360; j = j + 5)
			var dqb = ((Rz * j) * qb0) / qRef
			var x = ({0 0 0} + dqb%{0 0 1 0})*10
			if (p1) 
/*				if (x - p1 > -x);x = -x;endif */
				draw ID @{"dqb_"+i+"_" + j} @{offset + p1} @{offset + x} color red
			end if
			p1 = x
		end for
	end for
endif

# cross all normals
print "crossing all norms: " + nnorms

var diff = 0
var ndiff = 0
var avpt = {0 0 0}

for (var i = 1; i <= nnorms; i = i + 1)
	if (!aminoOK[normj[i]])
		continue
	endif
	for (var j = i + 1; j <= nnorms; j = j + 1)
		if (normj[i] == normj[j] || !aminoOK[normj[j]])
			continue
		endif
		var pt = cross(norms[i],norms[j])
		pt = (pt / pt) * (axisLength/2)
		if (pt - $aaa[2] > -pt - $aaa[2])
			pt = -pt
		endif
		diff = diff + angle(pt, {0 0 0}, $aaa[2])
		pt = (pt / pt) * (axisLength/2 + random(0.5))
		avpt = avpt + pt
		ndiff = ndiff + 1
		var c2 = normj[j] - normj[i];
		var c1 = c2
		var c = "[" + ((((1.0*i)/nnorms*256)\1)%256)+"," + (((c1 + 1) * 80)%256) + ","+ (((c2 + 1) * 80)%256) + "]"
		draw ID @{"pt_" + normj[i] + "_" + normj[j]} @pt @{aminoname[normj[i]] + aminoname[normj[j]]} color @c
	end for
end for

var ang = angle(avpt, {0 0 0}, A)

print "for " + ndiff + " crossings, random rot = " + n + "\naverage spherical distance/radians = " + (diff / ndiff * 3.14159 / 180) + "\naxis angle error/deg = " + ang

draw vResult vector {0 0 0} @{avpt / ndiff * 1.1} color yellow

end function

tau = 0
function doTau(ndeg)
	select *
	rotateselected x @{ndeg - tau}
	tau = ndeg
	doit(0, false)
end function


////////////////////////////////////////////////////////////////////////////

load ala20b.pdb;rotate axisAngle {1 0 0 180} 

data "thetaphi"
118.42,-29.74,[ALA]2     ,1,-1,1,-1,*
128.66,76.68,[ALA]3      ,1,-1,1,1,*
116.32,-178.93,[ALA]4    ,1,-1,1,-1,*
111.24,-84.96,[ALA]5     ,1,-1,1,-1,*
124.98,12.78,[ALA]6      ,1,-1,1,1,*
125.26,122.76,[ALA]7     ,1,-1,1,1,*
111.42,-139.2,[ALA]8     ,1,-1,1,-1,*
115.77,-45.31,[ALA]9     ,1,-1,1,-1,*
128.62,58.67,[ALA]10     ,1,-1,1,1,*
118.88,165.41,[ALA]11    ,1,-1,1,1,*
110.36,-100.01,[ALA]12   ,1,-1,1,-1,*
122.55,-4,[ALA]13        ,1,-1,1,-1,*
127.13,105.34,[ALA]14    ,1,-1,1,1,*
112.86,-154.26,[ALA]15   ,1,-1,1,-1,*
113.71,-60.65,[ALA]16    ,1,-1,1,-1,*
127.8,40.77,[ALA]17      ,1,-1,1,1,*
121.46,149.42,[ALA]18    ,1,-1,1,1,*
110.16,-115,[ALA]19      ,1,-1,1,-1,*
119.86,-20.1,[ALA]20     ,1,-1,1,-1,*
end "thetaphi"

select *

axislength = 30.0

set quaternionframe "N"
set drawpicking true
set perspectivedepth off
set drawhover
axes molecular;axes on
set picking center
select *.cb;label %r

SHOWDETAIL = false

//standard helix
q0 = quaternion( 0.049291, -0.788039, -0.609815,  0.068490) 

doTau(60)