# ------------------------------------------------------------------------------------------------------------------ #
#
#   DYNARDO example
#
# Damped oscillator
#
# Undamped eigenfrequency: omega_0 = sqrt(k/m)
# Damped eigenfrequency:     omega = sqrt(1-D^2)*omega_0
# Response:                   x(t) = exp(-D*omega_0*t)*sqrt(2*Ekin/m)/omega*sin(omega*t)
#
# Minimize maximum amplitude after 5 seconds
# Maximum eigenfrequencies as constraint: omega <= 8 1/s
#
# 0.1 kg <= m <= 5.0 kg
# 10 N/m <= k <= 50 N/m
# D = 0.02
# Ekin = 10 Nm
#
# For background information see ../doc/oszillator.pdf
# ------------------------------------------------------------------------------------------------------------------ #
from math import *

# values to default if not existing yet
try:
    test = m
except:
    m = 1.00000000
try:
    test = k
except:
    k = 20.00000000
try:
    test = D
except:
    D = 0.02000000
try:
    test = Ekin
except:
    Ekin = 10.00000000

#scalar values
v0 = sqrt(2*Ekin/m)
omega0 = sqrt(k/m)
omega_damped = omega0*sqrt(1-D*D)

#approximate maximum amplitude from envelope
x_max_env = exp(-5*D*omega0)*v0/omega_damped
numsteps=100

#time series
def CalcDisplacement(_t) :
    sin_value = sin(_t*omega_damped)
    envelope =  exp(-_t*D*omega0)*v0/omega_damped
    return sin_value*envelope

def CalcIndicator(_t,_x) :
    if _t > 5:
        return abs(_x)
    else:
        return 0

time = range(0,numsteps)
time = map(lambda t: 10.*t/(numsteps-1), time)
x_values = map(lambda t: CalcDisplacement(t), time)

#extract maximum amplitude from time series
x_max=max(map(lambda t, x: CalcIndicator(t,x), time, x_values))

print("omega_damped: " + str(omega_damped))
print("x_max_env: " + str(x_max_env))
print("x_max: " + str(x_max))

try :
    from pyvariant import List2VariantVector
    time = List2VariantVector(time)
    x_values = List2VariantVector(x_values)
except :
    print "Did not export result to optiSLang"
