""" Post4MOR in Python.

Evgenii Rudnyi (c) 2007, http://Evgenii.Rudnyi.Ru/

http://ModelReduction.com/soft/Post4MOR/ 

This software is a copyrighted work licensed under the terms, described
in the file "FREE_LICENSE". 

This is an implementation in Python functions from Post4MOR.m. The
design is the same, however as Python and Mathematica are quite
different from each other, there are some differences.

A good starting point is the documentation for Post4MOR. See
Post4MORmanual.nb.pdf in the parent folder.

One difference is that Python does not support function
overloading. The solution was to introduce new function names.

A related problem is that there is almost no checks for correctness of
parameters. If you see that Python complains about these functions,
first thing is to check that parameters are correct.

The function that reads Matrix Market format in SciPy assumes very
strict rules. It does not like empty lines that MOR for ANSYS writes.
The solution was to introduce a second parameter in ReadSystem to
remove empty lines. If

  ReadSystem('mor')

fails, please try

  ReadSystem('mor', True)

This way empty lines will be removed from files mor.*.

Other limits:
	Difference requires that XSeries are the same.
	TrasformResult does not take functions.
	TransientSimulation is not implemented. Use AnsysTransientSimulation.
	Input function can depend only on time/frequency.

Below variable __all__ describes publicly available function names. In
each public function there is short help. Use help(function_name).

Requirements:

	Python 2.5 - http://www.python.org
	SciPy and NumPy for linear algebra - http://www.scipy.org
	matplotlib for plotting - http://matplotlib.sourceforge.net/

Changes:
	13 September 2007
		ReadSystemOld, a new parameter to deal with phase angle
	20 September 2007
		Fix in StationarySolution
		ReadSystem(sys, , Format='MatrixMarket')
		ReadResultTable
	25 September 2007
		Fix in ReadSystem for Caspoc
	11 November 2007 (thanks Mr Huebner)
		Fix in AddDampingRatio 
		WriteResultAsDirichlet - comment out printing resume, removing debug printing
		PlotResult - more options
		PlotConvergence - more options
		WriteResult and WriteResultTable - fixed
		New functions ReadResultGrace, WriteResultGrace
		MakeSimulationResult - inserted transformation to array
	11 December 2007
		Change in ToFirstOrderSystem and AddOutputs - fixing a bug (r_ not defined)
	23 February 2008
	  Change in TakeSystem - optional parameter start
	24 April 2008
	  TransformResult - the name was written wrong
	19 May 2008
		Check for the sparse matrix in readMatrix (this is a new function)
	12 June 2008
	  Bug fix in checkConsistence (for MatrixD)
		Now MatrixD is taken into account in HarmonicSolution
	18 June 2008
		HarmonicSimulation without MatrixD - the implementation was not general
	30 July 2008
	  DeleteOutputs - bug fix
	6 Mai 2009	
		readMatrix - bug fix for the case Format = Table
		DeleteOutputs - bug fix
	6 August 2009
	  AnsysTransientSolution support option Velocity to return velocities
		WriteSystem support now Format=Table
	12 August 2009
	  WriteResult is corrected to work with the previous change
	16 August 2009
	  ToDiagonalSystem
		WriteSystem support now Format=Spice
	12 September 2009
	  ToPort
	20 September 2009
	  StationarySolution has now an option ReturnState
	26 September 2009
		PlotResultTogether
	3 October 2009
	  Now HarmonicSolution has SaveState and Title
	31 January 2010
	  PlotPoles
		ComputePoles
	26 September 2010
		Format = 'Simplorer' in WriteSystem (works only if ni == no)	  
	14 November 2010
		Format = 'VHDL-AMS' in WriteSystem (works only if ni == no)
	26 March 2011
		HarmonicSolutionTest E and F scales by w^2 for Klaus
	31 July 2011
	  SaveState in AnsysTransientSolution (only for the first order)
		ExpandMOR
	15 August 2011
	  Format = 'VHDL-AMS' for piezo
	30 November 2011
	  WriteCSV

"""

import numpy
import scipy.linalg
import scipy.io
import os
import operator
import pylab
import math

degree = math.pi/180.

__all__ = ['MakeDynamicSystem', 
		'MatrixMQ', 
		'MatrixEQ', 
		'MatrixDQ', 
		'FirstOrderSystemQ', 
		'SecondOrderSystemQ', 
		'ExplicitSystemQ', 
		'MatrixM', 
		'MatrixE', 
		'MatrixK', 
		'MatrixA', 
		'MatrixB', 
		'MatrixC', 
		'OutputNames', 
		'MatrixD', 
		'MakeFirstOrderSystem', 
		'MakeSecondOrderSystem', 
		'DeleteOutputs', 
		'AddOutputs', 
		'DeleteDamping', 
		'AddDampingMatrix', 
		'AddDamping', 
		'AddDampingRatio', 
		'AddDampingRatioMatrix',
		'ToSingleInputSystem',
		'ToFirstOrderSystem',
		'ToExplicitSystem',
		'ToDiagonalSystem',
		'ToPort',
		'TakeSystem',
		'ReadSystem',
		'WriteSystem',
		'MakeSimulationResult',
		'XSeries',
		'XName',
		'YSeries',
		'YNames',
		'TransformResult',
		'CompatibleResultQ',
		'Difference',
		'TakeResult',
		'ReadResult',
		'WriteResult',
		'ReadResultTable',
		'WriteResultTable',
		'WriteResultCSV',
		'ReadResultOld',
		'ReadResultGrace',
		'WriteResultGrace',
		'WriteResultAsDirichlet',
		'PlotResult',
		'PlotResultTogether',
		'ComputePoles',
		'PlotPoles',
		'StationarySolution',
		'HarmonicSolution',
		'HarmonicSolutionTest',
		'FrequencyConvergence',
		'RelativeError',
		'LocalError',
		'LocalErrorIndicator',
		'PlotConvergence',
		'AnsysTransientSolution',
		'ExpandMOR']

def checkConsistence(sys):
	dim = MatrixK(sys).shape
	if len(sys) != 7:
		raise BaseException('checkConsistence: The length of the list is not 7')
	elif MatrixMQ(sys) and originalMatrixM(sys).shape != dim:
		raise BaseException('checkConsistence: MatrixM and MatrixK are not compatible')
	elif MatrixEQ(sys) and originalMatrixE(sys).shape != dim:
		raise BaseException('checkConsistence: MatrixE and MatrixK are not compatible')
	elif dim[0] != MatrixB(sys).shape[0]:
		raise BaseException('checkConsistence: MatrixB and MatrixK are not compatible')
	elif dim[0] != MatrixC(sys).shape[1]:
		raise BaseException('checkConsistence: MatrixC and MatrixK are not compatible')
	elif len(MatrixC(sys)) != len(OutputNames(sys)):
		raise BaseException('checkConsistence: MatrixC and OutputNames are not compatible')
	elif MatrixDQ(sys) and originalMatrixD(sys).shape[0] !=  MatrixC(sys).shape[0]:
		raise BaseException('checkConsistence: MatrixD and MatrixC are not compatible')
	elif MatrixDQ(sys) and originalMatrixD(sys).shape[1] !=  MatrixB(sys).shape[1]:
		raise BaseException('checkConsistence: MatrixD and MatrixB are not compatible')
	else:
		return sys

def MakeDynamicSystem(matM, matE, matK, matB, matC, namesC, matD = None):
	"""gives a DynamicSystem object. 

	matM, matE, and matD should be matrices or None. matK, matB, and
	matC must be matrices. namesC is a list of strings. matD is optional
	and by default is None. """
	return checkConsistence([matM, matE, matK, matB, matC, namesC, matD])

def originalMatrixM(sys):
	return sys[0]

def originalMatrixE(sys):
	return sys[1]

def originalMatrixD(sys):
	return sys[6]

def MatrixMQ(sys):
	"""gives True if system matrix M is present in sys and False otherwise."""
	return sys[0] != None

def MatrixEQ(sys):
	"""gives True if system matrix E is present in sys and False otherwise."""
	return sys[1] != None

def MatrixDQ(sys):
	"""gives True if system matrix D is present in sys and False otherwise."""
	return sys[6] != None

def FirstOrderSystemQ(sys):
	"""gives True if sys represents a first order system and False otherwise."""
	return sys[0] == None

def SecondOrderSystemQ(sys):
	"""gives True if sys represents a second order system and False otherwise."""
	return sys[0] != None

def ExplicitSystemQ(sys):
	"""gives True if sys represents an explicit first order system and False otherwise."""
	return sys[0] == None and sys[1] == None

def MatrixM(sys):
	"""gives system matrix M of sys."""
	if MatrixMQ(sys):
		return sys[0]
	else:
		return numpy.zeros_like(MatrixK(sys))

def MatrixE(sys):
	"""gives system matrix E of sys."""
	if MatrixEQ(sys):
		return sys[1]
	else:
		if FirstOrderSystemQ(sys):
			return numpy.identity(len(MatrixK(sys)), MatrixK(sys).dtype)
		else:
			return numpy.zeros_like(MatrixK(sys))

def MatrixK(sys):
	"""gives system matrix K of sys."""
	return sys[2]

def MatrixA(sys):
	"""gives system matrix A of sys."""
	return -sys[2]

def MatrixB(sys):
	"""gives system matrix B of sys."""
	return sys[3]

def MatrixC(sys):
	"""gives system matrix C of sys."""
	return sys[4]

def OutputNames(sys):
	"""gives a list of output names of sys."""
	return sys[5]

def MatrixD(sys):
	"""gives system matrix D of sys."""
	if MatrixDQ(sys):
		return sys[6]
	else:
		return numpy.zeros((MatrixC(sys).shape[0], MatrixB(sys).shape[1]), MatrixK(sys).dtype)

def MakeFirstOrderSystem(matE, matA, matB, matC, namesC, matD = None):
	"""give a DynamicSystem object for a first order system. 
	
	matE, and matD should be matrices or None. matA, matB, and matC must
	be matrices. namesC is a list of strings. matD is optional and by
	default is None."""
	return MakeDynamicSystem(None, matE, -matA, matB, matC, namesC, matD)

def MakeSecondOrderSystem(matM, matE, matK, matB, matC, namesC, matD = None):
	"""give a DynamicSystem object for a second order system. 
	
	matM, matE, and matD should be matrices or None. matK, matB, and matC must
	be matrices. namesC is a list of strings. matD is optional and by
	default is None."""
	return MakeDynamicSystem(matM, matE, matK, matB, matC, namesC, matD)

def DeleteOutputs(sys, names):
	"""gives a new system.
	
	Outputs defined in a list of names are deleted."""
	ind = []
	out = OutputNames(sys)[:]
	for n in names:
		i = out.index(n)
		ind.append(i)
	for n in names:
		i = out.index(n)
		del out[i]
	return MakeDynamicSystem(   
		originalMatrixM(sys),
		originalMatrixE(sys),
		MatrixK(sys), 
		MatrixB(sys), 
		numpy.delete(MatrixC(sys), ind, 0), 
	  out, 
		originalMatrixD(sys)
	)

def AddOutputs(sys, matCadd, names):
	"""gives a new system.
	
	Additional outputs defined by rows of matCadd and a list of names
	will be inserted."""
	out = OutputNames(sys)[:]
	out.extend(names)
	return MakeDynamicSystem(
		originalMatrixM(sys),
		originalMatrixE(sys),
		MatrixK(sys), 
		MatrixB(sys), 
		numpy.r_[MatrixC(sys), matCadd], 
		out, 
		originalMatrixD(sys)
	) 

def DeleteDamping(sys):
	"""gives a new system without system matrix E. 
	
	sys must be a second order system."""
	if FirstOrderSystemQ(sys):
		return sys
	else:
		return MakeDynamicSystem(
			MatrixM(sys),
			None,
			MatrixK(sys),
			MatrixB(sys),
			MatrixC(sys),
			OutputNames(sys),
			originalMatrixD(sys)
		)

def AddDampingMatrix(sys, matE):
	"""gives a new system with matrix E inreased by matE.
	
	sys must be a second order system."""
	if FirstOrderSystemQ(sys):
		return sys
	else:
		return MakeDynamicSystem(
			MatrixM(sys),
			MatrixE(sys) + matE,
			MatrixK(sys),
			MatrixB(sys),
			MatrixC(sys),
			OutputNames(sys),
			originalMatrixD(sys)
		)

def AddDamping(sys, alpha, beta):
	"""gives a new system.
	
	Matrix E is increased by Rayleigh damping alpha*matM+beta*matK.
	sys must be a second order system."""
	if FirstOrderSystemQ(sys):
		return sys
	else:
		return MakeDynamicSystem(
			MatrixM(sys),
			MatrixE(sys) + alpha*MatrixM(sys) + beta*MatrixK(sys),
			MatrixK(sys),
			MatrixB(sys),
			MatrixC(sys),
			OutputNames(sys),
			originalMatrixD(sys)
		)

def AddDampingRatio(sys, dmprat):
	"""gives a new system.
	
	Matrix K is increased by 2j*dmprat*MatrixK. sys must be a second order
	system."""
	if FirstOrderSystemQ(sys):
		return sys
	else:
		return MakeDynamicSystem(
			MatrixM(sys),
			originalMatrixE(sys),
			MatrixK(sys) + 2j*dmprat*MatrixK(sys),
			MatrixB(sys),
			MatrixC(sys),
			OutputNames(sys),
			originalMatrixD(sys)
		)

def AddDampingRatioMatrix(sys, imatK):
	"""gives a new system with matrix K inreased by 1j*imatK.

	sys must be a second order system."""
	if FirstOrderSystemQ(sys):
		return sys
	else:
		return MakeDynamicSystem(
			MatrixM(sys),
			originalMatrixE(sys),
			MatrixK(sys) + j*imatK,
			MatrixB(sys),
			MatrixC(sys),
			OutputNames(sys),
			originalMatrixD(sys)
		)

def ToSingleInputSystem(sys):
	"""gives a new system in which there is only first input left."""
	return MakeDynamicSystem(
		originalMatrixM(sys),
		originalMatrixE(sys),
		MatrixK(sys), 
		MatrixB(sys)[:,:1], 
		MatrixC(sys), 
		OutputNames(sys), 
		originalMatrixD(sys)
	) 

def ToFirstOrderSystem(sys):
	"""gives a new system in the form of a first order system."""
	if FirstOrderSystemQ(sys):
		return sys
	else:
		zero = numpy.zeros_like(MatrixK(sys))
		one = numpy.identity(len(MatrixK(sys)), MatrixK(sys).dtype)
		r_ = numpy.r_
		return MakeFirstOrderSystem(
			r_[r_['1', MatrixE(sys), MatrixM(sys)], r_['1', -one, zero]],
			-r_[r_['1', MatrixK(sys), zero], r_['1', zero, one]],
			r_[MatrixB(sys), numpy.zeros_like(MatrixB(sys))],
			r_['1', MatrixC(sys), numpy.zeros_like(MatrixC(sys))],
			OutputNames(sys),
			originalMatrixD(sys)
		)

def ToFirstOrderSystemNew(sys):
	"""gives a new system in the form of a first order system."""
	if FirstOrderSystemQ(sys):
		return sys
	else:
		zero = numpy.zeros_like(MatrixK(sys))
		one = numpy.identity(len(MatrixK(sys)), MatrixK(sys).dtype)
		r_ = numpy.r_
		return MakeFirstOrderSystem(
			r_[r_['1', MatrixM(sys), zero], r_['1', zero, one]],
			-r_[r_['1', MatrixE(sys), MatrixK(sys)], r_['1', -one, zero]],
			r_[MatrixB(sys), numpy.zeros_like(MatrixB(sys))],
			r_['1', numpy.zeros_like(MatrixC(sys)), MatrixC(sys)],
			OutputNames(sys),
			originalMatrixD(sys)
		)

def ToExplicitSystem(sys):
	"""gives a new system in the form of an explicit first order system."""
	if FirstOrderSystemQ(sys):
		if MatrixEQ(sys):
			invE = scipy.linalg.inv(MatrixE(sys))
			return MakeFirstOrderSystem(
				None,
				numpy.dot(invE, MatrixA(sys)),
				numpy.dot(invE, MatrixB(sys)),
				MatrixC(sys),
				OutputNames(sys),
				originalMatrixD(sys)
			)
		else:
			return sys
	else:
		return ToExplicitSystem(ToFirstOrderSystem(sys))

def ToDiagonalSystem(sys):
	"""gives a new system: explicit and A is diagonal."""
	if FirstOrderSystemQ(sys):
		D, V = scipy.linalg.eig(MatrixK(sys), MatrixE(sys))
		E = numpy.dot(V.transpose(), numpy.dot(MatrixE(sys), V))
		scale = numpy.sqrt(1/E.diagonal())
		V = numpy.dot(V, numpy.diag(scale))
		K = numpy.dot(V.transpose(), numpy.dot(MatrixK(sys), V))
		K = numpy.diag(K.diagonal())
		B = numpy.dot(V.transpose(), MatrixB(sys))
		C = numpy.dot(MatrixC(sys), V)
		return MakeFirstOrderSystem(
			None,
			-K,
			B,
			C,
			OutputNames(sys),
			originalMatrixD(sys)
		)
	else:
		return ToDiagonalSystem(ToFirstOrderSystem(sys))

def ToPort(sys, names = []):
	"""gives a new system: number of inputs are equal to number of outputs."""
	nin = MatrixB(sys).shape[1]
	nout = MatrixC(sys).shape[0]
	dim = len(MatrixK(sys))
	matB = MatrixB(sys)
	matC = MatrixC(sys)
	namesC = OutputNames(sys)[:]
	namesC.extend(names)
	if nout > nin:
		matB = 	numpy.r_['1', MatrixB(sys), numpy.zeros((dim, nout - nin))]
	elif nin > nout:
		matC = numpy.r_[MatrixC(sys), numpy.zeros((nin - nout, dim))]
	return MakeDynamicSystem(originalMatrixM(sys), 
			originalMatrixE(sys), 
			MatrixK(sys), 
			matB, 
			matC, 
			namesC, 
			originalMatrixD(sys))

def TakeSystem(sys, dim, start = 0):
	"""gives a new system of dimension dim. 
	
	This operation makes sense only when an original system has been
	obtained by iterative model reduction method."""
	if MatrixMQ(sys):
		matM = MatrixM(sys)[start:dim, start:dim]
	else:
		matM = None
	if MatrixEQ(sys):
		matE = MatrixE(sys)[start:dim, start:dim]
	else:
		matE = None
	return MakeDynamicSystem(
		matM,
		matE,
		MatrixK(sys)[start:dim, start:dim], 
		MatrixB(sys)[start:dim, :], 
		MatrixC(sys)[:, start:dim], 
		OutputNames(sys), 
		originalMatrixD(sys)
	) 

def fileQ(name):
	return os.path.isfile(name)

def removeEmptyLines(name):
	f = open(name)
	s = []
	for line in f:
		if len(line) and (not line.isspace()):
			s.append(line)
	f.close()
	f = open(name, 'w')
	f.writelines(s)
	f.close()

def readMatrix(name, remove = False, Format = 'MatrixMarket'):
	if remove:
		removeEmptyLines(name)
	if Format == 'MatrixMarket':
		mat = scipy.io.mmread(name)
		if not isinstance(mat, type(numpy.array(1))):
			mat = mat.toarray()
	else:
		mat = scipy.io.read_array(name)
		if len(mat.shape) == 1:
			mat = numpy.array([mat])
			if name[-1] == 'B':
				mat = mat.transpose()
	return mat

def readMatrix1(name, remove = False, Format = 'MatrixMarket'):
	if fileQ(name):
		return readMatrix(name, remove, Format)
	else:
		raise BaseException('ReadSystem: no file ' + name)

def readMatrix2(name, remove = False, Format = 'MatrixMarket'):
	if fileQ(name):
		return readMatrix(name, remove, Format)
	else:
		return None

def readNames(name):
	if not fileQ(name):
		raise BaseException('ReadSystem: no file ' + name)
	f = open(name)
	n = f.read().split()
	f.close()
	return n

def readFirstOrder(baseName, remove, Format):
	return MakeFirstOrderSystem(
		readMatrix2(baseName + ".E", remove, Format),
		readMatrix1(baseName + ".A", remove, Format),
		readMatrix1(baseName + ".B", remove, Format),
		readMatrix1(baseName + ".C", remove, Format),
		readNames(baseName + ".C.names"),
		readMatrix2(baseName + ".D", remove, Format)
	)

def readSecondOrder(baseName, remove, Format):
	return MakeSecondOrderSystem(
		readMatrix1(baseName + ".M", remove, Format),
		readMatrix2(baseName + ".E", remove, Format),
		readMatrix1(baseName + ".K", remove, Format),
		readMatrix1(baseName + ".B", remove, Format),
		readMatrix1(baseName + ".C", remove, Format),
		readNames(baseName + ".C.names"),
		readMatrix2(baseName + ".D", remove, Format)
	)

def ReadSystem(baseName, remove = False, Format = 'MatrixMarket'):
	"""reads a system from files with the base name baseName.

	When remove is True, then empty lines are removed from files. This
	is necessary, as the function that reads Matrix Market format in
	SciPy assumes very strict rules. It does not like empty lines that
	MOR for ANSYS writes."""
	if fileQ(baseName + ".A"):
		return readFirstOrder(baseName, remove, Format)
	elif fileQ(baseName + ".K"):
		return readSecondOrder(baseName, remove, Format)
	else:
		raise BaseException("ReadSystem: no files with A or K matrix")

def writeNames(name, list):
	f = open(name, 'w')
	for si in list:
		f.write(si + '\n')
	f.close()

def writeMatrix(name, mat, Format):
	f = open(name, 'w')
	if Format == 'MatrixMarket':
		scipy.io.mmwrite(f, mat, precision=17)
	else:
		scipy.io.write_array(f, mat, precision=17)
	f.close()

def writeSpice(sys, basename, name):
	if not (ExplicitSystemQ(sys) and (numpy.diag(MatrixK(sys).diagonal()) == MatrixK(sys)).all()):
		writeSpice(ToDiagonalSystem(sys), basename, name)
		return
	f = open(basename, 'w')
# header
	f.write('.SUBCKT ' + name + ' ')
	for i in OutputNames(sys):
		f.write(i + ' ')
	f.write('pad\n\n')
# preparing information
	nz = len(MatrixA(sys))
	nu = MatrixB(sys).shape[1]
	ny = len(MatrixC(sys))
	if nu != ny:
		raise BaseException('Number of inputs is not equal to the number of outputs')
	R = 1./MatrixK(sys).diagonal()
# generate netlist
	iu = [['pad' for i in range(nu)]]
	for i in range(1, 2*nz):
		t = []
		for j in range(1, nu + 1):
			t.append(str((i-1)*nu+j))
		iu.append(t)
	t = [i for i in OutputNames(sys)]
	iu.append(t)

	# nodes per state variable 
	iz = [2*nz*nu + i for i in range(1, nz + 1)]
	f.write('* thermal channels *\n')
	for j in range(1, nu + 1):
		for i in range(1, nz + 1):
			f.write('\nF_' + str(i) + '_' + str(j) + ' ' + str(iz[i - 1]) + ' 0 VF_' + str(i) + '_' + str(j) + ' ' + repr(MatrixB(sys)[i - 1, j - 1]) + '\n')
			f.write('VF_' + str(i) + '_' + str(j) + ' ' + iu[i*2 - 2][j - 1] + ' ' + iu[i*2 - 1][j - 1] + ' 0V\n');
			f.write('E_' + str(j) + '_' + str(i) + ' ' + iu[i*2][j - 1] + ' ' + iu[i*2 - 1][j - 1] + ' ' + str(iz[i - 1]) + ' 0 ' + repr(MatrixC(sys)[j - 1, i -1]) + '\n')
	f.write('\n* RC cells *\n')
	for i in range(1, nz + 1):
		f.write('\nR_' + str(i) + ' ' + str(iz[i - 1]) + ' 0 ' + repr(R[i - 1]) + '\n')
		f.write('C_' + str(i) + ' ' + str(iz[i - 1]) + ' 0 1 IC=0\n')

	f.write('\n.ENDS\n')
	f.close()

def writeSimplorer(sys, basename):
	if not ExplicitSystemQ(sys):
		writeSimplorer(ToExplicitSystem(sys), basename)
		return
	f = open(basename + '.sml', 'w')
	nz = len(MatrixA(sys))
	nu = MatrixB(sys).shape[1]
	ny = len(MatrixC(sys))
	if nu != ny:
		raise BaseException('Number of inputs is not equal to the number of outputs')
	f.write('SMLDEF SimplorerThermalPort\n')
	f.write('{\n')
	for i in OutputNames(sys):
		f.write('PORT Thermal : ' + i + ';\n')
	f.write('PORT Thermal : T_REF;\n')
	f.write('PORT STRING in: InitVals =')
	for i in range(nz):
		f.write(' 0')
	f.write(';\n\n')
	f.write('INTERN ZStateSpaceModel sssm')
	for i in range(ny):
		f.write(' N' + str(i) + ':=' + OutputNames(sys)[i] + ',')
	f.write(' N' + str(ny) + ':=T_REF\n')
	f.write(' (RefType:= 1\n')
	f.write(', NumStates:= ' + str(nz) + '\n')
	f.write(', InitialValues:=InitVals\n')
	f.write(', MatrixA:=')
	for i in range(nz):
		for j in range(nz):
			f.write(' ' + repr(MatrixA(sys)[i, j]))
	f.write('\n')
	f.write(', MatrixB:=')
	for i in range(nz):
		for j in range(nu):
			f.write(' ' + repr(MatrixB(sys)[i, j]))
	f.write('\n')
	f.write(', MatrixC:=')
	for i in range(ny):
		for j in range(nz):
			f.write(' ' + repr(MatrixC(sys)[i, j]))
	f.write('\n')
	f.write(', MatrixD:=')
	for i in range(ny):
		for j in range(ny):
			f.write(' 0')
	f.write('\n')
	f.write(');\n\n')
	f.write('}\n')
	f.close()

def writeVHDLAMS(sys, basename):
	if FirstOrderSystemQ(sys):
		writeVHDLAMSFO(sys, basename)
	else:
		writeVHDLAMSSO(sys, basename)

def writeVHDLAMSFO(sys, basename):
	nz = len(MatrixA(sys))
	nu = MatrixB(sys).shape[1]
	ny = len(MatrixC(sys))
	if nu != ny:
		raise BaseException('Number of inputs is not equal to the number of outputs')
	f = open(basename + '.vhd', 'w')
	f.write('LIBRARY IEEE;\n')
	f.write('USE IEEE.THERMAL_SYSTEMS.ALL;\n\n')
	f.write('ENTITY thermal_port IS\n')
	f.write('  PORT(\n')
	f.write('    QUANTITY t0 : IN TEMPERATURE := 273.15;\n')
	f.write('    TERMINAL ')
	for i in range(nu):
		f.write('th' + str(i + 1))
		if i != nu - 1:
			f.write(',')
		f.write(' ')
	f.write(': THERMAL);\n')
	f.write('END ENTITY thermal_port;\n\n')
	f.write('ARCHITECTURE behav OF thermal_port IS\n')
	for i in range(1, nu + 1):
		f.write('  QUANTITY t' + str(i) + ' ACROSS h' + str(i) + ' THROUGH th' + str(i) + ';\n')
	for i in range(1, nz + 1): 
		f.write('  QUANTITY z' + str(i) + ': REAL;\n')
	f.write('BEGIN\n')
	f.write('	IF (domain = quiescent_domain) USE\n')
	for i in range(1, nz + 1): 
		f.write('    z' + str(i) + ' == 0.0;\n')
	f.write('	ELSE\n')
	for i in range(nz):
		f.write('    ')
		for j in range(nz):
			f.write(repr(MatrixE(sys)[i, j]) + '*z' + str(j + 1) + '\'dot')
			if (j != nz - 1):
				f.write(' + ')
			else:
				f.write(' == ')
		for j in range(nz):
			f.write(repr(MatrixA(sys)[i, j]) + '*z' + str(j + 1) + ' + ')
		for j in range(nu):
			f.write(repr(MatrixB(sys)[i, j]) + '*h' + str(j + 1))
			if (j != nu - 1):
				f.write(' + ')
		f.write(';\n')
	f.write('	END USE;\n')
	for i in range(nu):
		f.write('  t' + str(i + 1) + ' == t0 + ')
		for j in range(nz):
			f.write(repr(MatrixC(sys)[i, j]) + '*z' + str(j + 1))
			if (j != nz - 1):
				f.write(' + ')
		f.write(';\n')
	f.write('END ARCHITECTURE behav;\n')
	f.close()

def writeVHDLAMSSO(sys, basename):
	nz = len(MatrixA(sys))
	nu = MatrixB(sys).shape[1]
	ny = len(MatrixC(sys))
	if nu != ny:
		raise BaseException('Number of inputs is not equal to the number of outputs')
	if nu != 2:
		raise BaseException('Number of inputs is not equal to two')
	f = open(basename + '.vhd', 'w')
	f.write('LIBRARY IEEE;\n')
	f.write('USE IEEE.ELECTRICAL_SYSTEMS.ALL;\n')
	f.write('USE IEEE.MECHANICAL_SYSTEMS.ALL;\n\n')
	f.write('ENTITY piezo_port IS\n')
	f.write('  GENERIC(\n')
	f.write('    Vcon : REAL := 1.0;\n')
	f.write('    Scon : REAL := 1.0;\n')
	f.write('    Ccon : REAL := 1.0);\n')
	f.write('  PORT(\n')
	f.write('    TERMINAL pos1, pos2 : TRANSLATIONAL;\n')
	f.write('    TERMINAL p, m : ELECTRICAL);\n')
	f.write('END ENTITY piezo_port;\n\n')
	f.write('ARCHITECTURE behav OF piezo_port IS\n')
	f.write('  QUANTITY pos ACROSS f THROUGH pos1 TO pos2;\n')
	f.write('  QUANTITY v ACROSS i THROUGH p TO m;\n')
	f.write('  QUANTITY q: REAL;\n')
	for i in range(1, nz + 1): 
		f.write('  QUANTITY z' + str(i) + ': REAL;\n')
		f.write('  QUANTITY z' + str(i) + 'dot: REAL;\n')
	f.write('BEGIN\n')
	f.write('	IF (domain = quiescent_domain) USE\n')
	for i in range(1, nz + 1): 
		f.write('    z' + str(i) + ' == 0.0;\n')
		f.write('    z' + str(i) + 'dot == 0.0;\n')
	f.write('	  i == 0.0;\n')
	f.write('	ELSE\n')
	for i in range(1, nz + 1): 
		f.write('    z' + str(i) + 'dot == z' + str(i) + '\'dot;\n')
	f.write('	  q == i\'integ*Ccon;\n')
	for i in range(nz):
		f.write('    ')
		for j in range(nz):
			f.write(repr(MatrixM(sys)[i, j]) + '*z' + str(j + 1) + 'dot\'dot +')
		for j in range(nz):
			f.write(repr(MatrixE(sys)[i, j]) + '*z' + str(j + 1) + 'dot +')
		for j in range(nz):
			f.write(repr(MatrixK(sys)[i, j]) + '*z' + str(j + 1))
			if (j != nz - 1):
				f.write(' + ')
			else:
				f.write(' == ')
		f.write(repr(MatrixB(sys)[i, 0]) + '*f + ' + repr(MatrixB(sys)[i, 1]) + '*q;\n')
	f.write('	END USE;\n')
	f.write('  pos*Scon == ')
	for j in range(nz):
		f.write(repr(MatrixC(sys)[0, j]) + '*z' + str(j + 1))
		if (j != nz - 1):
			f.write(' + ')
	f.write(';\n')
	f.write('  v*Vcon == ')
	for j in range(nz):
		f.write(repr(MatrixC(sys)[1, j]) + '*z' + str(j + 1))
		if (j != nz - 1):
			f.write(' + ')
	f.write(';\n')	
	f.write('END ARCHITECTURE behav;\n')
	f.close()

def WriteSystem(sys, basename, Format = 'MatrixMarket', Circuit = 'circuit'):
	"""writes a system to files with the base name baseName."""
	if Format == 'Spice':
		writeSpice(sys, basename, Circuit)
		return
	if Format == 'Simplorer':
		writeSimplorer(sys, basename)
		return
	if Format == 'VHDL-AMS':
		writeVHDLAMS(sys, basename)
		return
	if MatrixMQ(sys):
		writeMatrix(basename + '.M', MatrixM(sys), Format)
	if MatrixEQ(sys):
		writeMatrix(basename + '.E', MatrixE(sys), Format)
	if FirstOrderSystemQ(sys):
		writeMatrix(basename + '.A', MatrixA(sys), Format)
	else:
		writeMatrix(basename + '.K', MatrixK(sys), Format)
	writeMatrix(basename + '.B', MatrixB(sys), Format)
	writeMatrix(basename + '.C', MatrixC(sys), Format)
	writeNames(basename + '.C.names', OutputNames(sys))
	if MatrixDQ(sys):
		writeMatrix(basename + '.D', MatrixD(sys), Format)

def checkConsistenceR(res):
	dim = YSeries(res).shape
	if len(res) != 4:
		raise BaseException('SimulationResult: The length of the list is not 4')
	elif len(XSeries(res)) != dim[1]:
		raise BaseException('SimulationResult: XSeries and YSeries are not compatible')
	elif len(YNames(res)) != dim[0]:
		raise BaseException('SimulationResult: YNames and YSeries are not compatible')
	return res

def MakeSimulationResult(xseries, xname, yseries, ynames):
	"""gives a SimulationResult object. 
	
	XSeries must be a vector, YSeries must be a matrix, XName must be a
	string and YNames must be a list of strings."""
	return checkConsistenceR([numpy.array(xseries), xname, numpy.array(yseries), ynames])

def XSeries(res):
	"""gives XSeries of res."""
	return res[0]

def XName(res):
	"""gives XName of res."""
	return res[1]

def YSeries(res):
	"""gives YSeries of res."""
	return res[2]

def YNames(res):
	"""gives YNames of res."""
	return res[3]

def TransformResult(res, names):
	"""gives a new res whose YSeries correspond to names."""
	seri = []
	nam = []
	for ni in names:
		pos = YNames(res).index(ni)
		nam.append(ni)
		seri.append(pos)
	return MakeSimulationResult(
		XSeries(res),
		XName(res),
		YSeries(res)[seri],
		nam
	)

def CompatibleResultQ(*arg):
	"""CompatibleResultQ(res1, res2, ...) 
	
	gives True if all objects are compatible and False otherwise."""
	if len(arg) == 1:
		return True
	elif len(arg) == 2:
		res1 = arg[0]
		res2 = arg[1]
		return XName(res1) == XName(res2) and YNames(res1) == YNames(res2)
	else:
		ret = True
		for i in arg[1:]:
			if not CompatibleResultQ(arg[0], i):
				ret = False
				break
		return ret

def diff1(res1, i, ErrorFunction):
	if (XSeries(res1) != XSeries(i)).all():
		raise BaseException('Difference: Different XSeries are not implemented')
	return numpy.array(map(lambda x, y: map(ErrorFunction, x, y), YSeries(res1), YSeries(i)))

def Difference(res1, resList, ErrorFunction = lambda x, y: x - y):
	"""gives a list of differences between res1 and each object in listres. 
	
	All SimulationResult must be compatible. Option is ErrorFunction, by
	default it is the difference."""
	if not CompatibleResultQ(res1, *resList):
		raise BaseException('Difference: Results are not compatible')
	return map(
		lambda i:
			MakeSimulationResult(
				XSeries(res1),
				XName(res1),
				diff1(res1, i, ErrorFunction),
				YNames(res1)
			)
		, resList)

def TakeResult(res, num):
	"""takes first num values for results."""
	return MakeSimulationResult(
			XSeries(res)[:num],
			XName(res),
			YSeries(res)[:, :num],
			YNames(res)
		)

def ReadResult(file):
	"""reads SimulationResult from file in Matrix Market format."""
	mat = readMatrix1(file)
	name = readNames(file + '.names')
	return MakeSimulationResult(
			mat[0],
			name[0],
			mat[1:],
			name[1:]
		)

def WriteResult(res, file):
	"""writes res to file in Matrix Market format."""
	writeMatrix(file, numpy.insert(YSeries(res), 0, XSeries(res), 0), 'MatrixMarket')
	nam = YNames(res)[:]
	nam.insert(0, XName(res))
	writeNames(file + '.names', nam)

def ReadResultTable(file):
	"""reads SimulationResult from file in Table format."""
	mat = readMatrix1(file, False, Format='Caspoc')
	mat = mat.transpose()
	name = readNames(file + '.names')
	return MakeSimulationResult(
			mat[0],
			name[0],
			mat[1:],
			name[1:]
		)

def WriteResultTable(res, file):
	"""writes res to file as a table."""
	scipy.io.write_array(file, numpy.insert(YSeries(res), 0, XSeries(res), 0).transpose(), precision=16)
	nam = YNames(res)[:]
	nam.insert(0, XName(res))
	writeNames(file + '.names', nam)

def WriteResultCSV(res, file):
	"""writes res to file in the CSV format."""
	nam = YNames(res)[:]
	nam.insert(0, XName(res))
	f = open(file, 'w')
	for n in nam[:-1]:
		f.write(n + '[],')
	f.write(nam[-1] + '[]\n')
	for i in range(len(XSeries(res))):
		f.write(repr(XSeries(res)[i]) + ',')
		for j in YSeries(res)[:-1, i]:
			f.write(repr(j) + ',')
		f.write(repr(YSeries(res)[-1, i]) + '\n')
	f.close()

def ReadResultOld(name, nameX, PhaseAngle=False):
	"""reads SimulationResult from file in my old format.
	
	nameX should be 'Time' or 'Frequency'."""
	f = open(name)
	tokens = f.read().split()
	f.close()
	ntime = int(float(tokens[0]))
	time = numpy.array([float(i) for i in tokens[1:ntime + 1]])
	itok = 1 + ntime
	names = []
	series = []
	step = 1
	if PhaseAngle:
		step = 2
	while itok < len(tokens):
		names.append(tokens[itok])
		series.append([float(i) for i in tokens[itok + 1:itok + step*ntime + 1:step]])
		itok += step*ntime + 1
	return MakeSimulationResult(time, nameX, numpy.array(series), names)

def ReadResultGrace(file):
	"""reads SimulationResult from file in Grace format."""
	f = open(file)
	lines = f.readlines()
	f.close()
	names = lines[0].split()
	if names[1] != 'FREQ':
		raise BaseException('ReadResultGrace: wrong first line - no FREQ')
	names = names[2:]
	n = len(names)
	what = lines[1].split()
	if len(what) != n*2 + 1:
		raise BaseException('ReadResultGrace: wrong second line - its length is inconsistent with the first one')
	if what[1] == 'AMPLITUDE' and what[2] == 'PHASE':
		convert = True
	elif what[1] == 'REAL' and what[2] == 'IMAGINARY':
		convert = False
	else:
		raise BaseException('ReadResultGrace: wrong second line - neither AMPLITUDE+PHASE nor REAL+IMAGINARY')
	m = 0
	for line in lines[2:]:
		if len(line.split()) != 2*n + 1:
			break
		else:
			m += 1
	freq = numpy.zeros((m))
	mat = numpy.zeros((n, m), dtype=complex)
	j = 0
	for line in lines[2:2+m]:
		val = line.split()
		freq[j] = float(val[0])
		if convert:
			for i in range(n):
				angle = float(val[2 + 2*i])*degree
				amp = float(val[1 + 2*i])
				mat[i, j] = complex(amp*math.cos(angle), amp*math.sin(angle))
		else:
			for i in range(n):
				mat[i, j] = complex(float(val[1 + 2*i]), float(val[2 + 2*i]))
		j += 1
	return MakeSimulationResult(
			freq,
			'Frequency',
			mat,
			names
		)
	
def WriteResultGrace(res, file, PhaseAngle = False):
	if XName(res) != 'Frequency':
		raise BaseException('WriteResultGrace: not a file for harmonic simulation')
	of = open(file, 'w')
	of.write('#  FREQ')
	n = len(YNames(res))
	for name in YNames(res):
		of.write('%28s' % name)
	of.write('\n')
	of.write('#              ')
	if PhaseAngle:
		for i in range(n):
			of.write('       AMPLITUDE       PHASE')
	else:
		for i in range(n):
			of.write('        REAL       IMAGINARY')
	of.write('\n')
	freq = XSeries(res)
	mat = YSeries(res).transpose()
	for i in range(len(XSeries(res))):
		of.write('%9.3f' % freq[i])
		of.write('      ')
		if PhaseAngle:
			for j in range(n):
				val = complex(mat[i, j])
				amp = abs(val)
				angle = math.atan2(val.imag, val.real)/degree
				of.write('%14G%14G'%(amp, angle))
		else:
			for j in range(n):
				val = complex(mat[i, j])
				of.write('%14G%14G'%(val.real, val.imag))
		of.write('\n')
	of.close()

def WriteResultAsDirichlet(res, file):
	"""writes res to file as ANSYS script to fix values as Dirichlet constraints."""
	of = open(file, 'w')
#	of.write("resume\n")
	of.write("/solu\n")
	of.write("allsel\n")
	of.write("antype,3\n")
	ser = YSeries(res).transpose()
	if numpy.iscomplexobj(ser):
		fun = lambda v: repr(v.real) + ',' + repr(v.imag)
	else:
		fun = repr
	for i in range(len(XSeries(res))):
		of.write("harfrq," + repr(XSeries(res)[i]) + "\n")
		for j in range(len(YNames(res))):
			val = fun(ser[i, j])
			s = YNames(res)[j].replace('Node_', '')
			s = s.replace('_', ',')
			of.write("D," + s + "," + val + "\n")
		of.write("solve\n")
	of.write("fini\n")		
	of.close()

def PlotResult(listres, 
		FileName = 'plot', 
		CommonTitle = '', 
		FunctionX = lambda x: x, 
		FunctionY = lambda x: x,
		LogX = False,
		LogY = False,
		LimX = False,
		LimY = False,
		Options = [{'color':(1,0,0)}, {'color':(0,1,0)}, {'color':(0,0,1)}],
		Grid = False,
		LabelX = False,
		LabelY = False,
		Legend = False):
	"""makes plots for the list of compatible SimulationResult. 
	
	Options: FileName for the base name, CommonTitle for a title,
	FunctionX and FunctionY to modify series values, LogX and LogY to
	turn on a logarithmic scale for X and Y axes, Options is a list of
	dictionaries with matplotlib properties for each series."""
	if not CompatibleResultQ(*listres):
		raise BaseException('PlotResult: Results are not compatible')
	for i in range(len(YNames(listres[0]))):
		pylab.title(YNames(listres[0])[i] + ': ' + CommonTitle)
		ax = pylab.gca()
		if LogX:
			ax.set_xscale("log")
		if LogY:
			ax.set_yscale("log")
		for j in range(len(listres)):
			x = [FunctionX(k) for k in XSeries(listres[j])]
			y = [FunctionY(k) for k in YSeries(listres[j])[i]]
			pylab.plot(x, y, **Options[j % len(Options)])
		pylab.grid(Grid)
		if LimX:
			pylab.xlim(LimX)
		if LimY:
			pylab.ylim(LimY)
		if LabelX:
			pylab.xlabel(LabelX)
		if LabelY:
			pylab.ylabel(LabelY)
		if Legend:
			pylab.legend(Legend)
		pylab.savefig(FileName + str(i + 1), dpi=125)
		pylab.clf()

def PlotResultTogether(listres, 
		FileName = 'plot', 
		CommonTitle = '', 
		FunctionX = lambda x: x, 
		FunctionY = lambda x: x,
		LogX = False,
		LogY = False,
		LimX = False,
		LimY = False,
		Options = [{'color':(1,0,0)}, {'color':(0,1,0)}, {'color':(0,0,1)}],
		Grid = False,
		LabelX = False,
		LabelY = False):
	"""makes plots for SimulationResult in the same figure. 
	
	Options: FileName for the base name, CommonTitle for a title,
	FunctionX and FunctionY to modify series values, LogX and LogY to
	turn on a logarithmic scale for X and Y axes, Options is a list of
	dictionaries with matplotlib properties for each series."""
	Legend = []
	pylab.title(CommonTitle)
	ax = pylab.gca()
	if LogX:
		ax.set_xscale("log")
	if LogY:
		ax.set_yscale("log")
	pylab.grid(Grid)
	if LimX:
		pylab.xlim(LimX)
	if LimY:
		pylab.ylim(LimY)
	if LabelX:
		pylab.xlabel(LabelX)
	if LabelY:
		pylab.ylabel(LabelY)
	kk = 0
	for j in range(len(listres)):
		for i in range(len(YNames(listres[j]))):
			x = [FunctionX(k) for k in XSeries(listres[j])]
			y = [FunctionY(k) for k in YSeries(listres[j])[i]]
			pylab.plot(x, y, **Options[kk])
			kk = kk + 1
			if kk == len(Options):
				kk = 0
			Legend.append(YNames(listres[j])[i])
	pylab.legend(Legend)
	pylab.savefig(FileName + str(i + 1), dpi=125)
	pylab.clf()

def ComputePoles(sys): 
	s = sys
	if SecondOrderSystemQ(sys):
		s = ToFirstOrderSystem(sys)
	return scipy.linalg.eigvals(MatrixA(s), MatrixE(s))

def PlotPoles(eig, 
		FileName = 'poles', 
		LimX = False,
		LimY = False,
		Options = [{'marker':'x', 'markerfacecolor':(0,1,0), 'linewidth':0}],
		Grid = False,
		LabelX = False,
		LabelY = False):
	ax = pylab.gca()
	pylab.grid(Grid)
	if LimX:
		pylab.xlim(LimX)
	if LimY:
		pylab.ylim(LimY)
	if LabelX:
		pylab.xlabel(LabelX)
	if LabelY:
		pylab.ylabel(LabelY)
	pylab.plot(eig.real, eig.imag, **Options[0])
	pylab.savefig(FileName, dpi=125)
	pylab.clf()	

def NumericQ(list):
	ret = True
	for i in list:
		if not operator.isNumberType(i):
			ret = False
			break
	return ret

def makeInputs(inpFun, n):
	if inpFun == 'Automatic':
		inp = [1. for i in range(n)]
	elif type(inpFun) == type([]) and len(inpFun) == n:
		inp = inpFun
	elif n == 1 and type(inpFun) != type([]):
		inp = [inpFun]
	else:
		raise BaseException('makeInputs: Input functions are not compatible with the system')
	if NumericQ(inp):
		return inp
#	state = [0. for i in range(n)]
	for i in range(n):
		if callable(inp[i]):
			if not operator.isNumberType(inp[i](1.)):
				raise BaseException('makeInputs: Not numeric values of input function')
		else:
			if operator.isNumberType(inp[i]):
				val = inp[i]
				inp[i] = lambda tau: val
			else:
				raise BaseException('makeInputs: Not numeric values with input functions')
	return inp

def StationarySolution(sys, InputFunction = 'Automatic', ReturnState = False):
	"""gives a stationary solution. 
	
	InputFunction is optionally a list of values to scale inputs."""
	inpFun = makeInputs(InputFunction, MatrixB(sys).shape[1])
	if not NumericQ(inpFun):
		raise BaseException('Harmonic Solution: Input functions must be numeric vales')
	state = scipy.linalg.solve(MatrixK(sys), numpy.dot(MatrixB(sys), inpFun))
	if ReturnState:
		return state
	else:
		return numpy.dot(MatrixC(sys), state)

def harmonicLoop(list, fun, rhs):
	return map(lambda w: scipy.linalg.solve(fun(w), rhs(w)), 2.*numpy.pi*numpy.array(list))

def checkFrequency(freq):
	if not NumericQ(freq):
		raise BaseException('Harmonic Solution: Not numeric values among frequencies')

def writeState(SaveState, freq, res, Title):
	if len(res) != len(Title[1]):
		raise BaseException('Harmonic Solution: SaveState: the number of frequencies in Title is not equal to number of computer frequencies')
	res = numpy.array(res).transpose()
	f = open(SaveState + '.st', 'w')
	f.write(str(res.shape[0]) + ' ' + str(res.shape[1]) + '\n')
	f.write(Title[0] + '"' + Title[2] + '"\n')
	for i in Title[1]:
		f.write(i + '"')
	f.write('\n')
	for i in freq:
		f.write(str(i) + ' ')
	f.write('\n')
	shape = list(res.shape)
	shape[1] = shape[1]*2
	r = numpy.zeros(shape)
	for i in range(res.shape[0]):
		for j in range(res.shape[1]):
			v = complex(res[i, j])
			r[i, 2*j] = v.real
			r[i, 2*j + 1] = v.imag
	scipy.io.write_array(f, r, precision=17)
	f.close()

def HarmonicSolution(freq, sys, InputFunction = 'Automatic', SaveState = False, Title = False):
	"""gives SimulationResult for a given list of freq.

	InputFunction is optionally a list of values or functions in w to
	modify inputs."""
	checkFrequency(freq)
	inpFun = makeInputs(InputFunction, MatrixB(sys).shape[1])
	if not NumericQ(inpFun):
		rhs = lambda w: numpy.dot(MatrixB(sys), numpy.array([i(w) for i in inpFun]))
	else:
		ret = numpy.dot(MatrixB(sys), numpy.array(inpFun))
		rhs = lambda w: ret
	if not MatrixMQ(sys):
		res = harmonicLoop(freq, lambda w : 1j*w*MatrixE(sys) + MatrixK(sys), rhs)
	elif not MatrixEQ(sys):
		res = harmonicLoop(freq, lambda w : -w**2*MatrixM(sys) + MatrixK(sys), rhs)
	else:
		res = harmonicLoop(freq, lambda w : -w**2*MatrixM(sys) + 1j*w*MatrixE(sys) + MatrixK(sys), rhs)
	if (SaveState != False):
		writeState(SaveState, freq, res, Title)
	return MakeSimulationResult(
			freq, 
			"Frequency", 
			numpy.dot(MatrixC(sys), numpy.array(res).transpose()),
			OutputNames(sys)
		)

def HarmonicSolutionTest(freq, sys, InputFunction = 'Automatic', SaveState = False, Title = False):
	"""gives SimulationResult for a given list of freq.

	InputFunction is optionally a list of values or functions in w to
	modify inputs."""
	checkFrequency(freq)
	inpFun = makeInputs(InputFunction, MatrixB(sys).shape[1])
	if not NumericQ(inpFun):
		rhs = lambda w: numpy.dot(MatrixB(sys), numpy.array([i(w) for i in inpFun]))
	else:
		ret = numpy.dot(MatrixB(sys), numpy.array(inpFun))
		rhs = lambda w: ret*w*w
	if not MatrixMQ(sys):
		res = harmonicLoop(freq, lambda w : 1j*w*MatrixE(sys) + MatrixK(sys), rhs)
	elif not MatrixEQ(sys):
		res = harmonicLoop(freq, lambda w : -w**2*MatrixM(sys) + MatrixK(sys), rhs)
	else:
		res = harmonicLoop(freq, lambda w : -w**2*MatrixM(sys) + 1j*w*w*MatrixE(sys) + MatrixK(sys), rhs)
	if (SaveState != False):
		writeState(SaveState, freq, res, Title)
	return MakeSimulationResult(
			freq, 
			"Frequency", 
			numpy.dot(MatrixC(sys), numpy.array(res).transpose()),
			OutputNames(sys)
		)

def FrequencyConvergence(freq, sys, start = 1, finish = numpy.inf, step = 1):
	"""gives an object for LocalError and LocalErrorIndicator.
	
	This object is computed for sys for frequency freq. start, finish,
	step allows us optionally to specify the dimension range."""
	fin = finish
	if fin > len(MatrixK(sys)):
		fin = len(MatrixK(sys))
	return [[i, YSeries(HarmonicSolution([freq], TakeSystem(sys, i))).transpose()[0]] for i in range(start, fin + 1, step)]

def RelativeError(x, y):
	"""gives a relative error for two numbers of two vectors."""
	return scipy.linalg.norm(x - y)/scipy.linalg.norm(x)

def LocalError(list, exact, ErrorFunction = RelativeError):
	"""gives a list of local errors.
	
	It is computed from the object obtained by FrequencyConvergence and
	a list of true values for the same frequency. Option is
	ErrorFunction (by default RelativeError)."""
	return [[i[0],  ErrorFunction(exact, i[1])] for i in list]

def LocalErrorIndicator(list, ErrorFunction = RelativeError):
	"""gives a list of local error estimates.
	
	It is computed from the object conv obtained by
	FrequencyConvergence. Option is ErrorFunction (by default
	RelativeError)."""
	return [[list[i][0],  ErrorFunction(list[i + 1][1], list[i][1])] for i in range(len(list) - 1)]

def PlotConvergence(errlist, 
		FileName = 'err', 
		Title = 'Relative Error vs. Dimension', 
		LogY = True,
		Options = [{'marker':'d', 'markerfacecolor':(1,0,0), 'linewidth':0}, {'marker':'+', 'markerfacecolor':(0,1,0), 'linewidth':0}, {'marker':'o', 'markerfacecolor':(0,0,1), 'linewidth':0}],
		LabelX = False,
		LabelY = False,
		Legend = False):
	"""plots a list of objects obtained by LocalError and LocalErrorIndicator.

	Options: FileName for the base file name, Title for plot title, LogY
	to turn on logarithmic scale, Options is a list of dictionaries with 
	matplotlib properties for each series, Legend is a list of strings
	for the legend.""" 
	pylab.title(Title)
	ax = pylab.gca()
	if LogY:
		ax.set_yscale("log")
	def convert2plot(list):
		return numpy.array(list).transpose()
	j = 0
	for err in errlist:
		x, y = convert2plot(err)
		pylab.plot(x, y, **Options[j % len(Options)])
		j += 1
	if LabelX:
		pylab.xlabel(LabelX)
	if LabelY:
		pylab.ylabel(LabelY)
	if Legend:
		pylab.legend(Legend)
	pylab.savefig(FileName, dpi=125)
	pylab.clf()

def AnsysTransientSolution(time, sys, InputFunction='Automatic', InitialState='Automatic', TINTP=[0.005,None,None,1], Velocity=False, SaveState = False):
	"""gives SimulationResult for a given list time.
	
	It uses integrators like in ANSYS: backward Euler for a first order
	system and Newmark for a second order system. InputFunction is
	optionally a list of values or functions in time to modify inputs.
	TINTP is a list of [gamma, alpha, delta, theta] (see ANSYS manual).
	InitialState is an optional vector of initial states, by default it
	is zero."""
	if FirstOrderSystemQ(sys):
		return firstAnsysTransientSolution(time, sys, InputFunction, InitialState, TINTP, Velocity, SaveState)
	else:
		return secondAnsysTransientSolution(time, sys, InputFunction, InitialState, TINTP, Velocity)

def writeStateTr(file, timeSaved, time, res):
	res = res.transpose()
	f = open(file, 'w')
	f.write(str(len(timeSaved)) + ' ' + str(res.shape[1]) + '\n')
	for i in timeSaved:
		f.write(str(i) + ' ')
	f.write('\n')
	seri = []
	for ti in timeSaved:
		pos = time.index(ti)
		seri.append(pos)
	scipy.io.write_array(f, res[seri], precision=17)
	f.close()

def firstAnsysTransientSolution(time, sys, InputFunction, InitialState, TINTP, Velocity, SaveState):
	theta = TINTP[3]
	init = InitialState
	if init == 'Automatic':
		init = numpy.zeros((len(MatrixK(sys))))
	inpFun = makeInputs(InputFunction, MatrixB(sys).shape[1])
	if not NumericQ(inpFun):
		rhs = lambda tau: numpy.dot(MatrixB(sys), numpy.array([i(tau) for i in inpFun]))
	else:
		ret = numpy.dot(MatrixB(sys), numpy.array(inpFun))
		rhs = lambda tau: ret
#	0 tprev = 0.
#	1 dtprev = 0.
#	2 un = init
#	3 vn = numpy.zeros((len(MatrixK(sys))))
	def fun(tau, par = [0., 0., init, numpy.zeros((len(MatrixK(sys)))), 0]):
		dt = tau - par[0]
		par[0] = tau
		if dt != 0:
			if dt != par[1]:
				par[4] = scipy.linalg.lu_factor(1./theta/dt*MatrixE(sys) + MatrixK(sys))
			un1 = scipy.linalg.lu_solve(par[4], rhs(tau) + numpy.dot(MatrixE(sys), 1./theta/dt*par[2] + (1 - theta)/theta*par[3]))
			par[3] = 1./theta/dt*(un1 - par[2] - (1 - theta)*dt*par[3])
			par[1] = dt
			par[2] = un1
		return par[2]
	res = numpy.array(map(fun, time)).transpose()
	if SaveState != False:
		writeStateTr(SaveState[0], SaveState[1], time, res)
	return MakeSimulationResult(
		time,
		'Time',
		numpy.dot(MatrixC(sys), res), 
		OutputNames(sys)
	)

def secondAnsysTransientSolution(time, sys, InputFunction, InitialState, TINTP, Velocity):
	if TINTP[0]	!= None:
		gamma = TINTP[0]
		alpha = 0.25*(1 + gamma)**2
		delta = 0.5 + gamma
	else:
		alpha = TINTP[1]
		delta = TINTP[2]
	init = InitialState
	if init == 'Automatic':
		init = numpy.zeros((len(MatrixK(sys))))
	inpFun = makeInputs(InputFunction, MatrixB(sys).shape[1])
	if not NumericQ(inpFun):
		rhs = lambda tau: numpy.dot(MatrixB(sys), numpy.array([i(tau) for i in inpFun]))
	else:
		ret = numpy.dot(MatrixB(sys), numpy.array(inpFun))
		rhs = lambda tau: ret
#0	tprev = 0
#1	dtprev = 0
#2	un = init
#3  vn = numpy.zeros((len(MatrixK(sys))))
#4	an = numpy.zeros((len(MatrixK(sys))))
	def fun(tau, par = [0., 0., init, numpy.zeros((len(MatrixK(sys)))), numpy.zeros((len(MatrixK(sys)))), 0]):
		dt = tau - par[0]
		par[0] = tau
		if dt != 0.:
			a0 = 1./alpha/dt**2
			a1 = delta/alpha/dt
			a2 = 1./alpha/dt
			a3 = 1./2/alpha - 1
			a4 = delta/alpha - 1.
			a5 = dt/2*(delta/alpha - 2)
			a6 = dt*(1 - delta)
			a7 = delta*dt
			if dt != par[1]:
				par[5] = scipy.linalg.lu_factor(a0*MatrixM(sys) + a1*MatrixE(sys) + MatrixK(sys))
			un1 = scipy.linalg.lu_solve(par[5], rhs(tau) + numpy.dot(MatrixM(sys), a0*par[2] + a2*par[3] + a3*par[4]) + numpy.dot(MatrixE(sys), a1*par[2] + a4*par[3] + a5*par[4]))
			an1 = a0*(un1 - par[2]) - a2*par[3] - a3*par[4]
			par[3] = par[3] + a6*par[4] + a7*an1
			par[4] = an1
			par[1] = dt
			par[2] = un1
# must be par[2], par[3] for velocities
		if Velocity:
			return par[3]
		else:
			return par[2]
	res = numpy.array(map(fun, time)).transpose()
	return MakeSimulationResult(
		time,
		'Time',
		numpy.dot(MatrixC(sys), res), 
		OutputNames(sys)
	)

def ExpandMOR(fileV, fileState, fileAnsys, Dirichlet):
	V = scipy.io.mmread(fileV)
	names = readNames(fileV + '.names')
	if len(names) != len(V):
		raise BaseException('Names and V are not compatible')
	f = open(fileState)
	n, m = f.readline().split()
	n = int(n)
	m = int(m)
	times = f.readline().split()
	st = scipy.io.read_array(f)
	f.close()
	if len(times) != n:
		raise BaseException('Times are not compatible with n')
	if len(Dirichlet) != n:
		raise BaseException('Times are not compatible with n')
	if st.shape[0] != n:
		raise BaseException('Shapes are not compatible with n')
	if st.shape[1] != m:
		raise BaseException('Shapes are not compatible with m')
	if m != V.shape[1]:
		raise BaseException('m is not compatible with V')
	res = numpy.dot(V, st.transpose())
	of = open(fileAnsys, 'w')
#	of.write("resume\n")
	of.write("/solu\n")
	of.write("allsel\n")
	of.write("antype,4\n")
	for i in range(n):
		of.write("time," + times[i] + "\n")
		of.write("D,all,TEMP," + repr(Dirichlet[i]) + "\n")
		for j in range(len(names)):
			val = repr(res[j, i])
			s = names[j].replace('Node_', '')
			s = s.replace('_', ',')
			of.write("D," + s + "," + val + "\n")
		of.write("solve\n")
	of.write("fini\n")		
	of.close()