# camera_travel - Laurence Pearl, November 2003
#     Pythonized by Lieven Buts, November 2003

import cmd
from math import sqrt,acos,sin

def quaternion(view):
	""" 
	Returns a quaternion representation of a view matrix.
	"""

        nxt = [1,2,0]
        q = [0.0,0.0,0.0,1.0]

        tr = view[0]+view[4]+view[8]

        if tr > 0.0 :
                s = sqrt(tr + 1.0)
                qw1 = s / 2.0
                s = 0.5 / s
		return ( (view[5] - view[7]) * s,
                         (view[6] - view[2]) * s,
                         (view[1] - view[3]) * s,
			 qw1                      )
        else :
                i = 0
                if (view[4] > view[0]):     i = 1
                if (view[8] > view[i+3*i]): i = 2
                j = nxt[i]
                k = nxt[j]
                s = sqrt ((view[i+i*3] - (view[j+j*3] + view[k+k*3])) + 1.0)
                q[i] = s * 0.5
                if (s != 0.0): s = 0.5 / s
                q[3] = (view[k+3*j] - view[j+3*k]) * s
                q[j] = (view[j+3*i] + view[i+3*j]) * s
                q[k] = (view[k+3*i] + view[i+3*k]) * s
		return q

def camera_travel(first,nframes=30,sel='(all)',zflag=0,zlevel=2):
	"""
	Generate progressive view matrices to move the camera smoothly
	from the current view to a view defined by a PyMol selection.

        first   - start frame
        nframes - duration
        sel     - atom selection that defines the orientation at the end of the
	          sequence
	zflag   - flag to indicate whether the final view should be 'zoomed'  
  	          or not
	zlevel  - buffer space for zooming final view (as in cmd.zoom)
	"""

	try:
        	first=int(first)
	except:
		print "camera_travel: first frame must be an integer"
		return -1
	
	try:
        	nframes=int(nframes)
	except:
		print "camera_travel: number of frames"

	ff=float(1.0/nframes)

        old_view = cmd.get_view(2)

#       print ( "view : (" + "%8.3f, "*17 + "%8.3f)" ) % (old_view)
#       print "oldtran : %8.3f %8.3f %8.3f" % (old_view[12], old_view[13], old_view[14])

#       do orient operation on selection
        cmd.orient(sel,0)

#       if zoom to selection is required add this into view matrix
        if zflag:
                cmd.zoom(sel,zlevel,0,1)

#       get new view
        new_view = cmd.get_view()

#       capture new zoom/clip parameters
        ozc1 = new_view[11]
        ozc2 = new_view[15]
        ozc3 = new_view[16]

#       calculate shift in zoom/clip parameters
        dzc1 = (ozc1 - old_view[11]) * ff
        dzc2 = (ozc2 - old_view[15]) * ff
        dzc3 = (ozc3 - old_view[16]) * ff

        ozc1 = old_view[11]
        ozc2 = old_view[15]
        ozc3 = old_view[16]

#       capture new translation vector component
        ox = new_view[12]
        oy = new_view[13]
        oz = new_view[14]

#       calculate shift vector
        dx = ox - old_view[12]
        dy = oy - old_view[13]
        dz = oz - old_view[14]

        dx = dx*ff
        dy = dy*ff
        dz = dz*ff

        ox = old_view[12]
        oy = old_view[13]
        oz = old_view[14]


#       capture old and new rotation matrix components in quaternion form

#       m[0][0] = v[0]  m[0][1] = v[1]  m[0][2] = v[2]
#       m[1][0] = v[3]  m[1][1] = v[4]  m[1][2] = v[5]
#       m[2][0] = v[6]  m[2][1] = v[7]  m[2][2] = v[8]

	qx1,qy1,qz1,qw1 = quaternion(old_view)
	qx2,qy2,qz2,qw2 = quaternion(new_view)

#       calc cosine
        cosom = qx1 * qx2 + qy1 * qy2 + qz1 * qz2 + qw1 * qw2

	limit = 0.001
	if cosom>1.0+limit:
		raise ValueError,"Cosine of omega way out of range (positive)"
	elif cosom>1.0:
		print "Warning: cosom corrected from ",cosom,"to",
		cosom = 1.0
		print cosom

	if cosom<-1.0-limit:
		raise ValueError,"Cosine of omega way out of range (negative)"
	elif cosom<-1.0:
		print "Warning: cosom corrected from ",cosom,"to",
		cosom = 1.0
		print cosom

#       adjust signs
        if (cosom < 0.0):
                cosom = -cosom
                to0 = -qx2
                to1 = -qy2
                to2 = -qz2
                to3 = -qw2
        else:
                to0 = qx2
                to1 = qy2
                to2 = qz2
                to3 = qw2

#       calc coefficients
        omega = acos(cosom)
        sinom = sin(omega)
	if sinom==0.0:
		sinom=limit
		print "Warning: sinom corrected!"

#       restore old view
        cmd.set_view( ("%8.3f, " * 17 + "%8.3f") % tuple(old_view) )

#       loop interpolating over nframes generating interpolated quaternion
	for a in range(nframes+1):
                scale0 = sin((1.0 - float(a*ff)) * omega) / sinom
                scale1 = sin(float(a*ff) * omega) / sinom
                rx = scale0 * qx1 + scale1 * to0;
                ry = scale0 * qy1 + scale1 * to1;
                rz = scale0 * qz1 + scale1 * to2;
                rw = scale0 * qw1 + scale1 * to3;

		# convert back to matrix
                x2 = rx + rx
                y2 = ry + ry
                z2 = rz + rz
                xx = rx * x2
                xy = rx * y2
                xz = rx * z2
                yy = ry * y2
                yz = ry * z2
                zz = rz * z2
                wx = rw * x2
                wy = rw * y2
                wz = rw * z2            

                nv0 = 1.0 - (yy + zz)
                nv3 = xy - wz
                nv6 = xz + wy

                nv1 = xy + wz
                nv4 = 1.0 - (xx + zz)
                nv7 = yz - wx

                nv2 = xz - wy
                nv5 = yz + wx
                nv8 = 1.0 - (xx + yy)

		# update translation vector
                ox = ox + dx
                oy = oy + dy
                oz = oz + dz

		# update zoom/clip parameters if required

                if zflag:
                        ozc1 = ozc1 + dzc1
                        ozc2 = ozc2 + dzc2
                        ozc3 = ozc3 + dzc3

                cmd.mdo("%d" % (first), ("set_view (" + "%8.3f, "*17 + "%8.3f)") %
                      (nv0,nv1,nv2,nv3,nv4,nv5,nv6,nv7,nv8,old_view[9],old_view[10],ozc1,
		       ox,oy,oz,ozc2,ozc3,old_view[17]))
                first += 1

        cmd.set_view( ("%8.3f, "*17 + "%8.3f") %  
                   (nv0,nv1,nv2,nv3,nv4,nv5,nv6,nv7,nv8,old_view[9],old_view[10],ozc1,
		    ox,oy,oz,ozc2,ozc3,old_view[17]))


def camera_view_travel(first,nframes=30,new_view=(\     1.000000000,    0.000000000,    0.000000000,\     0.000000000,    1.000000000,    0.000000000,\     0.000000000,    0.000000000,    1.000000000,\     0.000000000,    0.000000000,  -50.000000000,\     0.000000000,    0.000000000,    0.000000000,\    40.000000000,  100.000000000,    0.000000000 )):
	"""
	Generate progressive view matrices to move the camera smoothly
	from the current view to a new view provided as an argument.

        first   - start frame
        nframes - duration
        new_view  - PyMol view matrix that defines the view at the end of the sequence
	"""

	try:
        	first=int(first)
	except:
		print "camera_travel: first frame must be an integer"
		return -1
	
	try:
        	nframes=int(nframes)
	except:
		print "camera_travel: number of frames"

	ff=float(1.0/nframes)

        old_view = cmd.get_view(2)

#       print ( "view : (" + "%8.3f, "*17 + "%8.3f)" ) % (old_view)
#       print "oldtran : %8.3f %8.3f %8.3f" % (old_view[12], old_view[13], old_view[14])

#       capture new zoom/clip parameters
        ozc1 = new_view[11]
        ozc2 = new_view[15]
        ozc3 = new_view[16]

#       calculate shift in zoom/clip parameters
        dzc1 = (ozc1 - old_view[11]) * ff
        dzc2 = (ozc2 - old_view[15]) * ff
        dzc3 = (ozc3 - old_view[16]) * ff

        ozc1 = old_view[11]
        ozc2 = old_view[15]
        ozc3 = old_view[16]

#       capture new translation vector component
        ox = new_view[12]
        oy = new_view[13]
        oz = new_view[14]

#       calculate shift vector
        dx = ox - old_view[12]
        dy = oy - old_view[13]
        dz = oz - old_view[14]

        dx = dx*ff
        dy = dy*ff
        dz = dz*ff

        ox = old_view[12]
        oy = old_view[13]
        oz = old_view[14]


#       capture old and new rotation matrix components in quaternion form

#       m[0][0] = v[0]  m[0][1] = v[1]  m[0][2] = v[2]
#       m[1][0] = v[3]  m[1][1] = v[4]  m[1][2] = v[5]
#       m[2][0] = v[6]  m[2][1] = v[7]  m[2][2] = v[8]

	qx1,qy1,qz1,qw1 = quaternion(old_view)
	qx2,qy2,qz2,qw2 = quaternion(new_view)

#       calc cosine
        cosom = qx1 * qx2 + qy1 * qy2 + qz1 * qz2 + qw1 * qw2

	limit = 0.001
	if cosom>1.0+limit:
		raise ValueError,"Cosine of omega way out of range (positive)"
	elif cosom>1.0:
		print "Warning: cosom corrected from ",cosom,"to",
		cosom = 1.0
		print cosom

	if cosom<-1.0-limit:
		raise ValueError,"Cosine of omega way out of range (negative)"
	elif cosom<-1.0:
		print "Warning: cosom corrected from ",cosom,"to",
		cosom = 1.0
		print cosom

#       adjust signs
        if (cosom < 0.0):
                cosom = -cosom
                to0 = -qx2
                to1 = -qy2
                to2 = -qz2
                to3 = -qw2
        else:
                to0 = qx2
                to1 = qy2
                to2 = qz2
                to3 = qw2

#       calc coefficients
        omega = acos(cosom)
        sinom = sin(omega)
	if sinom==0.0:
		sinom=limit
		print "Warning: sinom corrected!"

#       restore old view
        cmd.set_view( ("%8.3f, " * 17 + "%8.3f") % tuple(old_view) )

#       loop interpolating over nframes generating interpolated quaternion
	for a in range(nframes+1):
                scale0 = sin((1.0 - float(a*ff)) * omega) / sinom
                scale1 = sin(float(a*ff) * omega) / sinom
                rx = scale0 * qx1 + scale1 * to0;
                ry = scale0 * qy1 + scale1 * to1;
                rz = scale0 * qz1 + scale1 * to2;
                rw = scale0 * qw1 + scale1 * to3;

		# convert back to matrix
                x2 = rx + rx
                y2 = ry + ry
                z2 = rz + rz
                xx = rx * x2
                xy = rx * y2
                xz = rx * z2
                yy = ry * y2
                yz = ry * z2
                zz = rz * z2
                wx = rw * x2
                wy = rw * y2
                wz = rw * z2            

                nv0 = 1.0 - (yy + zz)
                nv3 = xy - wz
                nv6 = xz + wy

                nv1 = xy + wz
                nv4 = 1.0 - (xx + zz)
                nv7 = yz - wx

                nv2 = xz - wy
                nv5 = yz + wx
                nv8 = 1.0 - (xx + yy)

		# update translation vector
                ox = ox + dx
                oy = oy + dy
                oz = oz + dz

		# update zoom/clip parameters 

                ozc1 = ozc1 + dzc1
                ozc2 = ozc2 + dzc2
                ozc3 = ozc3 + dzc3

                cmd.mdo("%d" % (first), ("set_view (" + "%8.3f, "*17 + "%8.3f)") %
                      (nv0,nv1,nv2,nv3,nv4,nv5,nv6,nv7,nv8,old_view[9],old_view[10],ozc1,
		       ox,oy,oz,ozc2,ozc3,old_view[17]))
                first += 1

        cmd.set_view( ("%8.3f, "*17 + "%8.3f") %  
                   (nv0,nv1,nv2,nv3,nv4,nv5,nv6,nv7,nv8,old_view[9],old_view[10],ozc1,
		    ox,oy,oz,ozc2,ozc3,old_view[17]))



def fly_sequence(object,first,last):
	"""
	Visit all residues in an object consecutively.
	
	Requires 15 movie frames per residue.

	object -
	first -
	last -
	"""
	
	fr=1
	for ires in range(int(first),int(last)+1):
        	res = "(" + object + (" and resi %d" % (ires,)) + ")"
        	camera_travel(fr,9,res,1,2)
        	fr = fr + 15


# Install new PyMol commands
cmd.extend("seqfly",fly_sequence)
cmd.extend("camtrav",camera_travel)
cmd.extend("camvtrav",camera_view_travel)


base_view = (\     0.627689123,    0.591809154,    0.505735278,\    -0.757547379,    0.314778954,    0.571870685,\     0.179243475,   -0.742075384,    0.645906866,\     0.000000000,   -5.000000000, -293.012847900,\    21.779773712,   60.280826569,   62.276878357,\   246.721801758,  352.805450439,    1.000000000 )

cmd.set_view(base_view)
camera_travel(1,29,'resi 300 and not n;c,n,o,ca',1,4)
camera_view_travel(30,29,base_view)