#This example repeats the example from Mathematica. Look at 
#SecondOrderSystem.nb.pdf in the parent directory: there are more comments there. 

#Setting Up Directory and Loading Functions

import os
os.chdir('.')
print 'current directory is', os.getcwd()

#in order for import to work post4mor.py should be in:
#  the directory in which python has started, or
#  Python25/Lib/site-packages, or
#  on the path

import sys
sys.path.append('.')

from post4mor import *

filePath = '../ex/'

#we need the second argument True to force removing of empty lines 
#produced by MOR for ANSYS
sys = ReadSystem(filePath + 'gyro', True)

#Transient Simulation

transient = ReadResult(filePath + 'gyro.transient')

sysDamped = AddDamping(sys, 0., 1e-6)

#subsection Using NDSolve with TransientSolution Function is skipped
#as there is no NDSolve in Python

#Using Ansys Integrator with AnsysTransientSolution Function

res10 = AnsysTransientSolution(XSeries(transient), TakeSystem(sysDamped, 10))

res5 = AnsysTransientSolution(XSeries(transient), TakeSystem(sysDamped, 5))

PlotResult([transient, res5, res10],
		FileName = 'gyroTransient',
		Options = [{'color':(1,0,0)}, {'color':(0,1,0)}, {'color':(0,0,1)}])

#Harmonic Simulation

harmonic = ReadResult(filePath + 'gyro.harmonic')

sysDamped = AddDamping(sys, 0., 1e-7)

res30 = HarmonicSolution(XSeries(harmonic), TakeSystem(sysDamped, 30))

res20 = HarmonicSolution(XSeries(harmonic), TakeSystem(sysDamped, 20))

res10 = HarmonicSolution(XSeries(harmonic), TakeSystem(sysDamped, 10))

from math import log10

PlotResult([harmonic, res10, res20, res30], 
		FileName = 'gyroHarmonic', 
		Options = [{'color':(1,0,0)}, {'color':(0,1,0)}, {'color':(0,0,1)}, {'color':(0,1,1)}], 
		LogY = True,
		FunctionX = lambda x: log10(x),
		FunctionY = lambda y: abs(y))

#Local Error Indicator

er1 = LocalErrorIndicator(FrequencyConvergence(1e3, sysDamped))
er2 = LocalErrorIndicator(FrequencyConvergence(1e4, sysDamped))
er3 = LocalErrorIndicator(FrequencyConvergence(1e5, sysDamped))

PlotConvergence([er1, er2, er3],
		FileName = 'gyroError1',
		Options = [
			{'marker':'o', 'markerfacecolor':(1,0,0), 'linewidth':0}, 
			{'marker':'+', 'markerfacecolor':(0,1,0), 'linewidth':0}, 
			{'marker':'d', 'markerfacecolor':(0,0,1), 'linewidth':0}], 
		Legend = ('1e3 Hz', '1e4 Hz', '1e5 Hz'))

freq = XSeries(harmonic)[50]
exact = YSeries(harmonic).transpose()[50]

from scipy.linalg import norm

error = lambda x, y: norm(abs(x) - abs(y))/norm(x)

conv = FrequencyConvergence(freq, sysDamped)

erEst = LocalErrorIndicator(conv, ErrorFunction = error)

erTrue = LocalError(conv, exact, ErrorFunction = error)

PlotConvergence([erTrue, erEst],
		FileName = 'gyroError2',
		Options = [
			{'marker':'o', 'markerfacecolor':(1,0,0), 'linewidth':0}, 
			{'marker':'+', 'markerfacecolor':(0,1,0), 'linewidth':0}], 
		Legend = ('True', 'Indicator'))