#!/usr/bin/env python
# coding: utf-8

# ### Package 内の Model のシミュレーションを実行する
# [Modelicaのクラスの概要](https://www.amane.to/wp-content/uploads/2018/12/8ec4f21241c98ee8413280240090c942.pdf) の ClassExample2.mo の中のモデル MassB　のシミュレーションを実行します。
# 
# このモデルは、質点の運動状態を表すレコード State を定義して、ClassExample1.mo の　ModelA を書き直したものです。
# 
# Modelica のコードをMagicコマンド %%writefile で出力します。

# In[1]:


get_ipython().run_cell_magic('writefile', 'ClassExample2.mo', 'package ClassExample2 // package\n\n  constant Acceleration g = -9.8;\n  \n  // type\n  type Acceleration = Real(quantity = "Acceleration", unit = "m/s2");\n  type Mass = Real(quantity = "Mass", unit = "kg");\n  type Force = Real(quantity = "Force", unit = "N");\n  type Velocity = Real(quantity = "Velocity", unit = "m/s");\n  type Position = Real(quantity = "Length", unit = "m");\n  \n  // record\n  record State\n    Position x;\n    Velocity v;\n  end State;\n  \n  // model\n  model MassB\n    parameter Mass m = 1.0;\n    parameter Force f = m * g;\n    parameter State s0(x = 0, v = 5);\n    State s(x(start = s0.x), v(start = s0.v));\n  equation\n    s.v = der(s.x);\n    f = m * der(s.v);\n  end MassB;\n\nend ClassExample2;\n')


# クラス MassB をコンパイルして FMU を作成します。

# In[2]:


# Import the function for compilation of models and the load_fmu method
from pymodelica import compile_fmu
from pyfmi import load_fmu
# Compile model
fmu_name = compile_fmu("ClassExample2.MassB","ClassExample2.mo")
# Load model
model = load_fmu(fmu_name)


# シミュレーションを実行します。オプション ncp は Number of communication points で、ソルバーとFMUが通信する回数を表しています。これを 0 にするとソルバー内部の時間ステップになります。オプションの詳細は [JModelica User Guide - Version 2.10](https://jmodelica.org/downloads/UsersGuide-2.10.pdf) P.26-2９ を参照してください。

# In[3]:


opts = model.simulate_options()
opts['ncp']=500
res = model.simulate(final_time=1.,options=opts)


# シミュレーション結果をプロットします。

# In[4]:


get_ipython().run_line_magic('matplotlib', 'notebook')
# Import the plotting library
import matplotlib.pyplot as plt 
x = res['s.x']
v = res['s.v']
t = res['time']

plt.figure(1)
plt.subplot(2,1,1)
plt.plot(t, x)
plt.ylabel('Position (m)')
plt.subplot(2,1,2)
plt.plot(t,v)
plt.ylabel('velocity (m/s)')
plt.xlabel('Time (s)')
plt.show()


# In[ ]: