(* 
	Evgenii Rudnyi (c) 2004, 2005, 2006, http://Evgenii.Rudnyi.Ru/
*)

(* This program is free software; *)
(* you can redistribute it and or modify it under the terms of the GNU General Public License *)
(* as published by the Free Software Foundation; either version 2 of the License,*)
(* or (at your option) any later version.This program is distributed in the \
hope that *)
(* it will be useful, but WITHOUT ANY WARRANTY; *)
(* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *)
(* See the GNU General Public License for more details. 
      You should have received a copy of *)
(* the GNU General Public License along with this program; if not, 
  write to the *)
(* Free Software Foundation,Inc.,59 Temple Place,Suite 330,Boston,
  MA 02111-1307 USA *)

BeginPackage["Imtek`Post4MOR`", {"LinearAlgebra`MatrixManipulation`", "Utilities`FilterOptions`", "Graphics`MultipleListPlot`"}];

DynamicSystem::usage = "Use MakeDynamicSystem, MakeFirstOrderSystem, and MakeSecondOrderSystem to construct the object; MatrixMQ, MatrixEQ, MatrixDQ, FirstOrderSystemQ, SecondOrderSystemQ to make a query; MatrixM, MatrixE, MatrixK, MatrixA, MatrixB, MatrixC, MatrixD, OutputNames to select a component.";

MakeDynamicSystem::usage = "MakeDynamicSystem[matM, matE, matK, matB, matC, namesC, matD] gives a DynamicSystem object. matM, matE, and matD should be matrices or Null. matK, matB, and matC must be matrices. namesC is a list of strings. matD is optional and by default is Null. MakeDynamicSystem[matM, matE, matK, matB, matD] makes matC equal to the identity matrix.";
MakeFirstOrderSystem::usage = "MakeFirstOrderSystem[matA, matB, matC, namesC, matD] and MakeFirstOrderSystem[matE, matA, matB, matC, namesC, matD] give a DynamicSystem object for a first order system. matE, and matD should be matrices or Null. matA, matB, and matC must be matrices. namesC is a list of strings. matD is optional and by default is Null. MakeFirstOrderSystem[matE, matA, matB, matD] makes matC equal to the identity matrix.";
MakeSecondOrderSystem::usage = "MakeSecondOrderSystem[matM, matK, matB, matC, namesC, matD] and MakeSecondOrderSystem[matM, matE, matK, matB, matC, namesC, matD] give a DynamicSystem object for a second order system. matM, matE, and matD should be matrices or Null. matK, matB, and matC must be matrices. namesC is a list of strings. matD is optional and by default is Null. MakeSecondOrderSystem[matM, matE, matK, matB, matD] makes matC equal to the identity matrix.";

MatrixM::usage = "MatrixM[sys] gives system matrix M of sys.";
MatrixE::usage = "MatrixE[sys] gives system matrix E of sys.";
MatrixA::usage = "MatrixA[sys] gives system matrix A of sys.";
MatrixK::usage = "MatrixK[sys] gives system matrix K of sys.";
MatrixB::usage = "MatrixB[sys] gives system matrix B of sys.";
MatrixC::usage = "MatrixC[sys] gives system matrix C of sys.";
MatrixD::usage = "MatrixD[sys] gives system matrix D of sys.";
OutputNames::usage = "OutputNames[sys] gives a list of output names of sys.";

MatrixMQ::usage = "MatrixMQ[sys] gives True if system matrix M is present in sys and False otherwise.";
MatrixEQ::usage = "MatrixEQ[sys] gives True if system matrix E is present in sys and False otherwise.";
MatrixDQ::usage = "MatrixMD[sys] gives True if system matrix D is present in sys and False otherwise.";
FirstOrderSystemQ::usage = "FirstOrderSystemQ[sys] gives True if sys represent a first order system and False otherwise.";
SecondOrderSystemQ::usage = "SecondOrderSystemQ[sys] gives True if sys represent a second order system and False otherwise.";
ExplicitSystemQ::usage = "ExplicitSystemQ[sys] gives True if sys represent an explicit first order system and False otherwise.";

DeleteOutputs::usage = "DeleteOutputs[sys, names] gives a new system when outputs defined in a list of names are deleted.";
AddOutputs::usage = "AddOutputs[sys, matCadd, names] gives a new system with additional outputs defined by rows of matCadd and a list of names.";
DeleteDamping::usage = "DeleteDamping[sys] gives a new system without system matrix E. sys must be a second order system.";
AddDamping::usage = "AddDamping[sys, matE] gives a new system with matrix E inreased by matE. AddDamping[sys, alpha, beta] gives a new system with matrix E increased by Rayleigh damping alpha*matM+beta*matK. sys must be a second order system.";
AddDampingRatio::usage = "AddDampingRatio[sys, imatK] gives a new system with matrix K inreased by I imatK. AddDampingRatio[sys, dmprat] gives a new system with matrix K increased by I 2 dmprat MatrixK. sys must be a second order system.";
ToSingleInputSystem::usage = "ToSingleInputSystem[sys] gives a new system in which there is only first input left.";
ToFirstOrderSystem::usage = "ToFirstOrderSystem[sys] gives a new system in the form of a first order system.";
ToExplicitSystem::usage = "ToExplicitSystem[sys] gives a new system in the form of an explicit first order system.";

TakeSystem::usage = "TakeSystem[sys, dim] gives a new system of dimension dim. This operation makes sense only when an original system has been obtained by iterative model reduction method.";
ReadSystem::usage = "ReadSystem[baseName] reads a system from files with the base name baseName. ReadSystem[file, type] reads a system from a file written in the old format of mor4ansys 1.5. type must be a string FirstOrderSystem or SecondOrderSystem.";
WriteSystem::usage = "WriteSystem[sys, baseName] writes a system to files with the base name baseName.";

SimulationResult::usage = "Use MakeSimulationResult to construct the object; XSeries, XName, YSeries, YNames to select a component.";

MakeSimulationResult::usage = "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.";

XSeries::usage = "XSeries[res] gives XSeries of res.";
XName::usage = "XName[res] gives XName of res.";
YSeries::usage = "YSeries[res] gives YSeries of res.";
YNames::usage = "YNames[res] gives YNames of res.";

TransformResult::usage = "TransformResult[res, names] gives a new res whose YSeries correspond to names. res also can be a list of SimulationResult. Option TransformFunction to define fun[y, x].";
TransformFunction::usage = "An option for TransformResult to define a function fun[y, x]. If it is a single function, it is applied to all YSeries. Alternatively, one can specify a list of pairs {name, fun}.";

CompatibleResultQ::usage = "CompatibleResultQ[res1, res2, ...] gives True if all objects are compatible and False otherwise.";

Difference::usage = "Difference[res1, {listres}, ops] gives a list of differences between res1 and each object in listres. All SimulationResult must be compatible. Option is ErrorFunction, by default ErrorFunction->Subtract.";

ErrorFunction::usage = "An option for Difference, LocalError and LocalErrroIndicator.";

TakeResult::usage = "TakeResult[res, num] takes first num values for results.";

ReadResult::usage = "ReadResult[file] reads SimulationResult from file in Matrix Market format. ReadResult[file, XName] reads SimulationResult from file in my old format, where XName should be \"Time\" or \"Frequency\".";
WriteResult::usage = "WriteResult[res, file] writes res to file in Matrix Market format. WriteResult[res, file, \"Table\"] writes res as a table.";

WriteResultOld::usage = "WriteResultOld[res, file] writes res to file in my old format.";
WriteResultAsDirichlet::usage = "WriteResultAsDirichlet[res, i, file] writes the i-th set of res to file as ANSYS script to fix values as Dirichlet constraints.";

Log10::usage = "Log10[x] gives Log[10, x].";

Log10RelativeError::usage = "Log10RelativeError[x, y] gives Log10[Norm[x - y]/Norm[x]].";

PlotResult::usage = "PlotResult[listres, ops] makes plots for the list of compatible SimulationResult. Options: CommonTitle for a title, FunctionX and FunctionY to modify series values, MultipleListPlot options.";

CommonTitle::usage = "An option for PlotResult.";
FunctionX::usage = "An option for PlotResult.";
FunctionY::usage = "An option for PlotResult.";

StationarySolution::usage = "StationorySolution[sys, ops] gives a stationary solution. Options: InputFunction, LinearSolve options.";
HarmonicSolution::usage = "HarmonicSolution[freq, sys, ops] gives SimulationResult for a given list of freq. Options: InputFunction, LinearSolve options.";
FrequencyConvergence::usage = "FrequencyConvergence[freq, sys, start, finish] gives an object for LocalError and LocalErrorIndicator computed for sys with dimensions from start to finish for frequency freq. FrequencyConvergence[freq, sys, start, finish, step] allows us to specify step.";
LocalError::usage = "LocalError[conv, true, ops] gives a list of local errors computed from the object conv obtained by FrequencyConvergence and a list of true values for the same frequency. Option is ErrorFunction (by default Log10RelativeError).";
LocalErrorIndicator::usage = "LocalErrorIndicator[conv, ops] gives a list of local error estimates computed from the object conv obtained by FrequencyConvergence. Option is ErrorFunction (by default Log10RelativeError).";
TransientSolution::usage = "TransientSolution[time, sys, ops] gives SimulationResult for a given list time. Options: Verbose, InputFunction, InitialState, NDSolve options. It uses different approaches for different types of DynamicSystem. Use ToFirstOrderSystem or ToExplicitSystem, if TransienSolution fails.";

AnsysTransientSolution::usage = "AnsysTransientSolution[time, sys, ops] gives SimulationResult for a given list time by using integrators like in ANSYS: backward Euler for a first order system and Newmark for a second order
system. Options: InputFunction, InitialState, TINTP.";
TINTP::usage = "An option for AnsysTransientSolution: TINTP->{gamma, alpha, delta, theta}";
InputFunction::usage = "An option for StationorySolution, TransientSolution, and HarmonicSolution.";
Verbose::usage = "An option for TransientSolution.";
InitialState::usage = "An option for TransientSolution.";

ModalSolution::usage = "ModalSolution[sys, n] gives solution for modal analysis for the first n lowest frequencies for the second order sys. n by default is equal to the system dimension.";

Begin["`Private`"];

(*SystemFirstOrder*)

Format[sys_DynamicSystem] := If[Length[sys] == 7 && MatrixQ[MatrixK[sys]], 
	"DynamicSystem[{"<>ToString[Length[MatrixK[sys]]]<>","<>ToString[Dimensions[MatrixB[sys]][[2]]]<>","<>ToString[Length[MatrixC[sys]]]<>"}, ...]",
	"- DynamicSystem - "
]

(*
SetAttributes[{
	DynamicSystem,
	SimulationResult,
	FunctionX,
	FunctionY,
	CommonTitle
}, Protected]
*)

checkConsistence[sys_DynamicSystem] := Module[{dim},
	dim = Dimensions[MatrixK[sys]];	
	Which[
		MatrixMQ[sys] && Dimensions[originalMatrixM[sys]] != dim,
			Throw["DynamicSystem: MatrixM and MatrixK are not compatible"],
		MatrixEQ[sys] && Dimensions[originalMatrixE[sys]] != dim,
			Throw["DynamicSystem: MatrixE and MatrixK are not compatible"],
		dim[[1]] != Length[MatrixB[sys]],
			Throw["DynamicSystem: MatrixB and MatrixK are not compatible"],
		dim[[1]] != Dimensions[MatrixC[sys]][[2]],
			Throw["DynamicSystem: MatrixC and MatrixK are not compatible"],
		Length[MatrixC[sys]] != Length[OutputNames[sys]],
			Throw["DynamicSystem: MatrixC and OutputNames are not compatible"],
		MatrixDQ[sys] && (Length[originalMatrixD[sys]] != Length[MatrixC[sys]]),
			Throw["DynamicSystem: MatrixD and MatrixC are not compatible"],
		MatrixDQ[sys] && (Dimensions[originalMatrixD[sys]][[2]] != Dimensions[MatrixB[sys]][[2]]),
			Throw["DynamicSystem: MatrixD and MatrixB are not compatible"]
	];
	sys
]

defaultMatrixC[n_Integer] := SparseArray[{{i_, i_}->1.}, {n, n}]

defaultOutputNames[n_Integer] := Table[ToString[i], {i, n}]

MakeDynamicSystem[
	matM_ ,
	matE_ ,
	matK_ ,
	matB_ ,
	matD_ : Null] := (
	checkConsistence[DynamicSystem[matM, matE, matK, matB, defaultMatrixC[Length[matK]], defaultOutputNames[Length[matK]], matD]]
) /; ((MatrixQ[matM] || matM === Null) 
	&& (MatrixQ[matE] || matE === Null)
	&& (MatrixQ[matK])
	&& (MatrixQ[matB])
	&& (MatrixQ[matD] || matD === Null))

MakeDynamicSystem[
	matM_ ,
	matE_ ,
	matK_ ,
	matB_ ,
	matC_ ,
	namesC:{___String},
	matD_ : Null] := (
	checkConsistence[DynamicSystem[matM, matE, matK, matB, matC, namesC, matD]]
) /; ((MatrixQ[matM] || matM === Null) 
	&& (MatrixQ[matE] || matE === Null)
	&& (MatrixQ[matK])
	&& (MatrixQ[matB])
	&& (MatrixQ[matC])
	&& (MatrixQ[matD] || matD === Null))

originalMatrixM[sys_DynamicSystem] := sys[[1]]

originalMatrixE[sys_DynamicSystem] := sys[[2]]

originalMatrixD[sys_DynamicSystem] := sys[[7]]

MatrixMQ[sys_DynamicSystem] := sys[[1]] =!= Null

MatrixEQ[sys_DynamicSystem] := sys[[2]] =!= Null

MatrixDQ[sys_DynamicSystem] := sys[[7]] =!= Null

FirstOrderSystemQ[sys_DynamicSystem] := sys[[1]] === Null

SecondOrderSystemQ[sys_DynamicSystem] := sys[[1]] =!= Null

ExplicitSystemQ[sys_DynamicSystem] := (sys[[1]] === Null && sys[[2]] === Null)

MatrixM[sys_DynamicSystem] := If[MatrixMQ[sys],
	sys[[1]],
	SparseArray[{},Dimensions[MatrixK[sys]]]
]

MatrixE[sys_DynamicSystem] := If[MatrixEQ[sys],
	sys[[2]],
	If[FirstOrderSystemQ[sys],
		SparseArray[{{i_, i_}->1.}, Dimensions[MatrixK[sys]]],
		SparseArray[{}, Dimensions[MatrixK[sys]]]
	]
]

MatrixK[sys_DynamicSystem] := sys[[3]]

MatrixA[sys_DynamicSystem] := -sys[[3]]

MatrixB[sys_DynamicSystem] := sys[[4]]

MatrixC[sys_DynamicSystem] := sys[[5]]

OutputNames[sys_DynamicSystem] := sys[[6]]

MatrixD[sys_DynamicSystem] := If[MatrixDQ[sys],
	sys[[7]],
	SparseArray[{},{Length[MatrixC[sys]], Dimensions[MatrixB[sys]][[2]]}]
]

MakeFirstOrderSystem[
	matA_,
	matB_,
	matC_,
	namesC:{___String},
	matD_ : Null] := 
MakeDynamicSystem[Null, Null, -matA, matB, matC, namesC, matD] /;
	((MatrixQ[matA])
	&& (MatrixQ[matB])
	&& (MatrixQ[matC])
	&& (MatrixQ[matD] || matD === Null))

MakeFirstOrderSystem[
	matE_, 
	matA_,
	matB_,
	matD_ : Null] := 
MakeDynamicSystem[Null, matE, -matA, matB, matD] /;
	((MatrixQ[matE] || matE === Null)
	&& (MatrixQ[matA])
	&& (MatrixQ[matB])
	&& (MatrixQ[matD] || matD === Null))

MakeFirstOrderSystem[
	matE_, 
	matA_,
	matB_,
	matC_,
	namesC:{___String},
	matD_ : Null] := 
MakeDynamicSystem[Null, matE, -matA, matB, matC, namesC, matD] /;
	((MatrixQ[matE] || matE === Null)
	&& (MatrixQ[matA])
	&& (MatrixQ[matB])
	&& (MatrixQ[matC])
	&& (MatrixQ[matD] || matD === Null))

MakeSecondOrderSystem[
	matM_,
	matK_,
	matB_,
	matC_,
	namesC:{___String},
	matD_ : Null] := 
MakeDynamicSystem[matM, Null, matK, matB, matC, namesC, matD] /;
	((MatrixQ[matM])
	&& (MatrixQ[matK])
	&& (MatrixQ[matB])
	&& (MatrixQ[matC])
	&& (MatrixQ[matD] || matD === Null))

MakeSecondOrderSystem[
	matM_,
	matE_,
	matK_,
	matB_,
	matD_ : Null] := 
MakeDynamicSystem[matM, matE, matK, matB, matD] /;
	((MatrixQ[matM])
	&& (MatrixQ[matE] || matE === Null)
	&& (MatrixQ[matK])
	&& (MatrixQ[matB])
	&& (MatrixQ[matD] || matD === Null))

MakeSecondOrderSystem[
	matM_,
	matE_,
	matK_,
	matB_,
	matC_,
	namesC:{___String},
	matD_ : Null] := 
MakeDynamicSystem[matM, matE, matK, matB, matC, namesC, matD] /;
	((MatrixQ[matM])
	&& (MatrixQ[matE] || matE === Null)
	&& (MatrixQ[matK])
	&& (MatrixQ[matB])
	&& (MatrixQ[matC])
	&& (MatrixQ[matD] || matD === Null))

DeleteOutputs[sys_DynamicSystem, name_String] := DeleteOutputs[sys, {name}] 

DeleteOutputs[sys_DynamicSystem, names:{___String}] := Module[{del,dd},
	del = Flatten[Map[Function[{n1},
		dd = {};
		MapIndexed[
			If[#1===n1,
				AppendTo[dd, #2]
			]&
		,OutputNames[sys]];
		dd		
	],names], 1];
	MakeDynamicSystem[
		originalMatrixM[sys],
		originalMatrixE[sys],
		MatrixK[sys], 
		MatrixB[sys], 
		Delete[MatrixC[sys], del], 
	  Delete[OutputNames[sys], del], 
		originalMatrixD[sys]
	]
]

AddOutputs[sys_DynamicSystem, matCadd_, name_String] := AddOutputs[sys, matCadd, {name}]

AddOutputs[sys_DynamicSystem, matCadd_, names:{___String}] := MakeDynamicSystem[
	originalMatrixM[sys],
	originalMatrixE[sys],
	MatrixK[sys], 
	MatrixB[sys], 
	Join[MatrixC[sys], matCadd], 
	Join[OutputNames[sys], names], 
	originalMatrixD[sys]
] /; MatrixQ[matCadd] && Dimensions[matCadd][[2]] == Length[MatrixK[sys]]

(*some functions for systemSecondOrder*)

DeleteDamping[sys_DynamicSystem] := MakeSecondOrderSystem[
	MatrixM[sys],
	MatrixK[sys],
	MatrixB[sys],
	MatrixC[sys],
	OutputNames[sys],
	originalMatrixD[sys]
]/; SecondOrderSystemQ[sys]

AddDamping[sys_DynamicSystem, matE_] :=
MakeSecondOrderSystem[
	MatrixM[sys],
	MatrixE[sys] + matE,
	MatrixK[sys],
	MatrixB[sys],
	MatrixC[sys],
	OutputNames[sys],
	originalMatrixD[sys]
]/; (SecondOrderSystemQ[sys] && MatrixQ[matE])

AddDamping[sys_DynamicSystem, alpha_, beta_ ] :=
MakeSecondOrderSystem[
	MatrixM[sys],
	MatrixE[sys] + N[alpha]*MatrixM[sys] + N[beta]*MatrixK[sys],
	MatrixK[sys],
	MatrixB[sys],
	MatrixC[sys],
	OutputNames[sys],
	originalMatrixD[sys]
] /; (SecondOrderSystemQ[sys] && NumericQ[alpha] && NumericQ[beta])

AddDampingRatio[sys_DynamicSystem, dmprat_] :=
MakeSecondOrderSystem[
	MatrixM[sys],
	MatrixE[sys],
	MatrixK[sys] + I*2*dmprat*MatrixK[sys],
	MatrixB[sys],
	MatrixC[sys],
	OutputNames[sys],
	originalMatrixD[sys]
]/; (SecondOrderSystemQ[sys] && NumericQ[dmprat])

AddDampingRatio[sys_DynamicSystem, imatK_] :=
MakeSecondOrderSystem[
	MatrixM[sys],
	MatrixE[sys],
	MatrixK[sys] + I imatK,
	MatrixB[sys],
	MatrixC[sys],
	OutputNames[sys],
	originalMatrixD[sys]
]/; (SecondOrderSystemQ[sys] && MatrixQ[matE])


(*functions for both systems*)

ToSingleInputSystem[sys_DynamicSystem] := (
	MakeDynamicSystem[
		originalMatrixM[sys],
		originalMatrixE[sys],
		MatrixK[sys], 
		Take[MatrixB[sys], Length[MatrixB[sys]], 1], 
		MatrixC[sys], 
	  OutputNames[sys], 
		originalMatrixD[sys]
  ]
)

ToFirstOrderSystem[sys_DynamicSystem] := Module[
	{zero = SparseArray[{}, Dimensions[MatrixK[sys]]], 
	one = SparseArray[{{i_, i_}->1.}, Dimensions[MatrixK[sys]]]},
	MakeFirstOrderSystem[
		BlockMatrix[{{MatrixE[sys], MatrixM[sys]},{-one, zero}}],
		-BlockMatrix[{{MatrixK[sys], zero},{zero, one}}],
		BlockMatrix[{{MatrixB[sys]},{SparseArray[{}, {Length[MatrixK[sys]],Dimensions[MatrixB[sys]][[2]]}] }}],
		BlockMatrix[{{MatrixC[sys], SparseArray[{}, {Length[MatrixC[sys]],Length[MatrixK[sys]]}] }}],
		OutputNames[sys],
		originalMatrixD[sys]
	]
] /; SecondOrderSystemQ[sys]

ToFirstOrderSystem[sys_DynamicSystem] := sys /; FirstOrderSystemQ[sys]

ToExplicitSystem[sys_DynamicSystem] := Module[{inv},
	If[MatrixEQ[sys],
			inv = Inverse[MatrixE[sys]];
			MakeFirstOrderSystem[
			inv.MatrixA[sys],
			inv.MatrixB[sys],
			MatrixC[sys],
			OutputNames[sys],
			originalMatrixD[sys]
		]
	,
		sys
	]
] /; FirstOrderSystemQ[sys]

ToExplicitSystem[sys_DynamicSystem] := ToExplicitSystem[ToFirstOrderSystem[sys]] /; SecondOrderSystemQ[sys]

TakeSystem[sys_DynamicSystem, dim_Integer] := (
	If[dim > Length[MatrixK[sys]],
		Throw["TakeSystem: dimension is smaller then required"]
	];
	DynamicSystem[
		If[MatrixMQ[sys], Take[MatrixM[sys], dim, dim], Null],
		If[MatrixEQ[sys], Take[MatrixE[sys], dim, dim], Null],
		Take[MatrixK[sys], dim, dim],
		Take[MatrixB[sys], dim],
		Take[MatrixC[sys], Length[MatrixC[sys]], dim],
		OutputNames[sys],
		If[MatrixDQ[sys], Take[MatrixD[sys], dim, dim], Null]
	]
)

fileQ[name_String] := Module[{in},
	Off[OpenRead::"noopen"];
	in = OpenRead[name];
	On[OpenRead::"noopen"];
	If [in === $Failed,
		False,
		Close[in]; True
	]
]

readMatrix1[name_String] := If[fileQ[name],
	Import[name, "MTX"],
	Throw["ReadSystem: No file "<>name]
]

readMatrix2[name_String] := If[fileQ[name],
	Import[name, "MTX"],
	Null
]

readNames[name_String] := Module[{in, Cnames},
	If[!fileQ[name], Throw["ReadSystem: No file "<>name]]; 
	in = OpenRead[name, DOSTextFormat->True];
	Cnames = ReadList[in, Word, WordSeparators->{" ", "\t", "\n", FromCharacterCode[10], FromCharacterCode[13]}];
	Close[in];
	Cnames 
]

readFirstOrder[baseName_String] := MakeFirstOrderSystem[
	readMatrix2[baseName<>".E"],
	readMatrix1[baseName<>".A"],
	readMatrix1[baseName<>".B"],
	readMatrix1[baseName<>".C"],
	readNames[baseName<>".C.names"],
	readMatrix2[baseName<>".D"]
]

readSecondOrder[baseName_String] := MakeSecondOrderSystem[
	readMatrix1[baseName<>".M"],
	readMatrix2[baseName<>".E"],
	readMatrix1[baseName<>".K"],
	readMatrix1[baseName<>".B"],
	readMatrix1[baseName<>".C"],
	readNames[baseName<>".C.names"],
	readMatrix2[baseName<>".D"]
]

ReadSystem[baseName_String] := (
	Which[
		fileQ[baseName<>".A"],
			readFirstOrder[baseName],
		fileQ[baseName<>".K"],
			readSecondOrder[baseName],
		True,
			Throw["ReadSystem: no files with A or K matrix"]
	]
)

ReadSystem[name_String, head_String] :=Module[
	{in, dimMax, Tini, redC, redK, ninp, funinp, bredin, nout, matout, names},
	in = OpenRead[name, DOSTextFormat->True];
	If[in == $Failed, Throw["ReadSystem: no such file"]];
	dimMax = Round[Read[in, Number]];
	Tini = N[Read[in, Number]];
	redC = N[Map[ReadList[in, Number, dimMax] &, Range[dimMax]]];
	redK = N[Map[ReadList[in, Number, dimMax] &, Range[dimMax]]];
	ninp = Round[Read[in, Number]];
	bredin = N[Map[ReadList[in, Number, dimMax] &, Range[ninp]]];
	nout = Round[Read[in, Number]];
	matout = {};
	names = {};
	Scan[
		(AppendTo[names, Read[in, Word, WordSeparators->{" ", "\t", "\n", FromCharacterCode[10], FromCharacterCode[13]}]]; 
		AppendTo[matout, N[ReadList[in, Number, dimMax]]])&, 
		Range[nout]
	];
	Close[in];
	Switch[
		head,
		"FirstOrderSystem",
			MakeFirstOrderSystem[redC, -redK, Transpose[bredin], matout, names],
		"SecondOrderSystem",
			MakeSecondOrderSystem[redC, redK, Transpose[bredin], matout, names],
		_,
			Throw["ReadSystem: wrong argument"]
	]
]

writeNames[name_String, list_List] := Module[{in},
	in = OpenWrite[name, FormatType -> OutputForm];
	Scan[
		Write[in, #]&
	,list];
	Close[in]
]

WriteSystem[sys_DynamicSystem, baseName_String] := (
	If[MatrixMQ[sys], Export[baseName<>".M", MatrixM[sys], "MTX"]];
	If[MatrixEQ[sys], Export[baseName<>".E", MatrixE[sys], "MTX"]];
	If[FirstOrderSystemQ[sys],
		Export[baseName<>".A", MatrixA[sys], "MTX"],
		Export[baseName<>".K", MatrixK[sys], "MTX"]
	];
	Export[baseName<>".B", MatrixB[sys], "MTX"];
	Export[baseName<>".C", MatrixC[sys], "MTX"];
	writeNames[baseName<>".C.names", OutputNames[sys]];
	If[MatrixDQ[sys], Export[baseName<>".D", MatrixD[sys], "MTX"]]
)

Format[res_SimulationResult] := "- SimulationResult - "

checkConsistence[res_SimulationResult] := Module[{dim},
	dim = Dimensions[YSeries[res]];	
	Which[
		Length[XSeries[res]] != dim[[2]],
			Throw["DynamicSystem: XSeries and YSeries are not compatible"],
		Length[YNames[res]] != dim[[1]],
			Throw["DynamicSystem: YNames and YSeries are not compatible"]
	];
	res
]

MakeSimulationResult[
	x_,
	xName_String,
	series_,
	names_List] := checkConsistence[SimulationResult[x, xName, series, names]] /;
	(VectorQ[x] && MatrixQ[series])

XSeries[res_SimulationResult] := res[[1]]

XName[res_SimulationResult] := res[[2]]

YSeries[res_SimulationResult] := res[[3]]

YNames[res_SimulationResult] := res[[4]]

Options[TransformResult] = {TransformFunction->Automatic};


TransformResult[res:{___SimulationResult}, names_List, ops___] := Map[
	TransformResult[#, names, ops]&
, res]

TransformResult[res_SimulationResult, names_List, ops___] := Module[{ser, nam, pos, flist, op},
	op = TransformFunction /. {ops} /. Options[TransformResult];
	If[VectorQ[op] && Dimensions[op]==={2}, op = {op}];
	flist = Which[
		op === Automatic,
			Table[#1&, {Length[names]}],
		MatrixQ[op] && Dimensions[op][[2]] == 2,
			Map[
			If[(pos = Position[Transpose[op][[1]], #, 1, 1]) == {},
				#1&
			,
				op[[pos[[1, 1]], 2]]
			]&
			,names],
		True,
			Table[op, {Length[names]}]
	];	
	ser = {};
	nam = {};
	MapThread[Function[{ni, fi},
		If[(pos = Position[YNames[res], ni, 1, 1]) == {},
			Print["Series "<>ni<>" is not present: ignored"]
		,
			AppendTo[ser, MapThread[Function[{y, x}, fi[y, x]], {Flatten[Extract[YSeries[res], pos]], XSeries[res]}]];
			AppendTo[nam, ni]
		]
	], {names, flist}];

	If[ser == {}, Throw["FindResults: empty SimulationResult"]];
	MakeSimulationResult[
		XSeries[res],
		XName[res],
		ser,
		nam
	]
]

compatibleQ[res1_SimulationResult, res2_SimulationResult] := 
	XName[res1] === XName[res2] && YNames[res1] === YNames[res2]

CompatibleResultQ[res1_SimulationResult] := True 

CompatibleResultQ[res1_SimulationResult, res2__] := 
Apply[And, Map[compatibleQ[res1, #]&, {res2}]] /; VectorQ[{res2}, MatchQ[#, _SimulationResult]&]

CompatibleResultQ[list_List] := Apply[CompatibleResultQ, list]

Options[Difference] = {ErrorFunction->Subtract};

Difference[res1_SimulationResult, res2_SimulationResult, ops___] := 
	Difference[res1, {res2}, ops]

diff1[res1_, resi_, fun_] := MapThread[
	fun[#1, #2]&
, {YSeries[res1], YSeries[resi]}]

diff2[res1_, resi_, fun_, ops___] := MapThread[
	fun[#1, Interpolation[Transpose[{XSeries[resi], #2}], FilterOptions[Interpolation, {ops}], InterpolationOrder->1][XSeries[res1]]]&
, {YSeries[res1], YSeries[resi]}]

Difference[res1_SimulationResult, resList_List, ops___] := Module[{
	fun = ErrorFunction /. {ops} /. Options[Difference]},
	Map[
		MakeSimulationResult[
			XSeries[res1],
			XName[res1],
			If[XSeries[res1] == XSeries[#],
				diff1[res1, #, fun],
				diff2[res1, #, fun, ops]
			],
			YNames[res1]
		]&
	, resList]
] /; CompatibleResultQ[Prepend[res1, resList]]

TakeResult[res_SimulationResult, num_Integer] := 	MakeSimulationResult[
		Take[XSeries[res], num],
		XName[res],
		Take[YSeries[res], Length[YSeries[res]], num],
		YNames[res]
	]

ReadResult[file_String] := Module[{mat, names},
	mat = Chop[Import[file, "MTX"], $MinNumber];
	names = readNames[file<>".names"];
	MakeSimulationResult[
		mat[[1]],
		names[[1]],
		Delete[mat, 1],
		Delete[names, 1]		
	]
]

WriteResult[res_SimulationResult, file_String] := ( 
  	Export[file, Prepend[YSeries[res], XSeries[res]], "MTX"];
	  writeNames[file<>".names", Prepend[YNames[res], XName[res]]]
)

WriteResult[res_SimulationResult, file_String, "Table"] := (
		Export[file, Transpose[Prepend[YSeries[res], XSeries[res]]], "Table"]
)

(*to read old format, the treatment of complex numbers is for Sam*)

ReadResult[name_String, nameX_String] := Module[{in, ntime, time, series, names, node, ser},
	Off[Read::"readn"];
	in = OpenRead[name];
	ntime = Round[Read[in, Number]];
	time = ReadList[in, Number, ntime];
	names = {};
	series = {};
	While[
    (node = Read[in, Word, WordSeparators->{" ", "\t", "\n", FromCharacterCode[10], FromCharacterCode[13]}]) =!= EndOfFile,
    AppendTo[names, node];
		ser = ReadList[in, Number, ntime*2];
		Which[
			Length[ser] === ntime*2,
				ser = Map[(#[[1]]*Cos[#[[2]]*Degree] +
I*#[[1]]*Sin[#[[2]]*Degree])&, Partition[ser, 2]],
			Length[ser] === ntime,
				ser,
			True,
				ser = Delete[ser, -1]
		];
    AppendTo[series, ser]
  ];
	Close[in];
	On[Read::"readn"];
	MakeSimulationResult[
		time,
		nameX,
		series,
		names
	]
]

(*for Sam*)

formatNumber[x_] := ToString[NumberForm[CForm[x],25,NumberPadding->{" ",""},SignPadding->False,ExponentFunction -> Identity]]

WriteResultOld[res_SimulationResult, file_String] := Module[{of},
	of = OpenWrite[file, FormatType -> OutputForm];
	Write[of, ToString[Length[XSeries[res]]]];
	Scan[
		Write[of, formatNumber[#]]&
	,XSeries[res]];
	MapIndexed[(
		Write[of, #1];
		Scan[
			If [Head[#] === Complex,
				Write[of, formatNumber[Abs[#]]<>" "<>formatNumber[Arg[#]/Degree]]				
			,
				Write[of, formatNumber[#]]
			]&
		, YSeries[res][[#2[[1]]]]]
  )&, YNames[res]];
	Close[of]		
]

(*tuned for harmonic simulation*)
WriteResultAsDirichlet[res_SimulationResult, file_String] := 
Module[{of,val},
	of = OpenWrite[file, FormatType -> OutputForm];
	Write[of, "resume"];
	Write[of, "/solu"];
	Write[of, "allsel"];
	Write[of, "antype,3"];
	MapThread[Function[{xi, yi},
		Write[of, "harfrq,"<>ToString[CForm[xi]]];
		MapThread[(
			If[Head[#2] === Complex,
				val = ToString[CForm[Re[#2]]]<>","<>ToString[CForm[Im[#2]]]
			,
				val = ToString[CForm[#2]]
			];
			Write[of, "D,"<>StringReplace[#1, {"Node_" -> "", "_" -> ","}]<>","<>val];
  	)&, {YNames[res], yi}];
		Write[of, "solve"]
	], {XSeries[res], Transpose[YSeries[res]]}];
	Write[of, "fini"];
	Close[of]	
]

(*does not handle complex-valued numbers correctly*)
WriteResultAsDirichlet[res_SimulationResult, i_Integer, file_String] := 
Module[{of},
	If[i > Length[XSeries[res]] || i < 1, Throw["The index i is not valid"]];
	of = OpenWrite[file, FormatType -> OutputForm];
	MapIndexed[(
		Write[of, "D,"<>StringReplace[#1, {"Node_" -> "", "_" -> ","}]<>","<>ToString[CForm[YSeries[res][[#2[[1]], i]]]]];
  )&, YNames[res]];
	Close[of]	
]

defaultPlotOptions = {PlotRange->All, ImageSize->300, Frame->True, Axes->False};

Log10 := Log[10, #] &

Options[PlotResults] = {CommonTitle->"", FunctionX->Identity, FunctionY->Identity};

PlotResult[res_SimulationResult, ops___] := PlotResult[{res}, ops]

PlotResult[listres_List, ops___] := Module[{
		options = Join[{ops}, defaultPlotOptions],
		title = CommonTitle /. {ops} /. Options[PlotResults],	
		funx = FunctionX /. {ops} /. Options[PlotResults],
		funy = FunctionY /. {ops} /. Options[PlotResults]
	},
	MapIndexed[Function[{name, i},
		Apply[MultipleListPlot, Join[Map[Function[{res},
			Transpose[{funx[XSeries[res]], funy[YSeries[res][[i[[1]]]]]}]]
			, listres],
 			{FilterOptions[MultipleListPlot, options]},
			{PlotLabel->name<>": "<>title,
			PlotJoined->True,
			SymbolShape->None}]
		]
	],YNames[listres[[1]]]]
] /; CompatibleResultQ[listres]

makeInputs[inpFun_, n_] := Module[{inp},
	inp = Which[
	inpFun === Automatic,
		Table[1., {n}],
	Head[inpFun] === List && Length[inpFun] === n,
		inpFun,
	n === 1 && Head[inpFun] =!= List,
		{inpFun},
	True,
		Throw["Input functions are not compatible with the system"] 
	];
	If[VectorQ[inp, NumericQ],
		inp
	,
		Map[Function[{u},
			If[NumericQ[u[0, Table[0, {n}]]], 
				u
			, 
				If[NumericQ[(u&)[0, Table[0, {n}]]], 
					u&
				,
					Throw[{"Not numeric values of input function ", u}]
				]
			]
		], inp]
	]
]

Options[StationarySolution] = {InputFunction->Automatic};

StationarySolution[sys_DynamicSystem, ops___] := Module[{inpFun},
	inpFun = InputFunction /. {ops} /. Options[StationarySolution];
	inpFun = makeInputs[inpFun, Last[Dimensions[MatrixB[sys]]]];
	If[!VectorQ[inpFun, NumericQ], Throw["Input functions must be numeric vales"]]; 
	MatrixC[sys].LinearSolve[MatrixK[sys], MatrixB[sys].inpFun, FilterOptions[LinearSolve, ops]]]

harmonicLoop[list_, fun_, rhs_, ops___] := Map[
 LinearSolve[fun[#], rhs[#], FilterOptions[LinearSolve, ops]]&, 
2.*N[Pi]*list]

Options[HarmonicSolution] = {InputFunction->Automatic};

HarmonicSolution[freq_List, sys_DynamicSystem, ops___] := Module[{inpFun, tmp, f},
	inpFun = InputFunction /. {ops} /. Options[HarmonicSolution];
	inpFun = makeInputs[inpFun, Last[Dimensions[MatrixB[sys]]]];
	If[VectorQ[inpFun, NumericQ], 
		f[w_] := MatrixB[sys].inpFun;
		,
		f[w_] := MatrixB[sys].Map[Function[el, el[w]], inpFun];
	];
MakeSimulationResult[
	freq,
	"Frequency",
	MatrixC[sys].Transpose[Which[
	!MatrixMQ[sys],
		harmonicLoop[freq, (I*#*MatrixE[sys] + MatrixK[sys])&, 
		f, ops], 
	!MatrixEQ[sys],
		harmonicLoop[freq, (-#^2*MatrixM[sys] + MatrixK[sys])&, 
		f, ops], 
	True,
		harmonicLoop[freq, (-#^2*MatrixM[sys] + I*#*MatrixE[sys] + MatrixK[sys])&, 
		f, ops]
	]], 
	OutputNames[sys]
]] /; VectorQ[freq, NumericQ]

FrequencyConvergence[freq_, sys_DynamicSystem, 
start_Integer: 1, 
finish_: Infinity,
step_Integer: 1] := Module[{fin},
	fin = finish;
  If[fin > Length[MatrixK[sys]], fin = Length[MatrixK[sys]]];
  Map[
  (
    {#, Transpose[YSeries[HarmonicSolution[{freq}, TakeSystem[sys, #]]]][[1]]}
   ) &,
   Range[start, fin, step]
   ]
]/; NumericQ[freq]

Log10RelativeError[x_, y_] := Log10[Norm[x - y]/Norm[x]]

Options[LocalErrorIndicator] = {ErrorFunction->Log10RelativeError};

LocalError[list_List, exact_List, ops___] := Module[{err},
	err = ErrorFunction /. {ops} /. Options[LocalErrorIndicator];
	Map[{#[[1]], err[exact, #[[2]]]}&, list]
]

LocalErrorIndicator[list_List, ops___] := Module[{l, prev, p, err},
	err = ErrorFunction /. {ops} /. Options[LocalErrorIndicator];
	l1 = Delete[list, -1];
	l2 = Delete[list, 1];
	MapThread[({#1[[1]], err[#2[[2]], #1[[2]]]})&, {l1, l2}]
]

Options[TransientSolution] = {Verbose->False, InputFunction->Automatic, InitialState->Automatic, TINTP->{0.005,Null,Null,1}};

AnsysTransientSolution[time_List, sys_DynamicSystem, ops___] := Module[
  {init, inpFun, f, theta, dt, dtprev, un, vn, res, tprev, un1, tintp, factor},
	tintp = TINTP /. {ops} /. Options[TransientSolution];
  theta = tintp[[4]];
	init = InitialState /. {ops} /. Options[TransientSolution];
	If[init === Automatic,
		init = Table[0., {Length[MatrixK[sys]]}]];
	inpFun = InputFunction /. {ops} /. Options[TransientSolution];
	inpFun = makeInputs[inpFun, Last[Dimensions[MatrixB[sys]]]];
	If[VectorQ[inpFun, NumericQ], 
		f[tau_, x_] := MatrixB[sys].inpFun,
		f[tau_, x_] := MatrixB[sys].Map[Function[el, el[tau, x]], inpFun]
	];
	tprev = 0;
	dtprev = 0;
	un = init;
	vn = Table[0., {Length[MatrixK[sys]]}];
	res = Map[
		(
		dt = # - tprev;
		tprev = #;
		If[dt != 0,
			If[dt != dtprev,
				factor = LinearSolve[1./theta/dt*MatrixE[sys] + MatrixK[sys], FilterOptions[LinearSolve, ops]];
			];
 	  	un1 = factor[f[#, un] + MatrixE[sys].(1./theta/dt*un + (1 - theta)/theta*vn)];
 	  	vn = 1./theta/dt((un1-un) - (1 - theta)*dt*vn);
			dtprev = dt;
			un = un1
		,
			un
		]
		)&, time];
	MakeSimulationResult[
		time,
		"Time",
		MatrixC[sys].Transpose[res], 
		OutputNames[sys]
	]
]/; FirstOrderSystemQ[sys]

AnsysTransientSolution[time_List, sys_DynamicSystem, ops___] := Module[
  {init, inpFun, f, gamma, alpha, beta, dt, dtprev, un, vn, an, res, tprev, a0, a1, a2, a3, a4, a5, a6, a7, vn1, un1, an1, tintp, factor},
	tintp = TINTP /. {ops} /. Options[TransientSolution];
	If[tintp[[1]] =!= Null,
	  gamma = tintp[[1]];
		alpha = 0.25*(1 + gamma)^2;
		delta = 0.5 + gamma
	,
		alpha = tintp[[2]];
		beta = tintp[[3]]
	];
	init = InitialState /. {ops} /. Options[TransientSolution];
	If[init === Automatic,
		init = Table[0., {Length[MatrixK[sys]]}]];
	inpFun = InputFunction /. {ops} /. Options[TransientSolution];
	inpFun = makeInputs[inpFun, Last[Dimensions[MatrixB[sys]]]];
	If[VectorQ[inpFun, NumericQ], 
		f[tau_, x_] := MatrixB[sys].inpFun,
		f[tau_, x_] := MatrixB[sys].Map[Function[el, el[tau, x]], inpFun]
	];

	tprev = 0;
	dtprev = 0;
	un = init;
  vn = Table[0., {Length[MatrixK[sys]]}];
	an = Table[0., {Length[MatrixK[sys]]}];
	res = Map[
		(
		dt = # - tprev;
		tprev = #;
		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 != dtprev,
	 			factor = LinearSolve[a0*MatrixM[sys] + a1*MatrixE[sys] + MatrixK[sys], FilterOptions[LinearSolve, ops]]
			];
	    un1 = factor[f[#, un] + MatrixM[sys].(a0*un + a2*vn + a3*an) + MatrixE[sys].(a1*un + a4*vn + a5*an)];
			an1 = a0*(un1 - un) - a2*vn - a3*an;
	    vn1 = vn + a6*an + a7*an1;
			vn = vn1;
			an = an1;
			dtprev = dt;
			un = un1
		,
			un
		]
		)&, time];
	MakeSimulationResult[
		time,
		"Time",
		MatrixC[sys].Transpose[res], 
		OutputNames[sys]
	]
]/; SecondOrderSystemQ[sys]

ModalSolution[sys_DynamicSystem, n_: Infinity] :=Module[{dim, values, vectors},
	dim = n;
  If[dim > Length[MatrixK[sys]], dim = Length[MatrixK[sys]]];
	{values, vectors}=Eigensystem[{MatrixK[sys],MatrixM[sys]}, -dim];
	MakeSimulationResult[
		Sqrt[Reverse[values]]/(2*Pi),
		"Modal",
		MatrixC[sys].Transpose[Reverse[vectors]], 
		OutputNames[sys]
	]
]/; SecondOrderSystemQ[sys] 

TransientSolution[time_List, sys_DynamicSystem, ops___] := Module[{sol, tim, ver, inpFun, f1, init, init2, initv},
	ver = Verbose /. {ops} /. Options[TransientSolution];
	init = InitialState /. {ops} /. Options[TransientSolution];
	inpFun = InputFunction /. {ops} /. Options[TransientSolution];
	inpFun = makeInputs[inpFun, Last[Dimensions[MatrixB[sys]]]];
	If[VectorQ[inpFun, NumericQ], 
		f1[tau_, x_?VectorQ] := -MatrixK[sys].x + MatrixB[sys].inpFun,
		f1[tau_, x_?VectorQ] := -MatrixK[sys].x + MatrixB[sys].Map[Function[el, el[tau, x]], inpFun]
	];
	If[init === Automatic,
		init = Table[0., {Length[MatrixK[sys]]}]];
	If[Length[init] != Length[MatrixK[sys]], Throw["Initial State is not compatible"]];
	tim = AbsoluteTiming[sol = Which[
	ExplicitSystemQ[sys],
		NDSolve[{T'[tau] == f1[tau, T[tau]], 
			T[0] == init}, T, {tau, 0, Last[time]}, FilterOptions[NDSolve, ops]],
	FirstOrderSystemQ[sys],
		init2 = LinearSolve[MatrixE[sys], f1[0, init]];
		NDSolve[{(MatrixE[sys].#)&[T'[tau]] == f1[tau, T[tau]], 
			T[0] == init, T'[0] == init2}, T, {tau, 0, Last[time]}, SolveDelayed->True, FilterOptions[NDSolve, ops]],
	True,
		initv = Table[0., {Length[MatrixK[sys]]}];
		init2 = LinearSolve[MatrixM[sys], f1[0, init] - MatrixE[sys].initv];
		NDSolve[{(MatrixM[sys].#)&[T''[tau]] + (MatrixE[sys].#)&[T'[tau]] == f1[tau, T[tau]], 
			T[0] == init, T'[0] == initv, T''[0] == init2}, T, {tau, 0, Last[time]}, SolveDelayed->True, FilterOptions[NDSolve, ops]]
	];];
If[ver, Print["NDSolve has made ", Length[sol[[1, 1, 2, 3, 1]]], " steps for ", tim[[1]]]];
MakeSimulationResult[
	time,
	"Time",
	MatrixC[sys].Transpose[T[time] /. First[sol]], 
	OutputNames[sys]
]
]	/; VectorQ[time, NumericQ]

End[];

EndPackage[];