#!/usr/bin/env python

'''
-------------------
Copyright (c) 2015 Computational Biomechanics (CoBi) Core, Department of
Biomedical Engineering, Cleveland Clinic

Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to permit
persons to whom the Software is furnished to do so, subject to the
following conditions:

The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
USE OR OTHER DEALINGS IN THE SOFTWARE.
-------------------

prescribe_motion_4_FEBio.py

DESCRIPTION:
Python script to read experimental kinetics/kinematics data and apply the
desired curves to an FEBio model.

REQUIREMENTS:
Python (http://www.python.org)
SciPy (http://www.scipy.org)
NumPy (http://www.numpy.org)

WRITTEN BY:
Craig Bennetts, MS
Computational Biomodeling (CoBi) Core
Department of Biomedical Engineering
Lerner Research Institute
Cleveland Clinic
Cleveland, OH
bennect2@ccf.org


LEGEND:
-------
<.> = variable name


INPUT FILES:
------------

<input>.cfg:

	DESCRIPTION:
	This is a user-defined input configuration file to specify all input
	parameters to prescribe anatomically defined kinematics/kinetics to a
	default FEBio model.

	FILE FORMAT (EXAMPLE):
	----------------------
	FEBIO INPUT FILE
	tf_joint_FEBio_v2.feb
	INPUT KINEMATICS FILENAME
	pivot_shift_kinematics_TIBIA.txt
	INPUT KINETICS FILENAME
	pivot_shift_kinetics_TIBIA_2_FEMUR_GLOBAL.txt
	KINEMATICS/KINETICS DURATION (s)
	1.0
	TIME STEPPER LOADCURVE SUB-STEP DURATION FACTORS
	0.5 1.0 0.5
	TIME STEPPER LOADCURVE SUB-STEP NUMBER OF INTERVALS
	1 10 1
	NUMBER OF SOLVER TIME STEPS
	20
	INPUT SPATIAL UNITS
	mm
	INPUT ROTATIONAL UNITS
	Degrees
	INPUT FORCE UNITS
	N
	INPUT TORQUE UNITS
	N*m
	MODEL SPATIAL UNITS
	mm
	MODEL ROTATIONAL UNITS
	Radians
	MODEL FORCE UNITS
	N
	MODEL TORQUE UNITS
	N*mm
	INPUT JOINT SIDE
	Right
	INPUT KINEMATICS/KINETICS BONE
	Femur
	INPUT KINEMATICS/KINETICS BONE SIDE
	Proximal
	MODEL JOINT SIDE
	Right
	MODEL KINEMATICS/KINETICS BONE
	Femur
	MODEL KINEMATICS/KINETICS BONE MAT NUMBER
	1
	MODEL FIXED BONE MAT NUMBER
	2
	MODEL CS TYPE
	Right
	MODEL CS DIMENSION (+x)
	Lateral
	MODEL CS DIMENSION (+y)
	Anterior
	SPECIFIED KINEMATIC DOFs
	'Flexion'
	SPECIFIED KINEMATIC DOFs SCALING FACTORS
	1.0
	SPECIFIED KINETIC DOFs
	'Medial Drawer' 'Posterior Drawer' 'Distraction' 'Varus Torque' 'Internal Rotation Torque'
	SPECIFIED KINETIC DOFs SCALING FACTORS
	1.0 1.0 1.0 1.0 1.0
	KINEMATIC OFFSETS DOFs
	'Medial' 'Anterior' 'Superior' 'Flexion' 'Valgus' 'Internal Rotation' 
	KINEMATIC OFFSETS (mm,radians)
	0.0 0.0 0.0 0.0 0.0 0.0

<input_FEBio_model>.feb:

	DESCRIPTION:
	Default FEBio model file which has pre-existing load curves and associated
	kinematic/kinetic constraints (so it is capable of running prior to
	the modifications applied in this script).

<kinematics>.txt:

	DESCRIPTION:
	Kinematics data file for the user-specified bone/joint of interest to be 
	applied appropriately in the FEBio model, using joint coordinate system (JCS) 
	offsets and the user-specified model coordinate system (CS).

	FORMAT:
	- Comma-separated columns.
	- First line: names for degree of freedom (arbitrary # of columns/DOFs).
	- All other lines: floating point values.

	DEFAULT DEGREE OF FREEDOM (DOF) HEADINGS (first line of file):
	'Medial,Anterior,Superior,Flexion,Valgus,Internal Rotation'

<kinetics>.txt

	DESCRIPTION:
	Kinematics data file for the user-specified bone/joint of interest to be 
	applied appropriately in the FEBio model, using joint coordinate system (JCS) 
	offsets and the user-specified model coordinate system (CS).
	
	FORMAT:
	- Comma-separated columns.
	- First line: names for degree of freedom (arbitrary # of columns/DOFs).
	- All other lines: floating point values.

	DEFAULT DEGREE OF FREEDOM (DOF) HEADINGS (first line of file):
	'Lateral Drawer,Anterior Drawer,Distraction,Extension Torque,Varus Torque,External Rotation Torque'


OUTPUT FILE:
------------

<input_FEBio_model>_MOTION.feb

	DESCRIPTION:
	Modified FEBio model, with specified kinematic/kinetic DOF data curves applied


ASSUMPTIONS:
------------
defined input units (experiment): mm, N, N*m
defined output units (model): mm, N, N*mm
 
'''


import sys
from xml.etree.ElementTree import ElementTree, parse, SubElement, dump, tostring
from numpy import pi


def USAGE(argv):

	print
	print 'USAGE: ' + argv[0] + ' <input.cfg>'
	print
	print '        <input.cfg>   = input configuration file (see script header for format)'
	print


# Function to reformat XML ElementTree
def indent_elements(elem, indent_string, level=0):
    i = "\n" + level*indent_string
    if len(elem):
        if not elem.text or not elem.text.strip():
            elem.text = i + indent_string
        if not elem.tail or not elem.tail.strip():
            elem.tail = i
        for elem in elem:
            indent_elements(elem, indent_string, level+1)
        if not elem.tail or not elem.tail.strip():
            elem.tail = i
    else:
        if level and (not elem.tail or not elem.tail.strip()):
            elem.tail = i


def prescribe_motion_4_FEBio(argv):

	# CHECK SCRIPT USAGE
	if len(argv)!=2:
		USAGE(argv)
		sys.exit(1)

	#--------------------------#
	# INPUT CONFIGURATION FILE #
	#--------------------------#

	# Set configuration filename
	config_filename = argv[1]

	# READ CONFIGURATION FILE
	ConfigFile = open(config_filename,'r')
	ConfigFileLines = ConfigFile.readlines()
	ConfigFile.close()


	# SET INPUT KINEMATICS/KINETICS FILES
	kinematics_filename = ConfigFileLines[ConfigFileLines.index('INPUT KINEMATICS FILENAME\n')+1].strip('\n')
	kinetics_filename = ConfigFileLines[ConfigFileLines.index('INPUT KINETICS FILENAME\n')+1].strip('\n')


	#-----------#
	# SET UNITS #
	#-----------#

	# Set equivalent values for different kinematic units
	# NOTE: if additional units are defined, specify their relative values here
	#       or, set up user-defined file.
	EQ_SPATIAL_UNITS = {}
	EQ_SPATIAL_UNITS['m'] = 1.0
	EQ_SPATIAL_UNITS['mm'] = 1000.0

	# Get valid spatial kinematic units
	DEFINED_SPATIAL_UNITS = []
	for key in EQ_SPATIAL_UNITS.keys():
		DEFINED_SPATIAL_UNITS.append( key )

	EQ_ROTATIONAL_UNITS = {}
	EQ_ROTATIONAL_UNITS['Radians'] = 1.0
	EQ_ROTATIONAL_UNITS['Degrees'] = 180.0/pi
	
	# Get valid rotational kinematic units
	DEFINED_ROTATIONAL_UNITS = []
	for key in EQ_ROTATIONAL_UNITS.keys():
		DEFINED_ROTATIONAL_UNITS.append( key )

	# Set equivalent values for different kinematic units
	# Set kinetic force equivalencies
	EQ_FORCE_UNITS = {}
	EQ_FORCE_UNITS['N'] = 1.0
	
	# Get valid spatial kinematic units
	DEFINED_FORCE_UNITS = []
	for key in EQ_FORCE_UNITS.keys():
		DEFINED_FORCE_UNITS.append( key )

	# Set kinetic torque equivalencies
	EQ_TORQUE_UNITS = {}
	EQ_TORQUE_UNITS['N*m'] = 1.0
	EQ_TORQUE_UNITS['N*mm'] = 1000.0

	# Get valid kinetic torque units
	DEFINED_TORQUE_UNITS = []
	for key in EQ_TORQUE_UNITS.keys():
		DEFINED_TORQUE_UNITS.append( key )


	# SET USER-DEFINED INPUT/MODEL UNITS

	# Set user-defined units for input kinematic/kinetic data
	INPUT_SPATIAL_UNITS = ConfigFileLines[ConfigFileLines.index('INPUT SPATIAL UNITS\n')+1].strip('\n')
	INPUT_ROTATIONAL_UNITS = ConfigFileLines[ConfigFileLines.index('INPUT ROTATIONAL UNITS\n')+1].strip('\n')
	INPUT_FORCE_UNITS = ConfigFileLines[ConfigFileLines.index('INPUT FORCE UNITS\n')+1].strip('\n')
	INPUT_TORQUE_UNITS = ConfigFileLines[ConfigFileLines.index('INPUT TORQUE UNITS\n')+1].strip('\n')

	# Set user-defined units for model
	MODEL_SPATIAL_UNITS = ConfigFileLines[ConfigFileLines.index('MODEL SPATIAL UNITS\n')+1].strip('\n')
	MODEL_ROTATIONAL_UNITS = ConfigFileLines[ConfigFileLines.index('MODEL ROTATIONAL UNITS\n')+1].strip('\n')
	MODEL_FORCE_UNITS = ConfigFileLines[ConfigFileLines.index('MODEL FORCE UNITS\n')+1].strip('\n')
	MODEL_TORQUE_UNITS = ConfigFileLines[ConfigFileLines.index('MODEL TORQUE UNITS\n')+1].strip('\n')


	#---------------------------------#
	# DEFINE JOINT/BONE RELATIONSHIPS #
	#---------------------------------#
	
	# NOTES:
	# - these relationships could be defined by the user an input configuration file, or other file.
	# - capitalized variables are either defined from user input (config) or constants.

	# Set labels for bones in joint
	# NOTE: this type of data-structure is used throughout this script!
	#       this is a two-way dictionary structure for associating opposite terms/DOFs/etc.

	# Set joint relationships	
	OPPOSITE_JOINT_BONE = {}
	OPPOSITE_JOINT_BONE['Femur'] = 'Tibia'
	OPPOSITE_JOINT_BONE['Tibia'] = 'Femur'

	# Set anatomically defined relationships
	OPPOSITE_JOINT_BONE_SIDE = {}
	OPPOSITE_JOINT_BONE_SIDE['Proximal'] = 'Distal'
	OPPOSITE_JOINT_BONE_SIDE['Distal'] = 'Proximal'

	# Cross-reference joint/anatomical reference
	JOINT_BONE_SIDE = {}
	JOINT_BONE_SIDE['Femur'] = 'Proximal'
	JOINT_BONE_SIDE['Tibia'] = 'Distal'

	# Set the side of the body for a bilaterally symmetric joint for the input kinematic/kinetic data
	INPUT_JOINT_SIDE = ConfigFileLines[ConfigFileLines.index('INPUT JOINT SIDE\n')+1].strip('\n')

	# Set the bone with the anatomical CS that the input kinematic/kinetic data describes
	INPUT_MOVING_BONE = ConfigFileLines[ConfigFileLines.index('INPUT KINEMATICS/KINETICS BONE\n')+1].strip('\n')
	INPUT_MOVING_BONE_SIDE = JOINT_BONE_SIDE[ INPUT_MOVING_BONE ]  # set from config file!

	# Set the input reference bone for the joint (based on the input moving bone)
	INPUT_FIXED_BONE = OPPOSITE_JOINT_BONE[ INPUT_MOVING_BONE ]
	INPUT_FIXED_BONE_SIDE = OPPOSITE_JOINT_BONE_SIDE[ INPUT_MOVING_BONE_SIDE ]
	
	# Set the side of the body for a bilaterally symmetric joint for the input kinematic/kinetic data
	MODEL_JOINT_SIDE = ConfigFileLines[ConfigFileLines.index('MODEL JOINT SIDE\n')+1].strip('\n')

	# Set the bone with the anatomical CS that the kinematics/kinetics should be applied to in the model
	MODEL_MOVING_BONE = ConfigFileLines[ConfigFileLines.index('MODEL KINEMATICS/KINETICS BONE\n')+1].strip('\n')
	MODEL_MOVING_BONE_SIDE = JOINT_BONE_SIDE[ MODEL_MOVING_BONE ]

	# Set the model reference bone for the joint (based on the model moving bone)
	MODEL_FIXED_BONE = OPPOSITE_JOINT_BONE[ INPUT_MOVING_BONE ]
	MODEL_FIXED_BONE_SIDE = OPPOSITE_JOINT_BONE_SIDE[ INPUT_MOVING_BONE_SIDE ]
	

	# relate opposite sides of body for bilaterally symmetric joints
	OPPOSITE_JOINT_SIDE = {}
	OPPOSITE_JOINT_SIDE['Right'] = 'Left'
	OPPOSITE_JOINT_SIDE['Left'] = 'Right'

	# model coordinate system (CS).  Assumed values: 'Right', 'Left'
	MODEL_CS_TYPE = ConfigFileLines[ConfigFileLines.index('MODEL CS TYPE\n')+1].strip('\n')


	#----------------------------------------#
	# DEFINE CS/ANATOMICAL DOF RELATIONSHIPS #
	#----------------------------------------#

	# ORIGINAL INPUT DATA FILE FORMAT:
	# Headings (first line) of kinematics.txt file:
	#      'Medial, Anterior, Superior, Flexion, Valgus, Internal Rotation\n'
	# Headings (first line) of kinetics.txt file:
	#      'Lateral Drawer, Anterior Drawer, Distraction, Extension Torque, Varus Torque, External Rotation Torque\n'
	#
	# OPPOSITE DOFs:
	#    - positive(+)/negative(-) [based on the original input file heading format
	#
	# Opposite kinematic DOFs:
	#    - medial/lateral
	#    - anterior/posterior
	#    - superior/inferior
	#    - flexion/extension
	#    - valgus/varus
	#    - internal rotation/external rotation
	#
	# Opposite kinetic DOFs:
	#    - lateral drawer/medial drawer
	#    - anterior drawer/posterior drawer
	#    - distraction/compression
	#    - extension torque/flexion torque
	#    - varus torque/valgus torque
	#    - internal/external
	#
	# NOTES:
	# - could also have the user define this themselves if input files change at a later date.
	# - each line can list opposing pairs.

	# KINEMATIC DOFs

	# SPATIAL DOFs

	SPATIAL_DOFs = []
	# kinematic spatial DOFs (default input kinematic file column headings)
	SPATIAL_DOFs.append('Medial')
	SPATIAL_DOFs.append('Anterior')
	SPATIAL_DOFs.append('Superior')
	# ... and their opposites
	SPATIAL_DOFs.append('Lateral')
	SPATIAL_DOFs.append('Posterior')
	SPATIAL_DOFs.append('Inferior')
	
	# ROTATIONAL DOFs

	ROTATIONAL_DOFs = []
	# kinematic rotational DOFs (default input kinematic file column headings)
	ROTATIONAL_DOFs.append('Flexion')
	ROTATIONAL_DOFs.append('Valgus')
	ROTATIONAL_DOFs.append('Internal Rotation')
	# ... and their opposites
	ROTATIONAL_DOFs.append('Extension')
	ROTATIONAL_DOFs.append('Varus')
	ROTATIONAL_DOFs.append('External Rotation')


	# KINETIC DOFs

	# Force DOFs

	FORCE_DOFs = []
	# kinetic force DOFs (default input kinetic file column headings)
	FORCE_DOFs.append('Lateral Drawer')
	FORCE_DOFs.append('Anterior Drawer')
	FORCE_DOFs.append('Distraction')
	# ... and their opposites
	FORCE_DOFs.append('Medial Drawer')
	FORCE_DOFs.append('Posterior Drawer')
	FORCE_DOFs.append('Compression')

	# Torque DOFs

	TORQUE_DOFs = []
	# kinetic torque DOFs (default input kinetic file column headings)
	TORQUE_DOFs.append('Extension Torque')
	TORQUE_DOFs.append('Varus Torque')
	TORQUE_DOFs.append('External Rotation Torque')
	# ... and their opposites
	TORQUE_DOFs.append('Flexion Torque')
	TORQUE_DOFs.append('Valgus Torque')
	TORQUE_DOFs.append('Internal Rotation Torque')


	# ASSOCIATE OPPOSITE DOFs

	OPPOSITE_DOFs = {}

	# kinematic DOF default keys
	OPPOSITE_DOFs['Medial'] = 'Lateral'
	OPPOSITE_DOFs['Anterior'] = 'Posterior'
	OPPOSITE_DOFs['Superior'] = 'Inferior'
	OPPOSITE_DOFs['Flexion'] = 'Extension'
	OPPOSITE_DOFs['Valgus'] = 'Varus'
	OPPOSITE_DOFs['Internal Rotation'] = 'External Rotation'
	# ... and their opposites
	OPPOSITE_DOFs['Lateral'] = 'Medial'
	OPPOSITE_DOFs['Posterior'] = 'Anterior'
	OPPOSITE_DOFs['Inferior'] = 'Superior'
	OPPOSITE_DOFs['Extension'] = 'Flexion'
	OPPOSITE_DOFs['Varus'] = 'Valgus'
	OPPOSITE_DOFs['External Rotation'] = 'Internal Rotation'

	# kinetic DOF default keys	
	OPPOSITE_DOFs['Lateral Drawer'] = 'Medial Drawer'
	OPPOSITE_DOFs['Anterior Drawer'] = 'Posterior Drawer'
	OPPOSITE_DOFs['Distraction'] = 'Compression'
	OPPOSITE_DOFs['Extension Torque'] = 'Flexion Torque'
	OPPOSITE_DOFs['Varus Torque'] = 'Valgus Torque'
	OPPOSITE_DOFs['External Rotation Torque'] = 'Internal Rotation Torque'
	# ... and their opposites
	OPPOSITE_DOFs['Medial Drawer'] = 'Lateral Drawer'
	OPPOSITE_DOFs['Posterior Drawer'] = 'Anterior Drawer'
	OPPOSITE_DOFs['Compression'] = 'Distraction'
	OPPOSITE_DOFs['Flexion Torque'] = 'Extension Torque'
	OPPOSITE_DOFs['Valgus Torque'] = 'Varus Torque'
	OPPOSITE_DOFs['Internal Rotation Torque'] = 'External Rotation Torque'


	# ASSOCIATE KINEMATIC/KINETIC DOFs

	# ASSUMPTION: input kinematic/kinetic data is for the distal body
	# => WARNING: if input kinematic/kinetic data is for the proximal body, make changes as described below!!!

	MATCHING_DOFs = {}	

	# kinematic DOF default keys 
	MATCHING_DOFs['Medial'] = 'Medial Drawer'
	MATCHING_DOFs['Anterior'] = 'Anterior Drawer'
	MATCHING_DOFs['Superior'] = 'Compression'  # WARNING: change to 'Distraction' if input kinematic/kinetic data is for proximal body
	MATCHING_DOFs['Flexion'] = 'Flexion Torque'
	MATCHING_DOFs['Valgus'] = 'Valgus Torque'
	MATCHING_DOFs['Internal Rotation'] = 'Internal Rotation Torque'
	# ... and their opposites
	MATCHING_DOFs['Lateral'] = 'Lateral Drawer'
	MATCHING_DOFs['Posterior'] = 'Posterior Drawer'
	MATCHING_DOFs['Inferior'] = 'Distraction'  # WARNING: change to 'Compression' if input kinematic/kinetic data is for proximal body
	MATCHING_DOFs['Extension'] = 'Extension Torque'
	MATCHING_DOFs['Varus'] = 'Varus Torque'
	MATCHING_DOFs['External Rotation'] = 'External Rotation Torque'

	# kinetic DOF default keys 
	MATCHING_DOFs['Lateral Drawer'] = 'Lateral'
	MATCHING_DOFs['Anterior Drawer'] = 'Anterior'
	MATCHING_DOFs['Distraction'] = 'Inferior'  # WARNING: change to 'Superior' if input kinematic/kinetic data is for proximal body
	MATCHING_DOFs['Extension Torque'] = 'Extension'
	MATCHING_DOFs['Varus Torque'] = 'Varus'
	MATCHING_DOFs['External Rotation Torque'] = 'External Rotation'
	# ... and their opposites
	MATCHING_DOFs['Medial Drawer'] = 'Medial'
	MATCHING_DOFs['Posterior Drawer'] = 'Posterior'
	MATCHING_DOFs['Compression'] = 'Superior'  # WARNING: change to 'Inferior' if input kinematic/kinetic data is for proximal body
	MATCHING_DOFs['Flexion Torque'] = 'Flexion'
	MATCHING_DOFs['Valgus Torque'] = 'Valgus'
	MATCHING_DOFs['Internal Rotation Torque'] = 'Internal Rotation'


	#---------------------------------#
	# ASSOCIATE INPUT & MODEL CS/DOFs #
	#---------------------------------#

	# TODO: convert this to build model CS/DOF relationships from user input!

	# SET MODEL CS/DOFs & DOFs/CS
	# NOTE: these are the CS/anatomical dimensions that we want to get the user-specified input DOFs to
	
	# associate 'positive CS dimensions' and 'kinematic/kinetic anatomical DOFs', which depends on:
	# - CS type: ['Right','Left']
	# - joint side: ['Right','Left']
	# - moving bone side: ['Proximal','Distal']

	MODEL_CS_DOFs = {}  # positive CS dimensions are dictionary keys, anatomical DOFs are values
	MODEL_DOFs_CS = {}  # anatomical DOFs are dictionary keys, positive CS dimensions are values


	if (MODEL_JOINT_SIDE=='Right' and MODEL_CS_TYPE=='Right') or \
	   (MODEL_JOINT_SIDE=='Left' and MODEL_CS_TYPE=='Left'):

		# associate positive CS->DOF labels
		if MODEL_MOVING_BONE_SIDE=='Distal':
			MODEL_CS_DOFs['x'] = ['Lateral','Lateral Drawer','Extension','Extension Torque']			
			MODEL_CS_DOFs['y'] = ['Anterior','Anterior Drawer','Varus','Varus Torque']
			MODEL_CS_DOFs['z'] = ['Superior','Compression','Internal Rotation','Internal Rotation Torque']
		elif MODEL_MOVING_BONE_SIDE=='Proximal':
			MODEL_CS_DOFs['x'] = ['Lateral','Lateral Drawer','Flexion','Flexion Torque']
			MODEL_CS_DOFs['y'] = ['Anterior','Anterior Drawer','Valgus','Valgus Torque']
			MODEL_CS_DOFs['z'] = ['Superior','Distraction','Internal Rotation','Internal Rotation Torque']
		#else:?

		# ... and their opposites (DOF->CS)
		for DIM in MODEL_CS_DOFs.keys():
			for DOF in MODEL_CS_DOFs[DIM]:
				MODEL_DOFs_CS[DOF] = DIM
					

	elif (MODEL_JOINT_SIDE=='Right' and MODEL_CS_TYPE=='Left') or \
	     (MODEL_JOINT_SIDE=='Left' and MODEL_CS_TYPE=='Right'):

		# ASSUMPTION: positive rotations in a left hand (LH) coordinate system (CS)			
		#             are applied by the left hand rotation rule.
		
		# associate 'positive CS->DOF labels
		if MODEL_MOVING_BONE_SIDE=='Distal':
			MODEL_CS_DOFs['x'] = ['Lateral','Lateral Drawer','Flexion','Flexion Torque']
			MODEL_CS_DOFs['y'] = ['Anterior','Anterior Drawer','Valgus','Valgus Torque']
			MODEL_CS_DOFs['z'] = ['Superior','Compression','External Rotation','External Rotation Torque']
		elif MODEL_MOVING_BONE_SIDE=='Proximal':
			MODEL_CS_DOFs['x'] = ['Lateral','Lateral Drawer','Extension','Extension Torque']			
			MODEL_CS_DOFs['y'] = ['Anterior','Anterior Drawer','Varus','Varus Torque']
			MODEL_CS_DOFs['z'] = ['Superior','Distraction','External Rotation','External Rotation Torque']
		#else:?

		# ... and their opposites (DOF->CS)
		for DIM in MODEL_CS_DOFs.keys():
			for DOF in MODEL_CS_DOFs[DIM]:
				MODEL_DOFs_CS[DOF] = DIM

	# else:


	#================================#
	# INPUT KINEMATICS/KINETICS DATA #
	#================================#

	# set delimiting string to parse input text files
	# NOTE: this could also be user-defined in config file if desired.
	DELIMITER = ','


	#-----------------------#
	# INPUT KINEMATICS DATA #
	#-----------------------#

	# Set the desired kinematic Degrees Of Freedom (DOF) to apply to the FEBio model
		
	line = ConfigFileLines[ConfigFileLines.index('SPECIFIED KINEMATIC DOFs\n')+1].strip('\n')
	line_list = line.split('\'')
	SPECIFIED_KINEMATIC_DOFs = filter(lambda string: string.strip(), line_list)

	#print SPECIFIED_KINEMATIC_DOFs

	# Open kinematics input file
	InputFile = open(kinematics_filename,'r')

	# Assumption: first line defines kinematic DOFs, and the order of the curve data
	line = InputFile.readline().strip('\r\n')  # works for all OS text file carriage returns
	INPUT_KINEMATIC_DOFs = line.split(DELIMITER)  # may be useful later

	#print INPUT_KINEMATIC_DOFs

	# Set the index order of the kinematic DOF keys (i.e. order of columns in input kinematics file)
	input_kinematic_DOFs_index = {}  # keys: INPUT_KINEMATIC_DOFs
	index = 0
	for DOF in INPUT_KINEMATIC_DOFs:
		input_kinematic_DOFs_index[DOF] = index
		index += 1

	#print input_kinematic_DOFs_index


	# NOTE: user-defined kinematic DOFs in configuration file must match the input kinematics file.

	APPLIED_INPUT_KINEMATIC_DOFs = []
	
	for DOF in SPECIFIED_KINEMATIC_DOFs:
		if DOF in INPUT_KINEMATIC_DOFs:
			APPLIED_INPUT_KINEMATIC_DOFs.append( DOF )
		else:
			print
			print 'ERROR: specified kinematic DOFs in configuration file must match input kinematics file!'
			print
			sys.exit(1)
			
	#print APPLIED_INPUT_KINEMATIC_DOFs

	'''
	# FUTURE FEATURE?: change so user can request the opposite DOF in input kinematic file, but use input file DOF.
	# Create DOFs that are actually applied (i.e. for headings found in the input file)
	# if user requests opposite of an input DOF, use input DOF and set inversion factor
	# Set kinematic DOFs that are actually applied (i.e. for headings found in the input
	# file), if the user specifies the opposite DOF as is defined in the input kinematics file.
	for DOF in SPECIFIED_KINEMATIC_DOFs:
		if DOF in INPUT_KINEMATIC_DOFs:
			APPLIED_INPUT_KINEMATIC_DOFs.append( DOF )
		elif OPPOSITE_DOFs[DOF] in INPUT_KINEMATIC_DOFs:
			APPLIED_INPUT_KINEMATIC_DOFs.append( OPPOSITE_DOFs[DOF] )
		else:
			print
			print 'ERROR: invalid desired or inversely defined kinematic DOFs!'
			print
			sys.exit(1)
	'''

	'''
	# INPUT KINEMATIC OFFSETS FOR INPUT KINEMATIC DATA

	# WARNING: conversions/units from user-defined config file should be in model units!
	#          model units: mm, radians, N, N*mm
	# NOTE: may want to make this less ambiguous in the configuration parameter name?
	
	line = ConfigFileLines[ConfigFileLines.index('KINEMATIC OFFSETS DOFs\n')+1].strip('\n')
	line_list = line.split('\'')
	SPECIFIED_KINEMATIC_OFFSETS_DOFs = filter(lambda string: string.strip(), line_list)

	line = ConfigFileLines[ConfigFileLines.index('KINEMATIC OFFSETS (mm,radians)\n')+1].strip('\n')
	line_list = line.split()
	float_line_list = map( float, filter(lambda string: string.strip(), line_list) )

	#print SPECIFIED_KINEMATIC_OFFSETS_DOFs

	# Create DOFs offsets that are actually applied (as above)
	APPLIED_KINEMATIC_OFFSETS_DOFs = []
	APPLIED_KINEMATIC_OFFSETS = {}  # keys: APPLIED_KINEMATIC_OFFSETS_DOFs
	
	for (DOF,offset) in zip(SPECIFIED_KINEMATIC_OFFSETS_DOFs,float_line_list):
		if DOF in INPUT_KINEMATIC_DOFs:
			APPLIED_KINEMATIC_OFFSETS_DOFs.append( DOF )
			APPLIED_KINEMATIC_OFFSETS[ DOF ] = offset
		elif OPPOSITE_DOFs[DOF] in INPUT_KINEMATIC_DOFs:
			APPLIED_KINEMATIC_OFFSETS_DOFs.append( OPPOSITE_DOFs[DOF] )
			APPLIED_KINEMATIC_OFFSETS[ OPPOSITE_DOFs[DOF] ] = -offset
		else:
			print
			print 'ERROR: invalid desired or inversely defined kinematic DOFs!'
			print
			sys.exit(1)

	#print APPLIED_KINEMATIC_OFFSETS

	# TODO: go from applied to flip all if necessary
	
	APPLIED_2_MODEL_KINEMATIC_OFFSETS = {}

	# Invert all applied input DOFs if transfering kinematics/kinetics from input bone to the opposing one
	if INPUT_MOVING_BONE_SIDE!=MODEL_MOVING_BONE_SIDE:
		for DOF in APPLIED_KINEMATIC_OFFSETS_DOFs:
			#KINEMATIC_INV_FACTOR[DOF] = -1.0 * KINEMATIC_INV_FACTOR[DOF]
			INPUT_2_MODEL_KINEMATIC_DOFs.append( OPPOSITE_DOFs[DOF] )
			APPLIED_2_MODEL_KINEMATIC_OFFSETS[ OPPOSITE_DOFs[DOF] ] = -
	else:
		INPUT_2_MODEL_KINEMATIC_DOFs = APPLIED_INPUT_KINEMATIC_DOFs

	for DOF in APPLIED_KINEMATIC_OFFSETS_DOFs:
		if DOF in 

	# TODO: flip from 
	'''


	# FLIP INPUT FILE DOFs TO MODEL COORDINATE SYSTEM (CS) SIDE DOFs
	# if model moving bone is opposite that of input file moving bone.

	INPUT_2_MODEL_KINEMATIC_DOFs = []	

	# Invert all applied input DOFs if transfering kinematics/kinetics from input bone to the opposing one
	if INPUT_MOVING_BONE_SIDE!=MODEL_MOVING_BONE_SIDE:
		for DOF in APPLIED_INPUT_KINEMATIC_DOFs:
			if DOF in ROTATIONAL_DOFs:
				if DOF=='Internal Rotation' or DOF=='External Rotation':
					INPUT_2_MODEL_KINEMATIC_DOFs.append( OPPOSITE_DOFs[DOF] )
				else:
					INPUT_2_MODEL_KINEMATIC_DOFs.append( DOF )
			else:
				INPUT_2_MODEL_KINEMATIC_DOFs.append( OPPOSITE_DOFs[DOF] )
	else:
		INPUT_2_MODEL_KINEMATIC_DOFs = APPLIED_INPUT_KINEMATIC_DOFs

	#print INPUT_2_MODEL_KINEMATIC_DOFs


	# SET INVERSION FACTORS FOR KINEMATIC CURVES

	# - this is accomplished through 2 succesive stages of potential inversion,
	#   which simplifies logic required for all combinations of inversions.
	# - inversions are independent, and can therefore be applied in any order.
	#
	# ORDER OF POTENTIAL DOF INVERSIONS:
	# 1) DOF/'joint bone side' (i.e. proximal/distal) matching:
	#    - invert all applied kinematic/kinetic DOFs if input/model body CSs are opposite
	# 2) DOF/'positive model coordinate system (CS) dimension' matching:
	#    - invert all DOFs if opposite to desired/specified DOF
	#
	# NOTE: don't need to invert DOFs for right/left joint since input DOFs 
	#       are specified in an anatomical reference frame.
	#       => left anatomical DOFs = right anatomical DOFs
	#

	# Initialize kinematic DOF inversion factors
	# relative to the model coordinate systems anatomical DOFs.

	KINEMATIC_INV_FACTOR = {}  # keys: INPUT_2_MODEL_KINEMATIC_DOFs
	for DOF in INPUT_2_MODEL_KINEMATIC_DOFs:
		KINEMATIC_INV_FACTOR[DOF] = 1.0

	# Create list to store kinematic DOFs applied to the model
	APPLIED_MODEL_KINEMATIC_DOFs = []

	for DOF in INPUT_2_MODEL_KINEMATIC_DOFs:
		if DOF in MODEL_DOFs_CS.keys():
			APPLIED_MODEL_KINEMATIC_DOFs.append( DOF )
		else:
			KINEMATIC_INV_FACTOR[DOF] = -1.0 * KINEMATIC_INV_FACTOR[DOF]
			APPLIED_MODEL_KINEMATIC_DOFs.append( OPPOSITE_DOFs[DOF] )

	#print APPLIED_MODEL_KINEMATIC_DOFs
	#print MODEL_DOFs_CS
	#print KINEMATIC_INV_FACTOR

	# Set kinematic DOF scaling factors
	line = ConfigFileLines[ConfigFileLines.index('SPECIFIED KINEMATIC DOFs SCALING FACTORS\n')+1].strip('\n')
	line_list = map(float,line.split())
	KINEMATIC_SCALING_FACTOR = {}  # keys: APPLIED_MODEL_KINEMATIC_DOFs
	for (DOF,scaling_factor) in zip(APPLIED_MODEL_KINEMATIC_DOFs,line_list):
		KINEMATIC_SCALING_FACTOR[DOF] = scaling_factor

	#print KINEMATIC_SCALING_FACTOR

	# Set index order of the kinetic input file headings
	input_2_model_kinematic_DOFs_index = {}  # keys: INPUT_2_MODEL_KINETIC_DOFs
	for INPUT_2_MODEL_DOF in INPUT_2_MODEL_KINEMATIC_DOFs:
		if INPUT_2_MODEL_DOF in input_kinematic_DOFs_index.keys():
			input_2_model_kinematic_DOFs_index[INPUT_2_MODEL_DOF] = input_kinematic_DOFs_index[INPUT_2_MODEL_DOF]
		elif OPPOSITE_DOFs[INPUT_2_MODEL_DOF] in input_kinematic_DOFs_index.keys():
			input_2_model_kinematic_DOFs_index[INPUT_2_MODEL_DOF] = input_kinematic_DOFs_index[ OPPOSITE_DOFs[INPUT_2_MODEL_DOF] ]			
		#else:

	#print input_kinematic_DOFs_index
	#print input_2_model_kinematic_DOFs_index

	# Set input/model unit conversion factors
	# for kinematic DOFs being applied to the model.
	KINEMATIC_UNIT_CONVERSION_FACTOR = {}  # keys: INPUT_2_MODEL_KINEMATIC_DOFs
	
	for DOF in INPUT_2_MODEL_KINEMATIC_DOFs:
		if DOF in SPATIAL_DOFs:
			KINEMATIC_UNIT_CONVERSION_FACTOR[DOF] = EQ_SPATIAL_UNITS[ MODEL_SPATIAL_UNITS ] / EQ_SPATIAL_UNITS[INPUT_SPATIAL_UNITS]
		elif DOF in ROTATIONAL_DOFs:
			KINEMATIC_UNIT_CONVERSION_FACTOR[DOF] = EQ_ROTATIONAL_UNITS[ MODEL_ROTATIONAL_UNITS ] / EQ_ROTATIONAL_UNITS[ INPUT_ROTATIONAL_UNITS ]


	# INPUT/INVERT KINEMATIC CURVES APPROPRIATELY

	# Build datastructure to store model kinematics curves
	kinematic_curves = {}  # keys: APPLIED_MODEL_KINEMATIC_DOFs
	for DOF in APPLIED_MODEL_KINEMATIC_DOFs:
		kinematic_curves[DOF] = []

	# read line by line
	while True:
		line = InputFile.readline()
		if not line:
			break

		# each line stored in a list with 6 floating point values
		# ORIGINAL ASSUMPTION: Medial, Anterior, Superior, Flexion, Valgus, Internal Rotation
		values = map(float,line.strip('\n').split(DELIMITER))

		for (INPUT_2_MODEL_DOF,MODEL_DOF) in zip(INPUT_2_MODEL_KINEMATIC_DOFs,APPLIED_MODEL_KINEMATIC_DOFs):
			# apply specified offsets to the current kinematic DOF
			shifted_kinematic_DOF_value = values[ input_2_model_kinematic_DOFs_index[INPUT_2_MODEL_DOF] ]  # + APPLIED_KINEMATIC_OFFSETS[MODEL_DOF]  # TODO: should this be positive or negative?
			# store the offset kinematic curve for the current DOF
			#kinematic_curves[MODEL_DOF].append( KINEMATIC_INV_FACTOR[INPUT_DOF] * KINEMATIC_SCALING_FACTOR[INPUT_DOF] * KINEMATIC_UNIT_CONVERSION_FACTOR[INPUT_DOF] * shifted_kinematic_DOF_value )
			kinematic_curves[MODEL_DOF].append( KINEMATIC_INV_FACTOR[INPUT_2_MODEL_DOF] * KINEMATIC_UNIT_CONVERSION_FACTOR[INPUT_2_MODEL_DOF] * shifted_kinematic_DOF_value )

	# close kinematics file
	InputFile.close()


	#---------------------#
	# INPUT KINETICS DATA #
	#---------------------#

	# Set the kinetic degrees of freedom (DOF) specified by the user
	line = ConfigFileLines[ConfigFileLines.index('SPECIFIED KINETIC DOFs\n')+1].strip('\r\n')
	line_list = line.split('\'')
	SPECIFIED_KINETIC_DOFs = filter(lambda string: string.strip(), line_list)

	# Open kinematics input file
	InputFile = open(kinetics_filename,'r')

	# Assumption: first line the input kinetic DOFs headings and the order of
	#             their associated data columns.
	line = InputFile.readline().strip('\r\n')
	INPUT_KINETIC_DOFs = line.split(DELIMITER)

	# Set the index order of the kinetic input file headings
	input_kinetic_DOFs_index = {}
	index = 0
	for DOF in INPUT_KINETIC_DOFs:
		input_kinetic_DOFs_index[DOF] = index
		index += 1

	# Create DOFs that are actually applied (i.e. for headings found in the input file).
	# NOTE: user-defined kinetic DOFs in configuration file must match the input kinetics file.
	# FUTURE FEATURE?: change so user can request the opposite DOF in input kinetic file, but use input file DOF.

	APPLIED_INPUT_KINETIC_DOFs = []

	for DOF in SPECIFIED_KINETIC_DOFs:
		if DOF in INPUT_KINETIC_DOFs:
			APPLIED_INPUT_KINETIC_DOFs.append( DOF )
		else:
			print
			print 'ERROR: specified kinetic DOFs in configuration file must match input kinematics file!'
			print
			sys.exit(1)

	'''
	for DOF in SPECIFIED_KINETIC_DOFs:
		if DOF in INPUT_KINETIC_DOFs:
			APPLIED_INPUT_KINETIC_DOFs.append( DOF )
		elif OPPOSITE_DOFs[DOF] in INPUT_KINETIC_DOFs:
			APPLIED_INPUT_KINETIC_DOFs.append( OPPOSITE_DOFs[DOF] )
		else:
			print
			print 'ERROR: invalid desired or inversely defined kinetic DOFs!'
			print
			sys.exit(1)
	'''


	# FLIP INPUT FILE DOFs TO MODEL COORDINATE SYSTEM (CS) SIDE DOFs
	# if model moving bone is opposite that of input file moving bone.

	INPUT_2_MODEL_KINETIC_DOFs = []	

	# Invert all applied input DOFs if transfering kinematics/kinetics from input bone to the opposing one
	if INPUT_MOVING_BONE_SIDE!=MODEL_MOVING_BONE_SIDE:
		for DOF in APPLIED_INPUT_KINETIC_DOFs:
			if DOF in TORQUE_DOFs:
				if DOF=='Internal Rotation Torque' or DOF=='External Rotation Torque':
					INPUT_2_MODEL_KINETIC_DOFs.append( OPPOSITE_DOFs[DOF] )
				else:
					INPUT_2_MODEL_KINETIC_DOFs.append( DOF )
			elif DOF in FORCE_DOFs:
				if DOF=='Compression' or DOF=='Distraction':
					INPUT_2_MODEL_KINETIC_DOFs.append( DOF )
				else:
					INPUT_2_MODEL_KINETIC_DOFs.append( OPPOSITE_DOFs[DOF] )
			#else:

	else:
		INPUT_2_MODEL_KINETIC_DOFs = APPLIED_INPUT_KINETIC_DOFs


	# DETERMINE KINETIC INVERSION FACTORS

	# Initialize kinetic DOF inversion factors
	# relative to the model coordinate systems anatomical DOFs.

	KINETIC_INV_FACTOR = {}  # keys: INPUT_2_MODEL_KINETIC_DOFs
	for DOF in INPUT_2_MODEL_KINETIC_DOFs:
		KINETIC_INV_FACTOR[DOF] = 1.0

	# Create list to store kinetic DOFs applied to the model
	APPLIED_MODEL_KINETIC_DOFs = []

	for DOF in INPUT_2_MODEL_KINETIC_DOFs:
		if DOF in MODEL_DOFs_CS.keys():
			APPLIED_MODEL_KINETIC_DOFs.append( DOF ) 
		else:
			KINETIC_INV_FACTOR[DOF] = -1.0 * KINETIC_INV_FACTOR[DOF]
			APPLIED_MODEL_KINETIC_DOFs.append( OPPOSITE_DOFs[DOF] )


	# Set kinematic DOF scaling factors
	line = ConfigFileLines[ConfigFileLines.index('SPECIFIED KINETIC DOFs SCALING FACTORS\n')+1].strip('\n')
	line_list = map(float,line.split())
	KINETIC_SCALING_FACTOR = {}  # keys: APPLIED_MODEL_KINETIC_DOFs
	for (DOF,scaling_factor) in zip(APPLIED_MODEL_KINETIC_DOFs,line_list):
		KINETIC_SCALING_FACTOR[DOF] = scaling_factor
		

	# Set index order of the kinetic input file headings
	# referenced by transformed 'INPUT_2_MODEL_KINETIC_DOFs' keys.
	input_2_model_kinetic_DOFs_index = {}  # keys: INPUT_2_MODEL_KINETIC_DOFs
	for INPUT_2_MODEL_DOF in INPUT_2_MODEL_KINETIC_DOFs:
		if INPUT_2_MODEL_DOF in input_kinetic_DOFs_index.keys():
			input_2_model_kinetic_DOFs_index[INPUT_2_MODEL_DOF] = input_kinetic_DOFs_index[INPUT_2_MODEL_DOF]
		elif not (INPUT_2_MODEL_DOF in input_kinetic_DOFs_index.keys()):
			input_2_model_kinetic_DOFs_index[INPUT_2_MODEL_DOF] = input_kinetic_DOFs_index[ OPPOSITE_DOFs[INPUT_2_MODEL_DOF] ]			
		#else:


	# set input/model unit conversion factors
	# for DOFs being applied to the model
	KINETIC_UNIT_CONVERSION_FACTOR = {}  # keys: INPUT_2_MODEL_KINETIC_DOFs
	
	for DOF in INPUT_2_MODEL_KINETIC_DOFs:
		if DOF in FORCE_DOFs:
			KINETIC_UNIT_CONVERSION_FACTOR[DOF] = EQ_FORCE_UNITS[ MODEL_FORCE_UNITS ] / EQ_FORCE_UNITS[ INPUT_FORCE_UNITS ]
		elif DOF in TORQUE_DOFs:
			KINETIC_UNIT_CONVERSION_FACTOR[DOF] = EQ_TORQUE_UNITS[ MODEL_TORQUE_UNITS ] / EQ_TORQUE_UNITS[ INPUT_TORQUE_UNITS ]


	# INPUT/INVERT KINETIC CURVES APPROPRIATELY

	# Build kinetics datastructure to store load curves
	kinetic_curves = {}  # keys: APPLIED_MODEL_KINETIC_DOFs
	for DOF in APPLIED_MODEL_KINETIC_DOFs:
		kinetic_curves[DOF] = []

	# read line by line
	while True:
		line = InputFile.readline()
		if not line:
			break

		# each line stored in a list with 6 floating point values
		# ASSUMPTION?: Lateral Drawer, Anterior Drawer, Distraction, Extension Torque, Varus Torque, External Rotation Torque
		values = map(float,line.strip('\n').split(DELIMITER))

		for (INPUT_2_MODEL_DOF,MODEL_DOF) in zip(INPUT_2_MODEL_KINETIC_DOFs,APPLIED_MODEL_KINETIC_DOFs):
			# store the offset kinetic curve for the current DOF
			#kinetic_curves[MODEL_DOF].append( KINETIC_INV_FACTOR[INPUT_2_MODEL_DOF] * KINETIC_SCALING_FACTOR[INPUT_2_MODEL_DOF] * KINETIC_UNIT_CONVERSION_FACTOR[INPUT_DOF] * values[ input_kinetic_DOFs_index[INPUT_DOF] ] )
			kinetic_curves[MODEL_DOF].append( KINETIC_INV_FACTOR[INPUT_2_MODEL_DOF] * KINETIC_UNIT_CONVERSION_FACTOR[INPUT_2_MODEL_DOF] * values[ input_2_model_kinetic_DOFs_index[INPUT_2_MODEL_DOF] ] )


	# close kinetics file
	InputFile.close()


	#===============================#
	# CREATE TIME STEPPER LOADCURVE #
	#===============================#

	# ASSUMPTION: user-specified time is the input kinematic/kinetic DOF curve duration
	
	line = ConfigFileLines[ConfigFileLines.index('KINEMATICS/KINETICS DURATION (s)\n')+1].strip('\n')
	curve_duration = float(line)

	line = ConfigFileLines[ConfigFileLines.index('TIME STEPPER LOADCURVE SUB-STEP DURATION FACTORS\n')+1].strip('\n')
	SUBSTEP_DURATION_FACTORS = map(float,line.split())

	# WARNING: may want to do some other error checking here if the number of
	#          substep may potentially be user defined in the future
	if len(SUBSTEP_DURATION_FACTORS)!=3:
		print
		print 'ERROR: there must be three sub-step timing factors in the input configuration file!'
		print
		sys.exit(1)

	# ASSUMPTION:
	# Full simulation step consists of three successive sub-steps:
	# 1) ramp to initial kinematic/kinetic DOF curve value, over time = T[0]
	# 2) apply kinematic/kinetic DOF curve values, over time = T[1] = curve_duration
	# 3) hold final kinematic/kinetic DOF curve value, for time = T[2]

	T = []
	T.append( SUBSTEP_DURATION_FACTORS[0] * curve_duration )  # time to ramp to initial curve value
	T.append( SUBSTEP_DURATION_FACTORS[1] * curve_duration )  # duration of the curve
	T.append( SUBSTEP_DURATION_FACTORS[2] * curve_duration )  # time to hold final curve value


	# CREATE SOLVER/OUTPUT TIMING LOADCURVE
	# default stepper load curve id number = "1"
	line = ConfigFileLines[ConfigFileLines.index('TIME STEPPER LOADCURVE SUB-STEP NUMBER OF INTERVALS\n')+1].strip('\n')
	#SUBSTEP_OUTPUT_INTERVAL_FACTORS = map(float,line.split())
	SUBSTEP_NUM_INTERVALS = map(int,line.split())

	#LOADCURVE_SUBSTEP_STEP_SIZE = SUBSTEP_OUTPUT_INTERVAL_FACTORS[1]  # assumption 2nd of 3 factors

	# Re-determine an integer number of intervals for each sub-step if the
	# specified substep output interval factors do not evenly divide the substep.
	
	#SUBSTEP_NUM_INTERVALS = []
	#for factor in SUBSTEP_OUTPUT_INTERVAL_FACTORS:
	#	SUBSTEP_NUM_INTERVALS.append( int(1.0/factor) )

	# initalize 
	time_stepper_loadcurve = [0.0]

	# Generate time points for solver/output requests load curve from sub-step definitions
	for substep_index in range(len(SUBSTEP_NUM_INTERVALS)):
		abs_time_last_substep = time_stepper_loadcurve[-1]
		for substep_interval in range(1,SUBSTEP_NUM_INTERVALS[substep_index]+1):
			time_stepper_loadcurve.append( abs_time_last_substep + (float(substep_interval)/SUBSTEP_NUM_INTERVALS[substep_index])*T[substep_index] )

	#print time_stepper_loadcurve


	#==================================#
	# ASSIGN APPLIED DOFs IN FEBio XML #
	#==================================#

	# Merge applied kinematic & kinetic DOF loadcurve datastructures
	# NOTE: This makes logic easier, but may be wasting more space.
	# Merge DOF lists
	APPLIED_INPUT_DOFs = APPLIED_INPUT_KINEMATIC_DOFs + APPLIED_INPUT_KINETIC_DOFs
	APPLIED_MODEL_DOFs = APPLIED_MODEL_KINEMATIC_DOFs + APPLIED_MODEL_KINETIC_DOFs

	#print APPLIED_MODEL_DOFs

	# Duplicate/merge curve dictionaries
	loadcurves = kinematic_curves.copy()
	loadcurves.update( kinetic_curves )

	del kinematic_curves, kinetic_curves


	#---------------------------#
	# CREATE LOAD CURVES IN DOM #
	#---------------------------#

	input_FEBio_filename = ConfigFileLines[ConfigFileLines.index('FEBIO INPUT FILE\n')+1].strip('\n')

	# INPUT FEBio XML DATASTRUCTURE FROM INPUT FILE

	tree = parse(input_FEBio_filename)
	root = tree.getroot()

	#print root.tag
	#print root.attrib

	# obtain reference to LoadData element
	LoadData = root.findall('LoadData')[0]  # WARNING: assumes only one Constraints element in input FEBio file


	# DELETE ALL DEFAULT LOAD CURVES
	# before adding new load curves, delete existing defaults from input FEBio DOM
	loadcurve_list = LoadData.findall('loadcurve')

	for loadcurve_element in loadcurve_list:
		LoadData.remove(loadcurve_element)
	

	# DEFINE SOLVER/OUTPUT TIMING LOADCURVE IN FEBio DOM
	
	# ASSUMPTION: default id number for solver/output timing loadcurve is "1"
	loadcurve = SubElement(LoadData,'loadcurve')
	loadcurve.set('id','1')
	loadcurve.set('type','smooth')

	# define first time point in solver/output timing loadcurve
	# NOTE: this is outside the loop to simplify dt calculation (iterative)

	# set initial time point (0.0?) as the last as initial condition prior to entering the loop.
	t_last = time_stepper_loadcurve[0]

	# create first time point for solver/output timing loadcurve
	# NOTE: this is outside the loop to simplify dt calculation (iterative)
	point = SubElement(loadcurve,'point')
	point.text = str(time_stepper_loadcurve[0]) + ',0.0'


	# Create remaining solver/output timing loadcurve points
	# while finding change in time from prior point (dt)
	for t in time_stepper_loadcurve[1:]:
		
		# determine time duration of the previous time interval
		dt = t-t_last

		# create current solver/output timing loadcurve point in FEBio tree
		point = SubElement(loadcurve,'point')
		point.text = str(t) + ',' + str(dt)

		# advance last time to current
		t_last = t

	#dump(LoadData)

	# increment load curve number for specified DOFs
	# Set initial value
	loadcurve_num = 2  # start number of next DOF loadcurve at "2", following solver/output timing loadcurve number, i.e. 1
	# Create dictionary to store
	loadcurve_DOF_num = {}  # keys: APPLIED_MODEL_DOFs

	# DOF LOADCURVE DEFINING LOOP
	# precondition: DOF loadcurve numbering starts at "2"
	for MODEL_DOF in APPLIED_MODEL_DOFs:

		# initialize time to zero for current DOF curve
		t = 0.0

		# associate the DOF name to the DOF curve number
		# this is done in order to ensure correct DOF curve
		loadcurve_DOF_num[MODEL_DOF] = loadcurve_num

		# create new load curve
		loadcurve = SubElement( LoadData, 'loadcurve' )
		loadcurve.set( 'id', str(loadcurve_DOF_num[MODEL_DOF]) )
		loadcurve.set( 'type', 'smooth' )


		# (1/3) CREATE FIRST SUB-STEP: ramp to initial value of load curve from zero over time t=T1

		point = SubElement(loadcurve,'point')
		point.text = str(t) + ',0.0'


		# (2/3) CREATE SECOND SUB-STEP: input load cuuve

		# LOADCURVE POINT LOOP
		# NOTE: time is calculated absolutely instead of relatively to avoid 
		#       best reduce accumulated errors in numerical precision that
		#       would occur through an iterative scheme.
		for (interval,value) in zip(range(len(loadcurves[MODEL_DOF])),loadcurves[MODEL_DOF]):

			# determine absolute time of current loadcurve value
			t = T[0] + ( float(interval)/(len(loadcurves[MODEL_DOF])-1) ) * T[1]

			# create current point for current loadcurve
			point = SubElement(loadcurve,'point')
			point.text = str(t) + ',' + str(value)

		# END POINT LOOP

		
		# (3/3) CREATE THIRD SUB-STEP: hold final load curve value
		# advance time (t = T1+T2+T3)
		# set absolute time to end of hold sub-step

		t = T[0] + T[1] + T[2]

		
		# Create final point in loadcurve for held final value of specified curve
		point = SubElement(loadcurve,'point')
		point.text = str(t) + ',' + str(loadcurves[MODEL_DOF][-1])

		# increment load curve counter
		loadcurve_num += 1		

	# END LOAD CURVE LOOP


	#========================#
	# SET SOLVER STEP VALUES #
	#========================#

	NUM_SOLVER_TIME_STEPS = int( ConfigFileLines[ConfigFileLines.index('NUMBER OF SOLVER TIME STEPS\n')+1].strip('\n') )

	Step = root.findall('Step')[0]

	Control = Step.findall('Control')[0]

	time_steps = Control.findall('time_steps')[0]
	time_steps.text = str( NUM_SOLVER_TIME_STEPS )

	step_size = Control.findall('step_size')[0]
	step_size.text = str( time_stepper_loadcurve[-1]/NUM_SOLVER_TIME_STEPS )


	#====================#
	# CREATE CONSTRAINTS #
	#====================#

	# Create FEBio DOF-to-Constraint pairings
	# NOTE: this doesn't indicate sign of dimension, just which dimension
	constraint_4_DOF = {}  # NOTE: this should be the same as MODEL_DOFs_CS dictionary

	# kinematic DOFs
	constraint_4_DOF['Medial'] = 'x'
	constraint_4_DOF['Anterior'] = 'y'
	constraint_4_DOF['Superior'] = 'z'
	constraint_4_DOF['Flexion'] = 'Rx'
	constraint_4_DOF['Valgus'] = 'Ry'
	constraint_4_DOF['Internal Rotation'] = 'Rz'
	# ... and their opposites
	constraint_4_DOF['Lateral'] = 'x'
	constraint_4_DOF['Posterior'] = 'y'
	constraint_4_DOF['Inferior'] = 'z'
	constraint_4_DOF['Extension'] = 'Rx'
	constraint_4_DOF['Varus'] = 'Ry'
	constraint_4_DOF['External Rotation'] = 'Rz'

	# kinetic DOFs	
	constraint_4_DOF['Lateral Drawer'] = 'x'
	constraint_4_DOF['Anterior Drawer'] = 'y'
	constraint_4_DOF['Distraction'] = 'z'
	constraint_4_DOF['Extension Torque'] = 'Rx'
	constraint_4_DOF['Varus Torque'] = 'Ry'
	constraint_4_DOF['External Rotation Torque'] = 'Rz'
	# ... and their opposites
	constraint_4_DOF['Medial Drawer'] = 'x'
	constraint_4_DOF['Posterior Drawer'] = 'y'
	constraint_4_DOF['Compression'] = 'z'
	constraint_4_DOF['Flexion Torque'] = 'Rx'
	constraint_4_DOF['Valgus Torque'] = 'Ry'
	constraint_4_DOF['Internal Rotation Torque'] = 'Rz'


	#------------------------------#
	# CREATE KINEMATIC CONSTRAINTS #
	#------------------------------#

	MODEL_MOVING_BONE_NUM = ConfigFileLines[ConfigFileLines.index('MODEL KINEMATICS/KINETICS BONE MAT NUMBER\n')+1].strip('\n')

	# obtain reference to Constraints element
	#Constraints = root.findall('.//Constraints')[0]  # WARNING: assumes only one Constraints element in input FEBio file
	Constraints = Step.findall('Constraints')[0]  # WARNING: assumes only one Constraints element in input FEBio file

	# DELETE ALL DEFAULT LOAD CURVES
	# before adding new load curves, delete existing defaults from input FEBio DOM
	rigid_body_list = Constraints.findall('rigid_body')

	for rigid_body_element in rigid_body_list:
		if rigid_body_element.attrib['mat']==MODEL_MOVING_BONE_NUM:
			Constraints.remove(rigid_body_element)
	

	# CREATE SPECIFIED KINEMATIC CONSTRAINTS (for moving bone)

	# create new load curve
	rigid_body = SubElement( Constraints, 'rigid_body' )
	rigid_body.set( 'mat', MODEL_MOVING_BONE_NUM )


	# create prescribed (kinematic) constraints in rigid body definition for moving bone in model
	for MODEL_DOF in APPLIED_MODEL_KINEMATIC_DOFs:

		prescribed_constraint = SubElement( rigid_body, 'prescribed' )
		prescribed_constraint.set( 'bc', constraint_4_DOF[MODEL_DOF] )
		prescribed_constraint.set( 'lc', str(loadcurve_DOF_num[MODEL_DOF]) )
		prescribed_constraint.text = str( KINEMATIC_SCALING_FACTOR[MODEL_DOF] )


	#----------------------------#
	# CREATE KINETIC CONSTRAINTS #
	#----------------------------#

	# add kinetic constraint definitions into same rigid body constraint element as kinematic constraints

	# create force (kinetic) constraints in rigid body definition for moving bone in model
	for MODEL_DOF in APPLIED_MODEL_KINETIC_DOFs:

		force_constraint = SubElement(rigid_body,'force')
		force_constraint.set( 'bc', constraint_4_DOF[MODEL_DOF] )
		force_constraint.set( 'lc', str(loadcurve_DOF_num[MODEL_DOF]) )
		force_constraint.text = str( KINETIC_SCALING_FACTOR[MODEL_DOF] )


	#===============================#
	# OUTPUT UPDATED FEBio XML FILE #
	#===============================#

	# OUTPUT FEBio FILE WITH PRESCRIBED EXPERIMENTAL MOTION
	output_FEBio_filename = input_FEBio_filename[:-4] + '_MOTION.feb'

	OUTPUT_XML_INDENT_STRING = "\t"
	
	#root_element = Element(root)

	indent_elements(root,OUTPUT_XML_INDENT_STRING)
	
	# TODO: get input XML file encoding to properly rewrite output file

	tree.write(output_FEBio_filename, xml_declaration=True, encoding='ISO-8859-1', method="xml")

	#xmlstr = ElementTree.tostring(tree, encoding='utf8', method='xml')


# END OF prescribe_motion_4_FEBio()


# Default calls image_2_mesh() if this script is not called as a function
if __name__ == "__main__":
	prescribe_motion_4_FEBio(sys.argv)