#
# Create a Chimera marker model of an actin filament lattice for fitting into tomography.
# The filaments are parallel to z and are hexagonally packed, and the phase of the twist is
# shown by having each full twist of actin modeled as two different colored segments.
# All units are Angstroms.
#
# Open this script in Chimera to create the marker model.
#
# The marker model can be saved using the Volume Tracer dialog (menu Tools / Volume Data)
# menu entry File / Save Current Marker Set.
#

lattice_spacing = 100	# Distance between parallel actin filaments, Angstroms.

from math import sqrt
plane_spacing = sqrt(3)/2 * lattice_spacing

helical_repeat = 20	# Pitch of actin filament helix.  Each full twist is shown as two colored segments.
                        # These can be aligned with bands where actin twist is in phase in tomogram.

radius = 10		# Radius of cylinder depicting actin filament, Angstroms.

# Dimensions of model.  
width = 1000		# x-axis size, Angstroms.
height = 1000		# y-axis size, Angstroms.
length = 2000		# Length of filaments in model.

nx = int(width / lattice_spacing)	# Number of actin filaments along x axis.
ny = int(height / plane_spacing)        # Number of actin filament planes along y axis.
nz = int(length / helical_repeat)	# Number of full actin twists.

# Colors
color1 = (1,1,.5,1)	# Light yellow. (red, green, blue, opacity) 0-1 scale
color2 = (.5,.5,1,1)	# Light blue.

from VolumePath import Marker_Set, Marker, Link
mset = Marker_Set('actin bundle')
id = 0
for i in range(nx):
    for j in range(ny):
        m2prev = None
        for k in range(nz):
            x = i * lattice_spacing
            if j % 2 == 1:
                x += 0.5 * lattice_spacing	# Every other plane is shifted to pack hexagonally
            y = j * plane_spacing
            z1 = k * helical_repeat
            z2 = z1 + 0.5 * helical_repeat	# Half way through twist
            id += 1
            m1 = Marker(mset, id, (x,y,z1), color1, radius)
            id += 1
            m2 = Marker(mset, id, (x,y,z2), color2, radius)
            Link(m1, m2, color1, radius)
            if m2prev:
                # Connect to previous marker in filament
                Link(m1, m2prev, color2, radius)
            m2prev = m2