diff --git a/LICENSE.txt b/LICENSE.txt
index 951a1e0..bee66a3 100644
--- a/LICENSE.txt
+++ b/LICENSE.txt
@@ -1,3 +1,25 @@
+ License Exception for "alphanumeric.csl"
+
+The "alphanumeric.csl" file located in the "doc/" directory of this project
+is sourced from an external repository and is subject to different
+licensing terms. Specifically, this file is released under the CREATIVE
+COMMONS ATTRIBUTION-SHAREALIKE 3.0 UNPORTED license.
+
+ Conditions for Use:
+
+ * Attribution: Any software using the "alphanumeric.csl" file must
+ include a clear mention of the CSL project and a link to
+ [https://citationstyles.org/].
+ * Redistribution: When redistributing the "alphanumeric.csl" file,
+ the listings of authors and contributors in the style metadata must
+ be preserved as is.
+
+This exception applies solely to the "alphanumeric.csl" file and does
+not affect the licensing of the rest of the project, which remains under
+the CeCILL FREE SOFTWARE LICENSE AGREEMENT, Version 2.1 dated 2013-06-21,
+included below.
+
+
CeCILL FREE SOFTWARE LICENSE AGREEMENT
diff --git a/doc/alphanumeric.csl b/doc/alphanumeric.csl
new file mode 100644
index 0000000..b9b2581
--- /dev/null
+++ b/doc/alphanumeric.csl
@@ -0,0 +1,361 @@
+
+
diff --git a/doc/data/sincos_discrete.csv b/doc/data/sincos_discrete.csv
new file mode 100644
index 0000000..92fd817
--- /dev/null
+++ b/doc/data/sincos_discrete.csv
@@ -0,0 +1,1001 @@
+0.0000000000e+00,0.0000000000e+00,1.0000000000e+00
+1.0000000000e-01,9.9880025199e-02,9.9550209979e-01
+2.0000000000e-01,1.9886154962e-01,9.8104841125e-01
+3.0000000000e-01,2.9595423026e-01,9.5677345681e-01
+4.0000000000e-01,3.9018561464e-01,9.2291006930e-01
+5.0000000000e-01,4.8061087966e-01,8.7978716289e-01
+6.0000000000e-01,5.6632230388e-01,8.2782654125e-01
+7.0000000000e-01,6.4645837847e-01,7.6753877409e-01
+8.0000000000e-01,7.2021246529e-01,6.9951818215e-01
+9.0000000000e-01,7.8684091515e-01,6.2443697999e-01
+1.0000000000e+00,8.4567056453e-01,5.4303863433e-01
+1.1000000000e+00,8.9610553520e-01,4.5613050345e-01
+1.2000000000e+00,9.3763326811e-01,3.6457583053e-01
+1.3000000000e+00,9.6982973037e-01,2.6928517038e-01
+1.4000000000e+00,9.9236374263e-01,1.7120733464e-01
+1.5000000000e+00,1.0050003824e+00,7.1319945516e-02
+1.6000000000e+00,1.0076034290e+00,-2.9380308005e-02
+1.7000000000e+00,1.0001368234e+00,-1.2988761419e-01
+1.8000000000e+00,9.8266512958e-01,-2.2919708378e-01
+1.9000000000e+00,9.5535298938e-01,-3.2631479607e-01
+2.0000000000e+00,9.1846357692e-01,-4.2026774534e-01
+2.1000000000e+00,8.7235606641e-01,-5.1011358816e-01
+2.2000000000e+00,8.1748213783e-01,-5.9495009404e-01
+2.3000000000e+00,7.5438155437e-01,-6.7392420442e-01
+2.4000000000e+00,6.8367685490e-01,-7.4624060926e-01
+2.5000000000e+00,6.0606721377e-01,-8.1116975496e-01
+2.6000000000e+00,5.2232152836e-01,-8.6805520293e-01
+2.7000000000e+00,4.3327080270e-01,-9.1632026467e-01
+2.8000000000e+00,3.3979990274e-01,-9.5547384625e-01
+2.9000000000e+00,2.4283876485e-01,-9.8511544308e-01
+3.0000000000e+00,1.4335314504e-01,-1.0049392341e+00
+3.1000000000e+00,4.2335000875e-02,-1.0147372334e+00
+3.2000000000e+00,-5.9207398179e-02,-1.0144014676e+00
+3.3000000000e+00,-1.6025953335e-01,-1.0039251546e+00
+3.4000000000e+00,-2.5981077170e-01,-9.8340287317e-01
+3.5000000000e+00,-3.5686447253e-01,-9.5302971876e-01
+3.6000000000e+00,-4.5044796407e-01,-9.1309945368e-01
+3.7000000000e+00,-5.3962229052e-01,-8.6400166945e-01
+3.8000000000e+00,-6.2349163182e-01,-8.0621798819e-01
+3.9000000000e+00,-7.0121230165e-01,-7.4031734023e-01
+4.0000000000e+00,-7.7200123329e-01,-6.6695036435e-01
+4.1000000000e+00,-8.3514386798e-01,-5.8684298553e-01
+4.2000000000e+00,-8.9000136638e-01,-5.0078923377e-01
+4.3000000000e+00,-9.3601707034e-01,-4.0964337487e-01
+4.4000000000e+00,-9.7272214957e-01,-3.1431143127e-01
+4.5000000000e+00,-9.9974037608e-01,-2.1574217701e-01
+4.6000000000e+00,-1.0167919777e+00,-1.1491769627e-01
+4.7000000000e+00,-1.0236965313e+00,-1.2843599585e-02
+4.8000000000e+00,-1.0203748655e+00,8.9461004982e-02
+4.9000000000e+00,-1.0068499537e+00,1.9097368558e-01
+5.0000000000e+00,-9.8324678657e-01,2.9067890375e-01
+5.1000000000e+00,-9.4979122441e-01,3.8757817287e-01
+5.2000000000e+00,-9.0680784059e-01,4.8070005635e-01
+5.3000000000e+00,-8.5471677567e-01,5.6910990543e-01
+5.4000000000e+00,-7.9402963321e-01,6.5191923896e-01
+5.5000000000e+00,-7.2534445714e-01,7.2829467105e-01
+5.6000000000e+00,-6.4933984005e-01,7.9746629696e-01
+5.7000000000e+00,-5.6676822042e-01,8.5873545272e-01
+5.8000000000e+00,-4.7844843486e-01,9.1148177048e-01
+5.9000000000e+00,-3.8525759934e-01,9.5516945817e-01
+6.0000000000e+00,-2.8812239955e-01,9.8935273999e-01
+6.1000000000e+00,-1.8800987715e-01,1.0136804026e+00
+6.2000000000e+00,-8.5917803329e-02,1.0278994006e+00
+6.3000000000e+00,1.7135264410e-02,1.0318574840e+00
+6.4000000000e+00,1.2012014321e-01,1.0255048214e+00
+6.5000000000e+00,2.2200730219e-01,1.0088946001e+00
+6.6000000000e+00,3.2177715358e-01,9.8218259794e-01
+6.7000000000e+00,4.1843025469e-01,9.4562572842e-01
+6.8000000000e+00,5.1099731874e-01,8.9957957387e-01
+6.9000000000e+00,5.9854893430e-01,8.4449492965e-01
+7.0000000000e+00,6.8020489577e-01,7.8091339308e-01
+7.1000000000e+00,7.5514305141e-01,7.0946204044e-01
+7.2000000000e+00,8.2260757980e-01,6.3084724397e-01
+7.3000000000e+00,8.8191661162e-01,5.4584769022e-01
+7.4000000000e+00,9.3246911976e-01,4.5530666839e-01
+7.5000000000e+00,9.7375100822e-01,3.6012370525e-01
+7.6000000000e+00,1.0053403381e+00,2.6124562953e-01
+7.7000000000e+00,1.0269116377e+00,1.5965715445e-01
+7.8000000000e+00,1.0382392522e+00,5.6371072256e-02
+7.9000000000e+00,1.0391996998e+00,-4.7581841873e-02
+8.0000000000e+00,1.0297730076e+00,-1.5116311569e-01
+8.1000000000e+00,1.0100430156e+00,-2.5333695304e-01
+8.2000000000e+00,9.8019664167e-01,-3.5308059056e-01
+8.3000000000e+00,9.4052211671e-01,-4.4939453456e-01
+8.4000000000e+00,8.9140620464e-01,-5.4131257551e-01
+8.5000000000e+00,8.3333043481e-01,-6.2791147974e-01
+8.6000000000e+00,7.6686638325e-01,-7.0832026139e-01
+8.7000000000e+00,6.9267004923e-01,-7.8172894122e-01
+8.8000000000e+00,6.1147538212e-01,-8.4739670443e-01
+8.9000000000e+00,5.2408702268e-01,-9.0465937519e-01
+9.0000000000e+00,4.3137233036e-01,-9.5293613262e-01
+9.1000000000e+00,3.3425277572e-01,-9.9173540022e-01
+9.2000000000e+00,2.3369478333e-01,-1.0206598490e+00
+9.3000000000e+00,1.3070011608e-01,-1.0394104637e+00
+9.4000000000e+00,2.6295896688e-02,-1.0477896301e+00
+9.5000000000e+00,-7.8475634285e-02,-1.0457032117e+00
+9.6000000000e+00,-1.8256752185e-01,-1.0331615947e+00
+9.7000000000e+00,-2.8493855746e-01,-1.0102796882e+00
+9.8000000000e+00,-3.8456369298e-01,-9.7727588070e-01
+9.9000000000e+00,-4.8044430346e-01,-9.3446995997e-01
+1.0000000000e+01,-5.7161819607e-01,-8.8228001820e-01
+1.0100000000e+01,-6.5716926492e-01,-8.2121837090e-01
+1.0200000000e+01,-7.3623669472e-01,-7.5188652987e-01
+1.0300000000e+01,-8.0802362109e-01,-6.7496927967e-01
+1.0400000000e+01,-8.7180516014e-01,-5.9122791557e-01
+1.0500000000e+01,-9.2693572663e-01,-5.0149271004e-01
+1.0600000000e+01,-9.7285556675e-01,-4.0665468214e-01
+1.0700000000e+01,-1.0090964394e+00,-3.0765675144e-01
+1.0800000000e+01,-1.0352863884e+00,-2.0548436428e-01
+1.0900000000e+01,-1.0511535570e+00,-1.0115568556e-01
+1.1000000000e+01,-1.0565290056e+00,4.2885463849e-03
+1.1100000000e+01,-1.0513485035e+00,1.0979540064e-01
+1.1200000000e+01,-1.0356532754e+00,2.1431026690e-01
+1.1300000000e+01,-1.0095896955e+00,3.1678739595e-01
+1.1400000000e+01,-9.7340792869e-01,4.1620036208e-01
+1.1500000000e+01,-9.2745953431e-01,5.1155234283e-01
+1.1600000000e+01,-8.7219405299e-01,6.0188611310e-01
+1.1700000000e+01,-8.0815461103e-01,6.8629365342e-01
+1.1800000000e+01,-7.3597258484e-01,7.6392527596e-01
+1.1900000000e+01,-6.5636137778e-01,8.3399817662e-01
+1.2000000000e+01,-5.7010937090e-01,8.9580432700e-01
+1.2100000000e+01,-4.7807211709e-01,9.4871762686e-01
+1.2200000000e+01,-3.8116385594e-01,9.9220024475e-01
+1.2300000000e+01,-2.8034843350e-01,1.0258080826e+00
+1.2400000000e+01,-1.7662971708e-01,1.0491953088e+00
+1.2500000000e+01,-7.1041600360e-02,1.0621179136e+00
+1.2600000000e+01,3.5362301644e-02,1.0644362501e+00
+1.2700000000e+01,1.4151916502e-01,1.0561165344e+00
+1.2800000000e+01,2.4636757201e-01,1.0372312899e+00
+1.2900000000e+01,3.4885812263e-01,1.0079587278e+00
+1.3000000000e+01,4.4796393673e-01,9.6858107190e-01
+1.3100000000e+01,5.4269094151e-01,9.1948184161e-01
+1.3200000000e+01,6.3208784132e-01,8.6114211912e-01
+1.3300000000e+01,7.1525566985e-01,7.9413583829e-01
+1.3400000000e+01,7.9135682876e-01,7.1912414021e-01
+1.3500000000e+01,8.5962352196e-01,6.3684885159e-01
+1.3600000000e+01,9.1936550048e-01,5.4812514997e-01
+1.3700000000e+01,9.6997703999e-01,4.5383348839e-01
+1.3800000000e+01,1.0109430803e+00,3.5491085945e-01
+1.3900000000e+01,1.0418444648e+00,2.5234148548e-01
+1.4000000000e+01,1.0623622263e+00,1.4714702726e-01
+1.4100000000e+01,1.0722808758e+00,4.0376408683e-02
+1.4200000000e+01,1.0714906602e+00,-6.6904641270e-02
+1.4300000000e+01,1.0599887648e+00,-1.7362422501e-01
+1.4400000000e+01,1.0378794492e+00,-2.7871498511e-01
+1.4500000000e+01,1.0053731112e+00,-3.8112477846e-01
+1.4600000000e+01,9.6278429084e-01,-4.7982720892e-01
+1.4700000000e+01,9.1052862946e-01,-5.7383191325e-01
+1.4800000000e+01,8.4911881659e-01,-6.6219449702e-01
+1.4900000000e+01,7.7915956184e-01,-7.4402602105e-01
+1.5000000000e+01,7.0134164215e-01,-8.1850194292e-01
+1.5100000000e+01,6.1643508275e-01,-8.8487042375e-01
+1.5200000000e+01,5.2528153904e-01,-9.4245991649e-01
+1.5300000000e+01,4.2878595488e-01,-9.9068595919e-01
+1.5400000000e+01,3.2790757988e-01,-1.0290571046e+00
+1.5500000000e+01,2.2365043477e-01,-1.0571799258e+00
+1.5600000000e+01,1.1705331981e-01,-1.0747630470e+00
+1.5700000000e+01,9.1794654397e-03,-1.0816201586e+00
+1.5800000000e+01,-9.8894071577e-02,-1.0776719843e+00
+1.5900000000e+01,-2.0608716086e-01,-1.0629471809e+00
+1.6000000000e+01,-3.1132739259e-01,-1.0375821597e+00
+1.6100000000e+01,-4.1356080530e-01,-1.0018198309e+00
+1.6200000000e+01,-5.1176244003e-01,-9.5600728162e-01
+1.6300000000e+01,-6.0494661502e-01,-9.0059241086e-01
+1.6400000000e+01,-6.9217681820e-01,-8.3611955291e-01
+1.6500000000e+01,-7.7257511796e-01,-7.6322413256e-01
+1.6600000000e+01,-8.4533099777e-01,-6.8262640432e-01
+1.6700000000e+01,-9.0970952576e-01,-5.9512433751e-01
+1.6800000000e+01,-9.6505877692e-01,-5.0158571728e-01
+1.6900000000e+01,-1.0108164329e+00,-4.0293953982e-01
+1.7000000000e+01,-1.0465154929e+00,-3.0016678718e-01
+1.7100000000e+01,-1.0717890369e+00,-1.9429067313e-01
+1.7200000000e+01,-1.0863739941e+00,-8.6366457058e-02
+1.7300000000e+01,-1.0901138762e+00,2.2529072551e-02
+1.7400000000e+01,-1.0829604484e+00,1.3130834045e-01
+1.7500000000e+01,-1.0649743200e+00,2.3888384552e-01
+1.7600000000e+01,-1.0363244473e+00,3.4417903174e-01
+1.7700000000e+01,-9.9728655298e-01,4.4613906071e-01
+1.7800000000e+01,-9.4824047696e-01,5.4374137777e-01
+1.7900000000e+01,-8.8966648340e-01,6.3600596605e-01
+1.8000000000e+01,-8.2214056043e-01,7.2200518546e-01
+1.8100000000e+01,-7.4632875811e-01,8.0087309808e-01
+1.8200000000e+01,-6.6298062061e-01,8.7181418597e-01
+1.8300000000e+01,-5.7292177708e-01,9.3411137385e-01
+1.8400000000e+01,-4.7704576454e-01,9.8713327564e-01
+1.8500000000e+01,-3.7630516385e-01,1.0303405917e+00
+1.8600000000e+01,-2.7170213651e-01,1.0632915917e+00
+1.8700000000e+01,-1.6427845644e-01,1.0856466285e+00
+1.8800000000e+01,-5.5105135725e-02,1.0971716347e+00
+1.8900000000e+01,5.4728252195e-02,1.0977405685e+00
+1.9000000000e+01,1.6412444562e-01,1.0873367818e+00
+1.9100000000e+01,2.7198945540e-01,1.0660532957e+00
+1.9200000000e+01,3.7724350401e-01,1.0340919806e+00
+1.9300000000e+01,4.7883183346e-01,9.9176164743e-01
+1.9400000000e+01,5.7573527399e-01,9.3947506691e-01
+1.9500000000e+01,6.6698046754e-01,8.7774494814e-01
+1.9600000000e+01,7.5164964349e-01,8.0717891305e-01
+1.9700000000e+01,8.2888984858e-01,7.2847351751e-01
+1.9800000000e+01,8.9792153804e-01,6.4240737736e-01
+1.9900000000e+01,9.5804644160e-01,5.4983346724e-01
+2.0000000000e+01,1.0086546249e+00,4.5167066844e-01
+2.0100000000e+01,1.0492306748e+00,3.4889464950e-01
+2.0200000000e+01,1.0793589463e+00,2.4252816994e-01
+2.0300000000e+01,1.0987278172e+00,1.3363090369e-01
+2.0400000000e+01,1.1071329071e+00,2.3288883149e-02
+2.0500000000e+01,1.1044792280e+00,-8.7396330585e-02
+2.0600000000e+01,1.0907822430e+00,-1.9731864374e-01
+2.0700000000e+01,1.0661678222e+00,-3.0537848208e-01
+2.0800000000e+01,1.0308710952e+00,-4.1049378909e-01
+2.0900000000e+01,9.8523420991e-01,-5.1161085996e-01
+2.1000000000e+01,9.2970301917e-01,-6.0771490308e-01
+2.1100000000e+01,8.6482272793e-01,-6.9784022307e-01
+2.1200000000e+01,7.9123254253e-01,-7.8107992324e-01
+2.1300000000e+01,7.0965937510e-01,-8.5659502998e-01
+2.1400000000e+01,6.2091066487e-01,-9.2362294728e-01
+2.1500000000e+01,5.2586638741e-01,-9.8148515628e-01
+2.1600000000e+01,4.2547033073e-01,-1.0295940820e+00
+2.1700000000e+01,3.2072072479e-01,-1.0674590579e+00
+2.1800000000e+01,2.1266031737e-01,-1.0946913277e+00
+2.1900000000e+01,1.0236599509e-01,-1.1110080332e+00
+2.2000000000e+01,-9.0619472952e-03,-1.1162351481e+00
+2.2100000000e+01,-1.2051078228e-01,-1.1103093263e+00
+2.2200000000e+01,-2.3086646029e-01,-1.0932786457e+00
+2.2300000000e+01,-3.3902474468e-01,-1.0653022396e+00
+2.2400000000e+01,-4.4390225975e-01,-1.0266488164e+00
+2.2500000000e+01,-5.4444734133e-01,-9.7769408359e-01
+2.2600000000e+01,-6.3965058123e-01,-9.1891709900e-01
+2.2700000000e+01,-7.2855495975e-01,-8.5089558542e-01
+2.2800000000e+01,-8.1026546475e-01,-7.7430025424e-01
+2.2900000000e+01,-8.8395810046e-01,-6.8988819393e-01
+2.3000000000e+01,-9.4888819532e-01,-5.9849538833e-01
+2.3100000000e+01,-1.0043979254e+00,-5.0102843894e-01
+2.3200000000e+01,-1.0499229768e+00,-3.9845557292e-01
+2.3300000000e+01,-1.0849982807e+00,-2.9179702613e-01
+2.3400000000e+01,-1.1092627611e+00,-1.8211489661e-01
+2.3500000000e+01,-1.1224630483e+00,-7.0502569448e-02
+2.3600000000e+01,-1.1244561199e+00,4.1926181624e-02
+2.3700000000e+01,-1.1152108405e+00,1.5404830744e-01
+2.3800000000e+01,-1.0948083846e+00,2.6474270037e-01
+2.3900000000e+01,-1.0634415381e+00,3.7290140316e-01
+2.4000000000e+01,-1.0214128826e+00,4.7744069748e-01
+2.4100000000e+01,-9.6913188053e-01,5.7731196133e-01
+2.4200000000e+01,-9.0711088879e-01,6.7151218638e-01
+2.4300000000e+01,-8.3596014044e-01,7.5909405001e-01
+2.4400000000e+01,-7.5638174230e-01,8.3917544061e-01
+2.4500000000e+01,-6.6916274855e-01,9.1094834071e-01
+2.4600000000e+01,-5.7516737806e-01,9.7368697816e-01
+2.4700000000e+01,-4.7532845268e-01,1.0267551635e+00
+2.4800000000e+01,-3.7063814112e-01,1.0696127391e+00
+2.4900000000e+01,-2.6213810042e-01,1.1018210746e+00
+2.5000000000e+01,-1.5090911271e-01,1.1230475534e+00
+2.5100000000e+01,-3.8060320643e-02,1.1330690036e+00
+2.5200000000e+01,7.5281831510e-02,1.1317740380e+00
+2.5300000000e+01,1.8798484078e-01,1.1191642801e+00
+2.5400000000e+01,2.9892146023e-01,1.0953544603e+00
+2.5500000000e+01,4.0698097242e-01,1.0605713822e+00
+2.5600000000e+01,5.1108030900e-01,1.0151517682e+00
+2.5700000000e+01,6.1017490496e-01,9.5953900270e-01
+2.5800000000e+01,7.0326917890e-01,8.9427880713e-01
+2.5900000000e+01,7.8942653410e-01,8.2001388699e-01
+2.6000000000e+01,8.6777878002e-01,7.3747760424e-01
+2.6100000000e+01,9.3753487936e-01,6.4748673715e-01
+2.6200000000e+01,9.9798893265e-01,5.5093339904e-01
+2.6300000000e+01,1.0485273198e+00,4.4877619585e-01
+2.6400000000e+01,1.0886349263e+00,3.4203071018e-01
+2.6500000000e+01,1.1179003910e+00,2.3175940631e-01
+2.6600000000e+01,1.1360203219e+00,1.1906105640e-01
+2.6700000000e+01,1.1428024372e+00,5.0597932721e-03
+2.6800000000e+01,1.1381675982e+00,-1.0910610140e-01
+2.6900000000e+01,1.1221507137e+00,-2.2229556143e-01
+2.7000000000e+01,1.0949005055e+00,-3.3337613974e-01
+2.7100000000e+01,1.0566781351e+00,-4.4123533721e-01
+2.7200000000e+01,1.0078547057e+00,-5.4479174345e-01
+2.7300000000e+01,9.4890766270e-01,-6.4300587795e-01
+2.7400000000e+01,8.8041612743e-01,-7.3489062294e-01
+2.7500000000e+01,8.0305520961e-01,-8.1952114324e-01
+2.7600000000e+01,7.1758935498e-01,-8.9604419349e-01
+2.7700000000e+01,6.2486479304e-01,-9.6368671898e-01
+2.7800000000e+01,5.2580115978e-01,-1.0217636636e+00
+2.7900000000e+01,4.2138237817e-01,-1.0696849057e+00
+2.8000000000e+01,3.1264688695e-01,-1.1069612522e+00
+2.8100000000e+01,2.0067731468e-01,-1.1332094299e+00
+2.8200000000e+01,8.6589701730e-02,-1.1481560223e+00
+2.8300000000e+01,-2.8477622542e-02,-1.1516403126e+00
+2.8400000000e+01,-1.4337539648e-01,-1.1436160038e+00
+2.8500000000e+01,-2.5695490353e-01,-1.1241517949e+00
+2.8600000000e+01,-3.6807945562e-01,-1.0934308101e+00
+2.8700000000e+01,-4.7563576782e-01,-1.0517488821e+00
+2.8800000000e+01,-5.7854511045e-01,-9.9951170809e-01
+2.8900000000e+01,-6.7577412686e-01,-9.3723090396e-01
+2.9000000000e+01,-7.6634520858e-01,-8.6551899606e-01
+2.9100000000e+01,-8.4934632344e-01,-7.8508339924e-01
+2.9200000000e+01,-9.2394019814e-01,-6.9671944027e-01
+2.9300000000e+01,-9.8937276258e-01,-6.0130249548e-01
+2.9400000000e+01,-1.0449807710e+00,-4.9977932040e-01
+2.9500000000e+01,-1.0901985229e+00,-3.9315865715e-01
+2.9600000000e+01,-1.1245636153e+00,-2.8250121281e-01
+2.9700000000e+01,-1.1477216687e+00,-1.6890910831e-01
+2.9800000000e+01,-1.1594299771e+00,-5.3514902805e-02
+2.9900000000e+01,-1.1595600466e+00,6.2529697218e-02
+3.0000000000e+01,-1.1480989935e+00,1.7806533156e-01
+3.0100000000e+01,-1.1251497890e+00,2.9193656787e-01
+3.0200000000e+01,-1.0909303458e+00,4.0300345559e-01
+3.0300000000e+01,-1.0457714546e+00,5.1015293669e-01
+3.0400000000e+01,-9.9011359083e-01,6.1230999893e-01
+3.0500000000e+01,-9.2450262057e-01,7.0844846006e-01
+3.0600000000e+01,-8.4958445000e-01,7.9760127462e-01
+3.0700000000e+01,-7.6609866852e-01,8.7887025996e-01
+3.0800000000e+01,-6.7487124944e-01,9.5143514355e-01
+3.0900000000e+01,-5.7680637980e-01,1.0145618406e+00
+3.1000000000e+01,-4.7287750005e-01,1.0676098784e+00
+3.1100000000e+01,-3.6411764268e-01,1.1100388924e+00
+3.1200000000e+01,-2.5160916532e-01,1.1414141275e+00
+3.1300000000e+01,-1.3647298058e-01,1.1614108905e+00
+3.1400000000e+01,-1.9857389729e-02,1.1698179049e+00
+3.1500000000e+01,9.7073368649e-02,1.1665395373e+00
+3.1600000000e+01,2.1315074070e-01,1.1515968684e+00
+3.1700000000e+01,3.2721353417e-01,1.1251275992e+00
+3.1800000000e+01,4.3811953331e-01,1.0873847915e+00
+3.1900000000e+01,5.4475693575e-01,1.0387344532e+00
+3.2000000000e+01,6.4605549677e-01,9.7965199282e-01
+3.2100000000e+01,7.4099726935e-01,9.1071757661e-01
+3.2200000000e+01,8.2862683207e-01,8.3261043390e-01
+3.2300000000e+01,9.0806090239e-01,7.4610216639e-01
+3.2400000000e+01,9.7849723825e-01,6.5204912748e-01
+3.2500000000e+01,1.0392227386e+00,5.5138394676e-01
+3.2600000000e+01,1.0896206609e+00,4.4510628348e-01
+3.2700000000e+01,1.1291768827e+00,3.3427290076e-01
+3.2800000000e+01,1.1574851435e+00,2.1998715911e-01
+3.2900000000e+01,1.1742512139e+00,1.0338803351e-01
+3.3000000000e+01,1.1792959485e+00,-1.4361236374e-02
+3.3100000000e+01,1.1725571923e+00,-1.3208475002e-01
+3.3200000000e+01,1.1540905189e+00,-2.4860568791e-01
+3.3300000000e+01,1.1240687926e+00,-3.6275807445e-01
+3.3400000000e+01,1.0827805577e+00,-4.7339844415e-01
+3.3500000000e+01,1.0306272703e+00,-5.7941729458e-01
+3.3600000000e+01,9.6811939767e-01,-6.7975021113e-01
+3.3700000000e+01,8.9587142501e-01,-7.7338855235e-01
+3.3800000000e+01,8.1459581664e-01,-8.5938958832e-01
+3.3900000000e+01,7.2509599221e-01,-9.3688599041e-01
+3.4000000000e+01,6.2825838646e-01,-1.0050945767e+00
+3.4100000000e+01,5.2504367129e-01,-1.0633242251e+00
+3.4200000000e+01,4.1647722686e-01,-1.1109828739e+00
+3.4300000000e+01,3.0363895641e-01,-1.1475835397e+00
+3.4400000000e+01,1.8765254582e-01,-1.1727492901e+00
+3.4500000000e+01,6.9674274745e-02,-1.1862171218e+00
+3.4600000000e+01,-4.9118509209e-02,-1.1878407039e+00
+3.4700000000e+01,-1.6753913849e-01,-1.1775919570e+00
+3.4800000000e+01,-2.8440347851e-01,-1.1555614525e+00
+3.4900000000e+01,-3.9854176704e-01,-1.1219576258e+00
+3.5000000000e+01,-5.0881032188e-01,-1.0771048107e+00
+3.5100000000e+01,-6.1410299946e-01,-1.0214401129e+00
+3.5200000000e+01,-7.1336228967e-01,-9.5550915417e-01
+3.5300000000e+01,-8.0558993567e-01,-8.7996072588e-01
+3.5400000000e+01,-8.8985697200e-01,-7.9554040727e-01
+3.5500000000e+01,-9.6531308007e-01,-7.0308320912e-01
+3.5600000000e+01,-1.0311951668e+00,-6.0350531624e-01
+3.5700000000e+01,-1.0868350800e+00,-4.9779501031e-01
+3.5800000000e+01,-1.1316663825e+00,-3.8700286284e-01
+3.5900000000e+01,-1.1652301157e+00,-2.7223129579e-01
+3.6000000000e+01,-1.1871794956e+00,-1.5462361326e-01
+3.6100000000e+01,-1.1972834911e+00,-3.5352613744e-02
+3.6200000000e+01,-1.1954292494e+00,8.4391104046e-02
+3.6300000000e+01,-1.1816233423e+00,2.0341102483e-01
+3.6400000000e+01,-1.1559918201e+00,3.2051667155e-01
+3.6500000000e+01,-1.1187790711e+00,4.3453551167e-01
+3.6600000000e+01,-1.0703454966e+00,5.4432469611e-01
+3.6700000000e+01,-1.0111640250e+00,6.4878251312e-01
+3.6800000000e+01,-9.4181549636e-01,7.4685944241e-01
+3.6900000000e+01,-8.6298296431e-01,8.3756869868e-01
+3.7000000000e+01,-7.7544497032e-01,9.1999615847e-01
+3.7100000000e+01,-6.8006785674e-01,9.9330957074e-01
+3.7200000000e+01,-5.7779719443e-01,1.0567669581e+00
+3.7300000000e+01,-4.6964840990e-01,1.1097241241e+00
+3.7400000000e+01,-3.5669670474e-01,1.1516411907e+00
+3.7500000000e+01,-2.4006636741e-01,1.1820880994e+00
+3.7600000000e+01,-1.2091958369e-01,1.2007490200e+00
+3.7700000000e+01,-4.4485709373e-04,1.2074256218e+00
+3.7800000000e+01,1.2015484536e-01,1.2020391741e+00
+3.7900000000e+01,2.3967410385e-01,1.1846314529e+00
+3.8000000000e+01,3.5691709302e-01,1.1553644433e+00
+3.8100000000e+01,4.7070954526e-01,1.1145188411e+00
+3.8200000000e+01,5.7991051063e-01,1.0624913653e+00
+3.8300000000e+01,6.8342379537e-01,9.9979090877e-01
+3.8400000000e+01,7.8020896449e-01,9.2703356313e-01
+3.8500000000e+01,8.6929179807e-01,8.4493656763e-01
+3.8600000000e+01,9.4977409598e-01,7.5431124057e-01
+3.8700000000e+01,1.0208427326e+00,6.5605496324e-01
+3.8800000000e+01,1.0817778701e+00,5.5114229563e-01
+3.8900000000e+01,1.1319602476e+00,4.4061531166e-01
+3.9000000000e+01,1.1708774718e+00,3.2557324990e-01
+3.9100000000e+01,1.1981292462e+00,2.0716158253e-01
+3.9200000000e+01,1.2134314844e+00,8.6560611108e-02
+3.9300000000e+01,1.2166192667e+00,-3.5026297127e-02
+3.9400000000e+01,1.2076486072e+00,-1.5638471536e-01
+3.9500000000e+01,1.1865970150e+00,-2.7630128584e-01
+3.9600000000e+01,1.1536628407e+00,-3.9357584998e-01
+3.9700000000e+01,1.1091634145e+00,-5.0703345868e-01
+3.9800000000e+01,1.0535319935e+00,-6.1553614257e-01
+3.9900000000e+01,9.8731354631e-01,-7.1799432449e-01
+4.0000000000e+01,9.1115941729e-01,-8.1337775955e-01
+4.0100000000e+01,8.2582092203e-01,-9.0072589311e-01
+4.0200000000e+01,7.3214193703e-01,-9.7915753243e-01
+4.0300000000e+01,6.3105055665e-01,-1.0478797347e+00
+4.0400000000e+01,5.2354989991e-01,-1.1061958217e+00
+4.0500000000e+01,4.1070815816e-01,-1.1535124405e+00
+4.0600000000e+01,2.9364798223e-01,-1.1893455978e+00
+4.0700000000e+01,1.7353531463e-01,-1.2133256079e+00
+4.0800000000e+01,5.1567777810e-02,-1.2252009020e+00
+4.0900000000e+01,-7.1037265870e-02,-1.2248406615e+00
+4.1000000000e+01,-1.9305486347e-01,-1.2122362465e+00
+4.1100000000e+01,-3.1326470881e-01,-1.1875014042e+00
+4.1200000000e+01,-4.3046334559e-01,-1.1508712544e+00
+4.1300000000e+01,-5.4347621431e-01,-1.1027000606e+00
+4.1400000000e+01,-6.5116942237e-01,-1.0434578077e+00
+4.1500000000e+01,-7.5246111942e-01,-9.7372562033e-01
+4.1600000000e+01,-8.4633236388e-01,-8.9419006409e-01
+4.1700000000e+01,-9.3183737150e-01,-8.0563638858e-01
+4.1800000000e+01,-1.0081130428e+00,-7.0894077635e-01
+4.1900000000e+01,-1.0743876735e+00,-6.0506167537e-01
+4.2000000000e+01,-1.1299887604e+00,-4.9503030043e-01
+4.2100000000e+01,-1.1743498226e+00,-3.7994039768e-01
+4.2200000000e+01,-1.2070161707e+00,-2.6093737381e-01
+4.2300000000e+01,-1.2276495639e+00,-1.3920689799e-01
+4.2400000000e+01,-1.2360317072e+00,-1.5963089878e-02
+4.2500000000e+01,-1.2320665537e+00,1.0756358857e-01
+4.2600000000e+01,-1.2157813874e+00,2.3013861671e-01
+4.2700000000e+01,-1.1873266732e+00,3.5053575179e-01
+4.2800000000e+01,-1.1469746766e+00,4.6754929499e-01
+4.2900000000e+01,-1.0951168636e+00,5.8000616452e-01
+4.3000000000e+01,-1.0322601069e+00,6.8677765459e-01
+4.3100000000e+01,-9.5902173446e-01,7.8679076272e-01
+4.3200000000e+01,-8.7612346920e-01,8.7903897139e-01
+4.3300000000e+01,-7.8438431865e-01,9.6259237600e-01
+4.3400000000e+01,-6.8471248548e-01,1.0366070571e+00
+4.3500000000e+01,-5.7809637807e-01,1.1003336023e+00
+4.3600000000e+01,-4.6559481033e-01,1.1531246923e+00
+4.3700000000e+01,-3.4832648801e-01,1.1944416739e+00
+4.3800000000e+01,-2.2745888573e-01,1.2238600529e+00
+4.3900000000e+01,-1.0419662544e-01,1.2410738517e+00
+4.4000000000e+01,2.0230528162e-02,1.2458987870e+00
+4.4100000000e+01,1.4457993550e-01,1.2382742329e+00
+4.4200000000e+01,2.6760849096e-01,1.2182639513e+00
+4.4300000000e+01,3.8808504883e-01,1.1860555788e+00
+4.4400000000e+01,5.0480274211e-01,1.1419588747e+00
+4.4500000000e+01,6.1659107094e-01,1.0864027471e+00
+4.4600000000e+01,7.2232763959e-01,1.0199310842e+00
+4.4700000000e+01,8.2094942434e-01,9.4319743315e-01
+4.4800000000e+01,9.1146345914e-01,8.5695857603e-01
+4.4900000000e+01,9.9295683162e-01,7.6206706860e-01
+4.5000000000e+01,1.0646058889e+00,6.5946281361e-01
+4.5100000000e+01,1.1256845603e+00,5.5016375267e-01
+4.5200000000e+01,1.1755717129e+00,4.3525576877e-01
+4.5300000000e+01,1.2137574658e+00,3.1588189944e-01
+4.5400000000e+01,1.2398483980e+00,1.9323096791e-01
+4.5500000000e+01,1.2535715975e+00,6.8525745064e-02
+4.5600000000e+01,1.2547775107e+00,-5.6989239649e-02
+4.5700000000e+01,1.2434415600e+00,-1.8206011713e-01
+4.5800000000e+01,1.2196645149e+00,-3.0543620323e-01
+4.5900000000e+01,1.1836716099e+00,-4.2588250415e-01
+4.6000000000e+01,1.1358104179e+00,-5.4219207737e-01
+4.6100000000e+01,1.0765474976e+00,-6.5319812467e-01
+4.6200000000e+01,1.0064638492e+00,-7.5778569588e-01
+4.6300000000e+01,9.2624922089e-01,-8.5490288606e-01
+4.6400000000e+01,8.3669532252e-01,-9.4357141372e-01
+4.6500000000e+01,7.3868801388e-01,-1.0228964736e+00
+4.6600000000e+01,6.3319854335e-01,-1.0920757647e+00
+4.6700000000e+01,5.2127392459e-01,-1.1504076034e+00
+4.6800000000e+01,4.0402654608e-01,-1.1972980375e+00
+4.6900000000e+01,2.8262311684e-01,-1.2322668920e+00
+4.7000000000e+01,1.5827305803e-01,-1.2549526825e+00
+4.7100000000e+01,3.2216456058e-02,-1.2651163476e+00
+4.7200000000e+01,-9.4288303027e-02,-1.2626437610e+00
+4.7300000000e+01,-2.1997709431e-01,-1.2475469973e+00
+4.7400000000e+01,-3.4359268482e-01,-1.2199643376e+00
+4.7500000000e+01,-4.6389730799e-01,-1.1801590138e+00
+4.7600000000e+01,-5.7968505623e-01,-1.1285167015e+00
+4.7700000000e+01,-6.8979396727e-01,-1.0655417880e+00
+4.7800000000e+01,-7.9311768347e-01,-9.9185244849e-01
+4.7900000000e+01,-8.8861656683e-01,-9.0817458094e-01
+4.8000000000e+01,-9.7532815821e-01,-8.1533465721e-01
+4.8100000000e+01,-1.0523768756e+00,-7.1425156227e-01
+4.8200000000e+01,-1.1189828535e+00,-6.0592750116e-01
+4.8300000000e+01,-1.1744698343e+00,-4.9143806413e-01
+4.8400000000e+01,-1.2182720324e+00,-3.7192154811e-01
+4.8500000000e+01,-1.2499399000e+00,-2.4856764080e-01
+4.8600000000e+01,-1.2691447373e+00,-1.2260557964e-01
+4.8700000000e+01,-1.2756820993e+00,4.7080963623e-03
+4.8800000000e+01,-1.2694739637e+00,1.3210208004e-01
+4.8900000000e+01,-1.2505696375e+00,2.5830298955e-01
+4.9000000000e+01,-1.2191453909e+00,3.8204809538e-01
+4.9100000000e+01,-1.1755028232e+00,5.0209795354e-01
+4.9200000000e+01,-1.1200659726e+00,6.1724881865e-01
+4.9300000000e+01,-1.0533772000e+00,7.2634471263e-01
+4.9400000000e+01,-9.7609188630e-01,8.2828902787e-01
+4.9500000000e+01,-8.8897199342e-01,9.2205554868e-01
+4.9600000000e+01,-7.9287855467e-01,1.0066987799e+00
+4.9700000000e+01,-6.8876316654e-01,1.0813634793e+00
+4.9800000000e+01,-5.7765856699e-01,1.1452932967e+00
+4.9900000000e+01,-4.6066839306e-01,1.1978384340e+00
+5.0000000000e+01,-3.3895621963e-01,1.2384622469e+00
+5.0100000000e+01,-2.1373398795e-01,1.2667467231e+00
+5.0200000000e+01,-8.6249939174e-02,1.2823967789e+00
+5.0300000000e+01,4.2223827032e-02,1.2852433322e+00
+5.0400000000e+01,1.7040404488e-01,1.2752451190e+00
+5.0500000000e+01,2.9700909912e-01,1.2524892335e+00
+5.0600000000e+01,4.2077183803e-01,1.2171903856e+00
+5.0700000000e+01,5.4045225467e-01,1.1696888829e+00
+5.0800000000e+01,6.5484990946e-01,1.1104473542e+00
+5.0900000000e+01,7.6281596963e-01,1.0400462474e+00
+5.1000000000e+01,8.6326474492e-01,9.5917814486e-01
+5.1100000000e+01,9.5518460352e-01,8.6864095280e-01
+5.1200000000e+01,1.0376481587e+00,7.6933003021e-01
+5.1300000000e+01,1.1098216237e+00,6.6222933626e-01
+5.1400000000e+01,1.1709732396e+00,5.4840168305e-01
+5.1500000000e+01,1.2204806927e+00,4.2897819034e-01
+5.1600000000e+01,1.2578374448e+00,3.0514704690e-01
+5.1700000000e+01,1.2826579122e+00,1.7814169025e-01
+5.1800000000e+01,1.2946814414e+00,4.9228522113e-02
+5.1900000000e+01,1.2937750395e+00,-8.0305717863e-02
+5.2000000000e+01,1.2799348314e+00,-2.0916679431e-01
+5.2100000000e+01,1.2532862276e+00,-3.3606590615e-01
+5.2200000000e+01,1.2140828000e+00,-4.5973257523e-01
+5.2300000000e+01,1.1627038755e+00,-5.7892736464e-01
+5.2400000000e+01,1.0996508697e+00,-6.9245429952e-01
+5.2500000000e+01,1.0255423970e+00,-7.9917286576e-01
+5.2600000000e+01,9.4110820365e-01,-8.9800946641e-01
+5.2700000000e+01,8.4718198473e-01,-9.8796822054e-01
+5.2800000000e+01,7.4469315394e-01,-1.0681409961e+00
+5.2900000000e+01,6.3465764885e-01,-1.1377165754e+00
+5.3000000000e+01,5.1816786185e-01,-1.1959888618e+00
+5.3100000000e+01,3.9638179687e-01,-1.2423640423e+00
+5.3200000000e+01,2.7051155925e-01,-1.2763666367e+00
+5.3300000000e+01,1.4181129341e-01,-1.2976443683e+00
+5.3400000000e+01,1.1564688166e-02,-1.3059718090e+00
+5.3500000000e+01,-1.1892782583e-01,-1.3012527594e+00
+5.3600000000e+01,-2.4836205874e-01,-1.2835213401e+00
+5.3700000000e+01,-3.7544309478e-01,-1.2529417805e+00
+5.3800000000e+01,-4.9889824582e-01,-1.2098069077e+00
+5.3900000000e+01,-6.1748979572e-01,-1.1545353476e+00
+5.4000000000e+01,-7.3002740784e-01,-1.0876674664e+00
+5.4100000000e+01,-8.3538007137e-01,-1.0098600908e+00
+5.4200000000e+01,-9.3248746648e-01,-9.2188005831e-01
+5.4300000000e+01,-1.0203706344e+00,-8.2459666215e-01
+5.4400000000e+01,-1.0981418445e+00,-7.1897306398e-01
+5.4500000000e+01,-1.1650135598e+00,-6.0605675979e-01
+5.4600000000e+01,-1.2203064095e+00,-4.8696919325e-01
+5.4700000000e+01,-1.2634560883e+00,-3.6289461949e-01
+5.4800000000e+01,-1.2940191127e+00,-2.3506832976e-01
+5.4900000000e+01,-1.3116773745e+00,-1.0476435429e-01
+5.5000000000e+01,-1.3162414469e+00,2.6717234541e-02
+5.5100000000e+01,-1.3076526062e+00,1.5806329197e-01
+5.5200000000e+01,-1.2859835497e+00,2.8796071432e-01
+5.5300000000e+01,-1.2514378006e+00,4.1510956511e-01
+5.5400000000e+01,-1.2043478044e+00,5.3823608276e-01
+5.5500000000e+01,-1.1451717347e+00,6.5610543963e-01
+5.5600000000e+01,-1.0744890386e+00,7.6753412455e-01
+5.5700000000e+01,-9.9299476647e-01,8.7140182491e-01
+5.5800000000e+01,-9.0149273887e-01,9.6666268875e-01
+5.5900000000e+01,-8.0088762078e-01,1.0523558539e+00
+5.6000000000e+01,-6.9217597898e-01,1.1276151380e+00
+5.6100000000e+01,-5.7643641209e-01,1.1916777919e+00
+5.6200000000e+01,-4.5481885075e-01,1.2438922275e+00
+5.6300000000e+01,-3.2853313392e-01,1.2837246426e+00
+5.6400000000e+01,-1.9883697501e-01,1.3107644750e+00
+5.6500000000e+01,-6.7023437350e-02,1.3247286293e+00
+5.6600000000e+01,6.5591956256e-02,1.3254644347e+00
+5.6700000000e+01,1.9768435132e-01,1.3129513017e+00
+5.6800000000e+01,3.2793279593e-01,1.2873010598e+00
+5.6900000000e+01,4.5503344923e-01,1.2487569721e+00
+5.7000000000e+01,5.7771263202e-01,1.1976914355e+00
+5.7100000000e+01,6.9473958901e-01,1.1346023867e+00
+5.7200000000e+01,8.0493883465e-01,1.0601084507e+00
+5.7300000000e+01,9.0720195887e-01,9.7494287763e-01
+5.7400000000e+01,1.0004987742e+00,8.7994632734e-01
+5.7500000000e+01,1.0838876919e+00,7.7605857379e-01
+5.7600000000e+01,1.1565252231e+00,6.6430920980e-01
+5.7700000000e+01,1.2176745087e+00,5.4580744484e-01
+5.7800000000e+01,1.2667127916e+00,4.2173109680e-01
+5.7900000000e+01,1.3031377564e+00,2.9331488687e-01
+5.8000000000e+01,1.3265726711e+00,1.6183815383e-01
+5.8100000000e+01,1.3367702785e+00,2.8612110144e-02
+5.8200000000e+01,1.3336153975e+00,-1.0503323338e-01
+5.8300000000e+01,1.3171262065e+00,-2.3776234388e-01
+5.8400000000e+01,1.2874541954e+00,-3.6824751128e-01
+5.8500000000e+01,1.2448827842e+00,-4.9518212820e-01
+5.8600000000e+01,1.1898246222e+00,-6.1729377225e-01
+5.8700000000e+01,1.1228175922e+00,-7.3335695971e-01
+5.8800000000e+01,1.0445195592e+00,-8.4220544270e-01
+5.8900000000e+01,9.5570191357e-01,-9.4274392655e-01
+5.9000000000e+01,8.5724197459e-01,-1.0339590897e+00
+5.9100000000e+01,7.5011432580e-01,-1.1149297949e+00
+5.9200000000e+01,6.3538117041e-01,-1.1848363897e+00
+5.9300000000e+01,5.1418180086e-01,-1.2429690011e+00
+5.9400000000e+01,3.8772128727e-01,-1.2887347418e+00
+5.9500000000e+01,2.5725849712e-01,-1.3216637535e+00
+5.9600000000e+01,1.2409356507e-01,-1.3414140270e+00
+5.9700000000e+01,-1.0445062223e-02,-1.3477749490e+00
+5.9800000000e+01,-1.4501387724e-01,-1.3406695387e+00
+5.9900000000e+01,-2.7826772660e-01,-1.3201553512e+00
+6.0000000000e+01,-4.0887325588e-01,-1.2864240366e+00
+6.0100000000e+01,-5.3552224996e-01,-1.2397995585e+00
+6.0200000000e+01,-6.5694473547e-01,-1.1807350880e+00
+6.0300000000e+01,-7.7192171395e-01,-1.1098086027e+00
+6.0400000000e+01,-8.7929739832e-01,-1.0277172341e+00
+6.0500000000e+01,-9.7799082961e-01,-9.3527041825e-01
+6.0600000000e+01,-1.0670067574e+00,-8.3338191653e-01
+6.0700000000e+01,-1.1454456743e+00,-7.2306078602e-01
+6.0800000000e+01,-1.2125129035e+00,-6.0540138794e-01
+6.0900000000e+01,-1.2675266473e+00,-4.8157253356e-01
+6.1000000000e+01,-1.3099249157e+00,-3.5280587488e-01
+6.1100000000e+01,-1.3392712639e+00,-2.2038365567e-01
+6.1200000000e+01,-1.3552592804e+00,-8.5625944401e-02
+6.1300000000e+01,-1.3577157809e+00,5.0122523634e-02
+6.1400000000e+01,-1.3466026719e+00,1.8550576394e-01
+6.1500000000e+01,-1.3220174671e+00,3.1917008632e-01
+6.1600000000e+01,-1.2841924482e+00,4.4977762905e-01
+6.1700000000e+01,-1.2334924778e+00,5.7601974824e-01
+6.1800000000e+01,-1.1704114847e+00,6.9663012866e-01
+6.1900000000e+01,-1.0955676559e+00,8.1039748444e-01
+6.2000000000e+01,-1.0096973807e+00,9.1617772250e-01
+6.2100000000e+01,-9.1364800863e-01,1.0129054464e+00
+6.2200000000e+01,-8.0836948956e-01,1.0996046849e+00
+6.2300000000e+01,-6.9490498063e-01,1.1753987377e+00
+6.2400000000e+01,-5.7438051183e-01,1.2395190385e+00
+6.2500000000e+01,-4.4799381281e-01,1.2913129455e+00
+6.2600000000e+01,-3.1700241181e-01,1.3302503821e+00
+6.2700000000e+01,-1.8271112491e-01,1.3559292575e+00
+6.2800000000e+01,-4.6459060101e-02,1.3680796147e+00
+6.2900000000e+01,9.0393734508e-02,1.3665664612e+00
+6.3000000000e+01,2.2647984509e-01,1.3513912532e+00
+6.3100000000e+01,3.6043815377e-01,1.3226920175e+00
+6.3200000000e+01,4.9092745097e-01,1.2807421090e+00
+6.3300000000e+01,6.1663986240e-01,1.2259476126e+00
+6.3400000000e+01,7.3631395627e-01,1.1588434176e+00
+6.3500000000e+01,8.4874739932e-01,1.0800879990e+00
+6.3600000000e+01,9.5280903477e-01,9.9045695934e-01
+6.3700000000e+01,1.0474502609e+00,8.9083539237e-01
+6.3800000000e+01,1.1317155956e+00,7.8220914522e-01
+6.3900000000e+01,1.2047523209e+00,6.6565506434e-01
+6.4000000000e+01,1.2658191098e+00,5.4233032212e-01
+6.4100000000e+01,1.3142935480e+00,4.1346092987e-01
+6.4200000000e+01,1.3496784748e+00,2.8032955118e-01
+6.4300000000e+01,1.3716070784e+00,1.4426273675e-01
+6.4400000000e+01,1.3798466924e+00,6.6177078084e-03
+6.4500000000e+01,1.3743012565e+00,-1.3123118039e-01
+6.4600000000e+01,1.3550124130e+00,-2.6790615976e-01
+6.4700000000e+01,1.3221592284e+00,-4.0203981854e-01
+6.4800000000e+01,1.2760565409e+00,-5.3228878060e-01
+6.4900000000e+01,1.2171519491e+00,-6.5734715824e-01
+6.5000000000e+01,1.1460214703e+00,-7.7595964367e-01
+6.5100000000e+01,1.0633639114e+00,-8.8693410796e-01
+6.5200000000e+01,9.6999400556e-01,-9.8915358111e-01
+6.5300000000e+01,8.6683438471e-01,-1.0815874927e+00
+6.5400000000e+01,7.5490646412e-01,-1.1633020603e+00
+6.5500000000e+01,6.3532033108e-01,-1.2334697204e+00
+6.5600000000e+01,5.0926373688e-01,-1.2913775073e+00
+6.5700000000e+01,3.7799030143e-01,-1.3364342951e+00
+6.5800000000e+01,2.4280704771e-01,-1.3681768278e+00
+6.5900000000e+01,1.0506138980e-01,-1.3862744790e+00
+6.6000000000e+01,-3.3872295742e-02,-1.3905326890e+00
+6.6100000000e+01,-1.7260638155e-01,-1.3808950460e+00
+6.6200000000e+01,-3.0975384726e-01,-1.3574439881e+00
+6.6300000000e+01,-4.4394214510e-01,-1.3204001184e+00
+6.6400000000e+01,-5.7382693474e-01,-1.2701201378e+00
+6.6500000000e+01,-6.9810554982e-01,-1.2070934155e+00
+6.6600000000e+01,-8.1553006148e-01,-1.1319372299e+00
+6.6700000000e+01,-9.2491980768e-01,-1.0453907261e+00
+6.6800000000e+01,-1.0251732627e+00,-9.4830764920e-01
+6.6900000000e+01,-1.1152791276e+00,-8.4164792471e-01
+6.7000000000e+01,-1.1943265293e+00,-7.2646816896e-01
+6.7100000000e+01,-1.2615142268e+00,-6.0391122379e-01
+6.7200000000e+01,-1.3161587299e+00,-4.7519481861e-01
+6.7300000000e+01,-1.3577012498e+00,-3.4159947262e-01
+6.7400000000e+01,-1.3857134090e+00,-2.0445575724e-01
+6.7500000000e+01,-1.3999016545e+00,-6.5131045442e-02
+6.7600000000e+01,-1.4001103270e+00,7.4984120029e-02
+6.7700000000e+01,-1.3863233547e+00,2.1448990368e-01
+6.7800000000e+01,-1.3586645536e+00,3.5199116110e-01
+6.7900000000e+01,-1.3173965300e+00,4.8611138984e-01
+6.8000000000e+01,-1.2629181940e+00,6.1550650793e-01
+6.8100000000e+01,-1.1957609085e+00,7.3887832211e-01
+6.8200000000e+01,-1.1165833098e+00,8.5498755082e-01
+6.8300000000e+01,-1.0261648514e+00,9.6266627125e-01
+6.8400000000e+01,-9.2539813283e-01,1.0608296656e+00
+6.8500000000e+01,-8.1528009064e-01,1.1484869485e+00
+6.8600000000e+01,-6.9690213679e-01,1.2247513648e+00
+6.8700000000e+01,-5.7143934334e-01,1.2888491584e+00
+6.8800000000e+01,-4.4013877978e-01,1.3401274195e+00
+6.8900000000e+01,-3.0430711905e-01,1.3780607325e+00
+6.9000000000e+01,-1.6529763531e-01,1.4022565556e+00
+6.9100000000e+01,-2.4496722934e-02,1.4124592775e+00
+6.9200000000e+01,1.1668992911e-01,1.4085529099e+00
+6.9300000000e+01,2.5685136959e-01,1.3905623864e+00
+6.9400000000e+01,3.9458548395e-01,1.3586534543e+00
+6.9500000000e+01,5.2851301907e-01,1.3131311586e+00
+6.9600000000e+01,6.5729139346e-01,1.2544369320e+00
+6.9700000000e+01,7.7962815474e-01,1.1831443189e+00
+6.9800000000e+01,8.9429394948e-01,1.0999533741e+00
+6.9900000000e+01,1.0001348753e+00,1.0056837914e+00
+7.0000000000e+01,1.0960840908e+00,9.0126682948e-01
+7.0100000000e+01,1.1811725676e+00,7.8773611460e-01
+7.0200000000e+01,1.2545388742e+00,6.6621741035e-01
+7.0300000000e+01,1.3154378953e+00,5.3791745655e-01
+7.0400000000e+01,1.3632483960e+00,4.0411198738e-01
+7.0500000000e+01,1.3974793563e+00,2.6613304784e-01
+7.0600000000e+01,1.4177750091e+00,1.2535573463e-01
+7.0700000000e+01,1.4239185325e+00,-1.6815506590e-02
+7.0800000000e+01,1.4158343558e+00,-1.5896089103e-01
+7.0900000000e+01,1.3935890564e+00,-2.9965947194e-01
+7.1000000000e+01,1.3573908363e+00,-4.3750334361e-01
+7.1100000000e+01,1.3075875828e+00,-5.7111172816e-01
+7.1200000000e+01,1.2446635305e+00,-6.9914480531e-01
+7.1300000000e+01,1.1692345574e+00,-8.2031714653e-01
+7.1400000000e+01,1.0820421597e+00,-9.3341061892e-01
+7.1500000000e+01,9.8394616594e-01,-1.0372866293e+00
+7.1600000000e+01,8.7591625960e-01,-1.1308975854e+00
+7.1700000000e+01,7.5902239635e-01,-1.2132974590e+00
+7.1800000000e+01,6.3442420858e-01,-1.2836513442e+00
+7.1900000000e+01,5.0335950320e-01,-1.3412439144e+00
+7.2000000000e+01,3.6713196641e-01,-1.3854866930e+00
+7.2100000000e+01,2.2709819765e-01,-1.4159240622e+00
+7.2200000000e+01,8.4654201605e-02,-1.4322379508e+00
+7.2300000000e+01,-5.8778527159e-02,-1.4342511512e+00
+7.2400000000e+01,-2.0176718833e-01,-1.4219292318e+00
+7.2500000000e+01,-3.4288198716e-01,-1.3953810242e+00
+7.2600000000e+01,-4.8071043005e-01,-1.3548576781e+00
+7.2700000000e+01,-6.1387146154e-01,-1.3007502936e+00
+7.2800000000e+01,-7.4102930106e-01,-1.2335861515e+00
+7.2900000000e+01,-8.6090684111e-01,-1.1540235788e+00
+7.3000000000e+01,-9.7229847218e-01,-1.0628454990e+00
+7.3100000000e+01,-1.0740822059e+00,-9.6095173006e-01
+7.3200000000e+01,-1.1652309743e+00,-8.4935010728e-01
+7.3300000000e+01,-1.2448229918e+00,-7.2914651618e-01
+7.3400000000e+01,-1.3120510746e+00,-6.0153393612e-01
+7.3500000000e+01,-1.3662308245e+00,-4.6778060211e-01
+7.3600000000e+01,-1.4068075929e+00,-3.2921740246e-01
+7.3700000000e+01,-1.4333621552e+00,-1.8722463760e-01
+7.3800000000e+01,-1.4456150368e+00,-4.3218271684e-02
+7.3900000000e+01,-1.4434294467e+00,1.0136418609e-01
+7.4000000000e+01,-1.4268127876e+00,2.4507802961e-01
+7.4100000000e+01,-1.3959167263e+00,3.8648579026e-01
+7.4200000000e+01,-1.3510358217e+00,5.2417161355e-01
+7.4300000000e+01,-1.2926047234e+00,6.5675543385e-01
+7.4400000000e+01,-1.2211939671e+00,7.8290680579e-01
+7.4500000000e+01,-1.1375044070e+00,9.0135825331e-01
+7.4600000000e+01,-1.0423603406e+00,1.0109180027e+00
+7.4700000000e+01,-9.3670139225e-01,1.1104819715e+00
+7.4800000000e+01,-8.2157323557e-01,1.1990448930e+00
+7.4900000000e+01,-6.9811724701e-01,1.2757104642e+00
+7.5000000000e+01,-5.6755919198e-01,1.3397004141e+00
+7.5100000000e+01,-4.3119705626e-01,1.3903624017e+00
+7.5200000000e+01,-2.9038814321e-01,1.4271766632e+00
+7.5300000000e+01,-1.4653556523e-01,1.4497613401e+00
+7.5400000000e+01,-1.0742637060e-03,1.4578764342e+00
+7.5500000000e+01,1.4454330321e-01,1.4514263489e+00
+7.5600000000e+01,2.8886166216e-01,1.4304609893e+00
+7.5700000000e+01,4.3043687088e-01,1.3951754084e+00
+7.5800000000e+01,5.6785096374e-01,1.3459080031e+00
+7.5900000000e+01,6.9972615204e-01,1.2831372747e+00
+7.6000000000e+01,8.2473863696e-01,1.2074771856e+00
+7.6100000000e+01,9.4163189660e-01,1.1196711578e+00
+7.6200000000e+01,1.0492293137e+00,1.0205847711e+00
+7.6300000000e+01,1.1464460177e+00,9.1119723239e-01
+7.6400000000e+01,1.2322998204e+00,7.9259170103e-01
+7.6500000000e+01,1.3059211379e+00,6.6594456554e-01
+7.6600000000e+01,1.3665617949e+00,5.3251377718e-01
+7.6700000000e+01,1.4136026258e+00,3.9362635684e-01
+7.6800000000e+01,1.4465597927e+00,2.5066519888e-01
+7.6900000000e+01,1.4650897575e+00,1.0505530329e-01
+7.7000000000e+01,1.4689928563e+00,-4.1750426875e-02
+7.7100000000e+01,1.4582154393e+00,-1.8828568112e-01
+7.7200000000e+01,1.4328505532e+00,-3.3308538574e-01
+7.7300000000e+01,1.3931371577e+00,-4.7470035027e-01
+7.7400000000e+01,1.3394578828e+00,-6.1171176988e-01
+7.7500000000e+01,1.2723353479e+00,-7.4274543847e-01
+7.7600000000e+01,1.1924270774e+00,-8.6648553022e-01
+7.7700000000e+01,1.1005190628e+00,-9.8168781131e-01
+7.7800000000e+01,9.9751803454e-01,-1.0871921492e+00
+7.7900000000e+01,8.8444251871e-01,-1.1819341939e+00
+7.8000000000e+01,7.6241276745e-01,-1.2649561129e+00
+7.8100000000e+01,6.3263966248e-01,-1.3354162729e+00
+7.8200000000e+01,4.9641270142e-01,-1.3925977692e+00
+7.8300000000e+01,3.5508718634e-01,-1.4359157165e+00
+7.8400000000e+01,2.1007074166e-01,-1.4649232281e+00
+7.8500000000e+01,6.2809295497e-02,-1.4793160205e+00
+7.8600000000e+01,-8.5227335854e-02,-1.4789355987e+00
+7.8700000000e+01,-2.3256011667e-01,-1.4637709855e+00
+7.8800000000e+01,-3.7771556739e-01,-1.4339589794e+00
+7.8900000000e+01,-5.1924049945e-01,-1.3897829346e+00
+7.9000000000e+01,-6.5571656203e-01,-1.3316700755e+00
+7.9100000000e+01,-7.8577445506e-01,-1.2601873696e+00
+7.9200000000e+01,-9.0810766621e-01,-1.1760360002e+00
+7.9300000000e+01,-1.0214855939e+00,-1.0800444911e+00
+7.9400000000e+01,-1.1247659246e+00,-9.7316055185e-01
+7.9500000000e+01,-1.2169061401e+00,-8.5644172391e-01
+7.9600000000e+01,-1.2969740387e+00,-7.3104491856e-01
+7.9700000000e+01,-1.3641571638e+00,-5.9821495180e-01
+7.9800000000e+01,-1.4177710455e+00,-4.5927218874e-01
+7.9900000000e+01,-1.4572661706e+00,-3.1559942052e-01
+8.0000000000e+01,-1.4822336108e+00,-1.6862810398e-01
+8.0100000000e+01,-1.4924092513e+00,-1.9824101194e-02
+8.0200000000e+01,-1.4876765751e+00,1.2932693926e-01
+8.0300000000e+01,-1.4680679764e+00,2.7733441340e-01
+8.0400000000e+01,-1.4337645849e+00,4.2271765735e-01
+8.0500000000e+01,-1.3850946046e+00,5.6402075838e-01
+8.0600000000e+01,-1.3225301797e+00,6.9982713331e-01
+8.0700000000e+01,-1.2466828193e+00,8.2877372837e-01
+8.0800000000e+01,-1.1582974235e+00,9.4956469825e-01
+8.0900000000e+01,-1.0582449712e+00,1.0609844268e+00
+8.1000000000e+01,-9.4751393968e-01,1.1619097592e+00
+8.1100000000e+01,-8.2720054051e-01,1.2513213212e+00
+8.1200000000e+01,-6.9849786994e-01,1.3283138136e+00
+8.1300000000e+01,-5.6268407905e-01,1.3921051754e+00
+8.1400000000e+01,-4.2110968221e-01,1.4420445253e+00
+8.1500000000e+01,-2.7518412936e-01,1.4776187986e+00
+8.1600000000e+01,-1.2636177577e-01,1.4984580144e+00
+8.1700000000e+01,2.3872611131e-02,1.5043391172e+00
+8.1800000000e+01,1.7401866344e-01,1.4951883529e+00
+8.1900000000e+01,3.2257539523e-01,1.4710821564e+00
+8.2000000000e+01,4.6805620614e-01,1.4322465371e+00
+8.2100000000e+01,6.0900375625e-01,1.3790549694e+00
+8.2200000000e+01,7.4400456323e-01,1.3120248073e+00
+8.2300000000e+01,8.7170317576e-01,1.2318122561e+00
+8.2400000000e+01,9.9081578104e-01,1.1392059523e+00
+8.2500000000e+01,1.1001431098e+00,1.0351192125e+00
+8.2600000000e+01,1.1985825089e+00,9.2058102802e-01
+8.2700000000e+01,1.2851390606e+00,7.9672589523e-01
+8.2800000000e+01,1.3589356359e+00,6.6478257990e-01
+8.2900000000e+01,1.4192217798e+00,5.2606192864e-01
+8.3000000000e+01,1.4653813406e+00,3.8194384745e-01
+8.3100000000e+01,1.4969387626e+00,2.3386357692e-01
+8.3200000000e+01,1.5135639814e+00,8.3297400558e-02
+8.3300000000e+01,1.5150758681e+00,-6.8252071440e-02
+8.3400000000e+01,1.5014441894e+00,-2.1927089632e-01
+8.3500000000e+01,1.4727900607e+00,-3.6824892119e-01
+8.3600000000e+01,1.4293848864e+00,-5.1369488266e-01
+8.3700000000e+01,1.3716477980e+00,-6.5415133281e-01
+8.3800000000e+01,1.3001416115e+00,-7.8820924202e-01
+8.3900000000e+01,1.2155673453e+00,-9.1452213242e-01
+8.4000000000e+01,1.1187573510e+00,-1.0318196002e+00
+8.4100000000e+01,1.0106671245e+00,-1.1389200910e+00
+8.4200000000e+01,8.9236587719e-01,-1.2347428000e+00
+8.4300000000e+01,7.6502596255e-01,-1.3183185764e+00
+8.4400000000e+01,6.2991125948e-01,-1.3887997234e+00
+8.4500000000e+01,4.8836463013e-01,-1.4454685933e+00
+8.4600000000e+01,3.4179457524e-01,-1.4877448914e+00
+8.4700000000e+01,1.9166122011e-01,-1.5151916141e+00
+8.4800000000e+01,3.9461770468e-02,-1.5275195609e+00
+8.4900000000e+01,-1.1328441687e-01,-1.5245903730e+00
+8.5000000000e+01,-2.6505099974e-01,-1.5064180672e+00
+8.5100000000e+01,-4.1431990131e-01,-1.4731690485e+00
+8.5200000000e+01,-5.5959649343e-01,-1.4251605990e+00
+8.5300000000e+01,-6.9942456078e-01,-1.3628578570e+00
+8.5400000000e+01,-8.3240089599e-01,-1.2868693156e+00
+8.5500000000e+01,-9.5718937950e-01,-1.1979408833e+00
+8.5600000000e+01,-1.0725344028e+00,-1.0969485654e+00
+8.5700000000e+01,-1.1772735004e+00,-9.8488983707e-01
+8.5800000000e+01,-1.2703490635e+00,-8.6287379398e-01
+8.5900000000e+01,-1.3508190164e+00,-7.3211017729e-01
+8.6000000000e+01,-1.4178663502e+00,-5.9389738137e-01
+8.6100000000e+01,-1.4708074143e+00,-4.4960956342e-01
+8.6200000000e+01,-1.5090988839e+00,-3.0068298287e-01
+8.6300000000e+01,-1.5323433316e+00,-1.4860170627e-01
+8.6400000000e+01,-1.5402933463e+00,5.1171799463e-03
+8.6500000000e+01,-1.5328541565e+00,1.5893870163e-01
+8.6600000000e+01,-1.5100847300e+00,3.1132532299e-01
+8.6700000000e+01,-1.4721973384e+00,4.6075231363e-01
+8.6800000000e+01,-1.4195555890e+00,6.0572300297e-01
+8.6900000000e+01,-1.3526709408e+00,7.4478376935e-01
+8.7000000000e+01,-1.2721977403e+00,8.7653861393e-01
+8.7100000000e+01,-1.1789268229e+00,9.9966317307e-01
+8.7200000000e+01,-1.0737777448e+00,1.1129180287e+00
+8.7300000000e+01,-9.5778971894e-01,1.2151611826e+00
+8.7400000000e+01,-8.3211134682e-01,1.3053595702e+00
+8.7500000000e+01,-6.9798924625e-01,1.3825994954e+00
+8.7600000000e+01,-5.5675568784e-01,1.4460958843e+00
+8.7700000000e+01,-4.0981536295e-01,1.4952002615e+00
+8.7800000000e+01,-2.5863141455e-01,1.5294073687e+00
+8.7900000000e+01,-1.0471086974e-01,1.5483603591e+00
+8.8000000000e+01,5.0410380994e-02,1.5518545131e+00
+8.8100000000e+01,2.0518290800e-01,1.5398394362e+00
+8.8200000000e+01,3.5805921745e-01,1.5124197181e+00
+8.8300000000e+01,5.0750922237e-01,1.4698540414e+00
+8.8400000000e+01,6.5203555523e-01,1.4125527507e+00
+8.8500000000e+01,7.9018856870e-01,1.3410739017e+00
+8.8600000000e+01,9.2058087447e-01,1.2561178310e+00
+8.8700000000e+01,1.0419012742e+00,1.1585202974e+00
+8.8800000000e+01,1.1529279427e+00,1.0492442632e+00
+8.8900000000e+01,1.2525407313e+00,9.2937039521e-01
+8.9000000000e+01,1.3397324666e+00,8.0008638010e-01
+8.9100000000e+01,1.4136191314e+00,6.6267515889e-01
+8.9200000000e+01,1.4734488252e+00,5.1850219768e-01
+8.9300000000e+01,1.5186094120e+00,3.6900192074e-01
+8.9400000000e+01,1.5486347796e+00,2.1566344058e-01
+8.9500000000e+01,1.5632096447e+00,6.0015727143e-02
+8.9600000000e+01,1.5621728561e+00,-9.6387636317e-02
+8.9700000000e+01,1.5455191589e+00,-2.5198395858e-01
+8.9800000000e+01,1.5133994039e+00,-4.0521705242e-01
+8.9900000000e+01,1.4661191950e+00,-5.5455279715e-01
+9.0000000000e+01,1.4041359898e+00,-6.9849449614e-01
+9.0100000000e+01,1.3280546783e+00,-8.3559787564e-01
+9.0200000000e+01,1.2386216840e+00,-9.6448557451e-01
+9.0300000000e+01,1.1367176438e+00,-1.0838609797e+00
+9.0400000000e+01,1.0233487393e+00,-1.1925212680e+00
+9.0500000000e+01,8.9963676451e-01,-1.2893695242e+00
+9.0600000000e+01,7.6680802755e-01,-1.3734258115e+00
+9.0700000000e+01,6.2618119690e-01,-1.4438370844e+00
+9.0800000000e+01,4.7915421199e-01,-1.4998858430e+00
+9.0900000000e+01,3.2719038837e-01,-1.5409974409e+00
+9.1000000000e+01,1.7180385543e-01,-1.5667459724e+00
+9.1100000000e+01,1.4544471623e-02,-1.5768586788e+00
+9.1200000000e+01,-1.4301763253e-01,-1.5712188280e+00
+9.1300000000e+01,-2.9930772962e-01,-1.5498670378e+00
+9.1400000000e+01,-4.5276223211e-01,-1.5130010269e+00
+9.1500000000e+01,-6.0184433347e-01,-1.4609737961e+00
+9.1600000000e+01,-7.4505939728e-01,-1.3942902546e+00
+9.1700000000e+01,-8.8096994022e-01,-1.3136023248e+00
+9.1800000000e+01,-1.0082100586e+00,-1.2197025728e+00
+9.1900000000e+01,-1.1254991541e+00,-1.1135164263e+00
+9.2000000000e+01,-1.2316548199e+00,-9.9609305661e-01
+9.2100000000e+01,-1.3256047591e+00,-8.6859501499e-01
+9.2200000000e+01,-1.4063976131e+00,-7.3228672455e-01
+9.2300000000e+01,-1.4732125935e+00,-5.8852194290e-01
+9.2400000000e+01,-1.5253678168e+00,-4.3873031897e-01
+9.2500000000e+01,-1.5623272599e+00,-2.8440317780e-01
+9.2600000000e+01,-1.5837062643e+00,-1.2707867460e-01
+9.2700000000e+01,-1.5892755328e+00,3.1673534183e-02
+9.2800000000e+01,-1.5789635766e+00,1.9026795005e-01
+9.2900000000e+01,-1.5528575884e+00,3.4711906562e-01
+9.3000000000e+01,-1.5112027289e+00,5.0065721376e-01
+9.3100000000e+01,-1.4543998347e+00,6.4934427422e-01
+9.3200000000e+01,-1.3830015669e+00,7.9168908061e-01
+9.3300000000e+01,-1.2977070386e+00,9.2626237348e-01
+9.3400000000e+01,-1.1993549726e+00,1.0517111495e+00
+9.3500000000e+01,-1.0889154575e+00,1.1667722626e+00
+9.3600000000e+01,-9.6748038144e-01,1.2702851407e+00
+9.3700000000e+01,-8.3625263937e-01,1.3612034898e+00
+9.3800000000e+01,-6.9653421959e-01,1.4386058670e+00
+9.3900000000e+01,-5.4971328793e-01,1.5017050168e+00
+9.4000000000e+01,-3.9725039750e-01,1.5498558745e+00
+9.4100000000e+01,-2.4066396106e-01,1.5825621572e+00
+9.4200000000e+01,-8.1515130442e-02,1.5994814730e+00
+9.4300000000e+01,7.8607766308e-02,1.6004288982e+00
+9.4400000000e+01,2.3810507510e-01,1.5853789831e+00
+9.4500000000e+01,3.9538179501e-01,1.5544661657e+00
+9.4600000000e+01,5.4886352695e-01,1.5079835884e+00
+9.4700000000e+01,6.9701223239e-01,1.4463803258e+00
+9.4800000000e+01,8.3834164430e-01,1.3702570521e+00
+9.4900000000e+01,9.7143217613e-01,1.2803601880e+00
+9.5000000000e+01,1.0949451790e+00,1.1775745854e+00
+9.5100000000e+01,1.2076364041e+00,1.0629148204e+00
+9.5200000000e+01,1.3083685351e+00,9.3751518111e-01
+9.5300000000e+01,1.3961226639e+00,8.0261844912e-01
+9.5400000000e+01,1.4700085944e+00,6.5956358458e-01
+9.5500000000e+01,1.5292738699e+00,5.0977243794e-01
+9.5600000000e+01,1.5733114326e+00,3.5473561973e-01
+9.5700000000e+01,1.6016658374e+00,1.9599766878e-01
+9.5800000000e+01,1.6140379564e+00,3.5141666621e-02
+9.5900000000e+01,1.6102881253e+00,-1.2622654884e-01
+9.6000000000e+01,1.5904376991e+00,-2.8649441295e-01
+9.6100000000e+01,1.5546689998e+00,-4.4405874714e-01
+9.6200000000e+01,1.5033236550e+00,-5.9734179408e-01
+9.6300000000e+01,1.4368993417e+00,-7.4480701484e-01
+9.6400000000e+01,1.3560449685e+00,-8.8497448968e-01
+9.6500000000e+01,1.2615543392e+00,-1.0164357684e+00
+9.6600000000e+01,1.1543583635e+00,-1.1378680209e+00
+9.6700000000e+01,1.0355158882e+00,-1.2480473465e+00
+9.6800000000e+01,9.0620324065e-01,-1.3458611071e+00
+9.6900000000e+01,7.6770258761e-01,-1.4303191607e+00
+9.7000000000e+01,6.2138922417e-01,-1.5005638816e+00
+9.7100000000e+01,4.6871791914e-01,-1.5558788664e+00
+9.7200000000e+01,3.1120845234e-01,-1.5956962361e+00
+9.7300000000e+01,1.5043048750e-01,-1.6196024617e+00
+9.7400000000e+01,-1.2012068502e-02,-1.6273426523e+00
+9.7500000000e+01,-1.7449706454e-01,-1.6188232618e+00
+9.7600000000e+01,-3.3540030233e-01,-1.5941131851e+00
+9.7700000000e+01,-4.9311177034e-01,-1.5534432324e+00
+9.7800000000e+01,-6.4605175200e-01,-1.4972039837e+00
+9.7900000000e+01,-7.9268664731e-01,-1.4259420443e+00
+9.8000000000e+01,-9.3154434919e-01,-1.3403547370e+00
+9.8100000000e+01,-1.0612290206e+00,-1.2412832821e+00
+9.8200000000e+01,-1.1804351238e+00,-1.1297045324e+00
+9.8300000000e+01,-1.2879605616e+00,-1.0067213442e+00
+9.8400000000e+01,-1.3827187968e+00,-8.7355167876e-01
+9.8500000000e+01,-1.4637498293e+00,-7.3151654221e-01
+9.8600000000e+01,-1.5302299193e+00,-5.8202688397e-01
+9.8700000000e+01,-1.5814799576e+00,-4.2656958223e-01
+9.8800000000e+01,-1.6169723992e+00,-2.6669265680e-01
+9.8900000000e+01,-1.6363366880e+00,-1.0398985586e-01
+9.9000000000e+01,-1.6393631183e+00,5.9915229766e-02
+9.9100000000e+01,-1.6260050919e+00,2.2338536661e-01
+9.9200000000e+01,-1.5963797472e+00,3.8478603108e-01
+9.9300000000e+01,-1.5507669520e+00,5.4250175129e-01
+9.9400000000e+01,-1.4896066684e+00,6.9495227479e-01
+9.9500000000e+01,-1.4134947155e+00,8.4060840038e-01
+9.9600000000e+01,-1.3231769691e+00,9.7800731548e-01
+9.9700000000e+01,-1.2195420558e+00,1.1057672852e+00
+9.9800000000e+01,-1.1036126131e+00,1.2226015456e+00
+9.9900000000e+01,-9.7653520048e-01,1.3273312614e+00
+1.0000000000e+02,-8.3956896276e-01,1.4188974183e+00
diff --git a/doc/remarques.typ b/doc/remarques.typ
new file mode 100644
index 0000000..28a8858
--- /dev/null
+++ b/doc/remarques.typ
@@ -0,0 +1,228 @@
+La phrase "it is less computationally intensive... the compiler". est mal dite;
+tu peux dire que c'est la methode classique, par exemple, implementee pour
+simulink et aussi Zelus.
+
+Rmq:
+
+- Les methodes numeriques (centralisees) les premieres ont ete faites
+ pour resoudre un probleme initial. dx/dt = f (x(t)) t et x(0) = x0 ou x est un
+ vecteur.
+
+- La "simulation distribuee" (cad faisant intervenir plusieurs solveurs)
+ existe aussi mais pose d'autres problemes: efficacite, convergence. Elle est
+ plus difficile a implementer.
+
+- Sundials CVODE a une methode centralisee mais avec une implementation de
+ manipulation des vecteurs qui peut etre parallele, cad que tous les calculs de
+ vecteur/matrice peuvent se faire en parallele. On ne l'utilise pas dans Zelus.
+
+p2.
+"we cannot necessarily use "pre"...": a mon avis, le lecteur ne peut pas
+comprendre ce que tu ecris ici. A mon avis, tu peux dire ce qui marche (si je ne
+me trompe pas, pouvoir internaliser le solveur au modele et donc le composer
+avec un observateur qui contient son propre solveur) et ce qui est a faire (si
+je comprends, pouvoir ecrire directement assert dans le code et la compilation
+pour pouvoir co-simuler le modele et ses differents observateurs/assertions.
+Pour le "pre", c'est un point que tu peux mentionner plus loin dans le rapport
+lorsque le lecteur a plus de contexte et de matiere sur ce que tu as fait.
+
+Il n'y a pas de numero en bas des pages dans ton texte. Il en faut. Impossible
+de savoir ou on est sinon.
+
+section 2. Tu parles de streams (Nat -> V) alors que tu donnes une semantique
+avec des etats. C'est bizarre. Tu ecris pre(s) alors que l'on ne sait pas ce que
+c'est. le
+
+```ocaml
+let rec (o, s) = M.f_step(i, M.s0 -> pre(s))
+```
+
+n'est pas du caml valide.
+
+Que veux-tu dire ? Quel est le type de Simulate(M) ? Je pensais (cf. page
+d'apres) que tu voudrais ecrire:
+
+```ocaml
+simulate(M): list(I) -> list(O)
+```
+
+Le listing 1 ne correspond pas a la page precedente (ou l'inverse). Pourquoi ?
+
+Page d'apres (5 ?).
+
+Je ne suis pas fan de l'exemple sincos_discrete() tel que tu l'as ecrit car il
+diverge. Il laisse soupconner que l'on ne peut pas ecrire de schemas numeriques
+avec des streams. C'est faux, en general.
+
+En somme, le code listing 2 est un mauvais programme. Il y a une raison
+mathematique liee a l'analyse des poles; je te montrerai a mon retour. Si tu
+veux garder cet exemple de sin/cos, il faut mettre deux integrateurs differents,
+un forward et un backward.
+
+Ou alors, il faut faire un discours (plus long) pour expliquer pourquoi c'est
+mieux d'ecrire la forme a temps continu, que justement, elle permet de ne pas se
+casser la tete pour savoir ou et quand mettre un pre. Mais c'est un peu
+orthogonal a ce que tu veux raconter.
+
+Discrete models -> Discrete-time models.
+Continuous models -> Continuous-time models.
+
+Tu pourras ensuite dire que tu confonds les deux.
+
+2.2 la raison de la divergence n'est pas la. Si tu prends un schema backward +
+forward, ca ne diverge pas. Par contre, si tu prends Van der Pol en discret, pas
+fixe, avec un pas trop petit ou trop petit, tu tombes sur nan. Regarde. J'avais
+fait l'exemple (regarde sur le depot; ou refais le).
+
+Rmq: sin/cos n'est pas un super exemple pour justifier l'interet de Zelus (ou du
+temps continu en general) car ce systeme est lineaire. Un schema numerique pas
+fixe (donc programmable directement en synchrone) marche bien. C'est plutot
+quand la dynamique est compliquee, non lineaire, qu'il te faut un schema
+numerique plus elabore et que l'expression directe d'un modele a temps continu
+se justifie. La modelisation d'evenement en etat (zero-crossing) le justifie
+aussi. Le cas classique, c'est le modele bang-bang. Tu veux detecter les
+instants ou la temperature passe un seuil. Avec un modele de temps pas fixe, tu
+ne peux pas.
+
+avant de donner CNode(...), dit ce qu'est le probleme.
+Qu'est ce qu'un modele M ?
+
+```ocaml
+model (i) = o where
+ der s = f_der i s
+ der o = f_out i s
+```
+
+2.3. Tu n'as pas dit/defini ce qu'etait un "initial value problem". Le faire
+(sinon, on ne comprends pas).
+
+"the integr node from Listing 2 is another example of a numerical solver
+(albeit not a very good one)"
+
+A mon avis, inutile. Que veux-tu dire que tu n'as pas deja dit (on peut calculer
+une approximation de `(der(x)/dt)(t)` a des instants `n * h` (`n in Nat`) par
+`(x - pre(x))/h`?
+
+
+"Of particular interest is the fact that numerical ODE solvers compute
+approximations sequentially. "
+
+Ca depend. Les solveurs a memoire (dits "multi-pas"). Les solveurs sans memoire
+tels que RK, non. A mon avis, tu peux commencer par exempliquer ce que calcule
+un solveur: etant donne le IVP, l'objectif est de calculer une solution
+[0, tmax] -> V. En pratique, il ne calcule pas toute la solution. Il calcule une
+suite de solution pour [0, h0][h0, h1][h1, h2], etc. Relis l'article de
+Chapoutot et al., de memoire, c'est tres bien explique. Sinon, il y a les livres
+jaunes dans l'etagere de mon bureau.
+
+"This sequential process gives way to a synchronous interpretation of an ODE
+solver as a discrete node."
+
+Qu'est-ce qu'un "discrete node" ?
+
+IVP(Y, Y'). Aurait du arriver bien plus tot.
+
+dy /dt (t) = f t (y(t))
+y(0) = y0
+
+donner l'horizon de temps (par defaut, +infty ?).
+
+Je prendrais une autre lettre que "s". (on note s pour "state" en general). "o",
+pour "output" ?
+
+La notation pointee est un peu troublante (en math, on la confond avec la
+multiplication. Utiliser une police ad-hoc pour h et u de sorte que s_0.u ne
+soit pas ambigu.
+
+mets des numeros de ligne dans ton texte stp. Il y a un style latex pour ca. Ou
+au moins des no. de page.
+
+"The ODE solver does not itself introduce discontinuities; the only
+discontinuities in the system are those introduced by the input signal."
+
+Avant de dire cela, tu ne parles par d'entrees. Pour le moment, ton systeme
+n'est pas "hybride". C'est une ODE. Tu ne donnes pas d'hypotheses sur f. En
+faut-il ? Lesquelles ?
+
+"The simulation of a continuous-time system with an ODE solver is now itself a
+synchronous node:"
+
+c'est un point de vue. Puisque la simulation d'un modele a temps continu calcule
+une suite de solutions approchees, ne pourrait-on pas la voir elle meme comme
+une machine synchrone qui...
+
+page suivante.
+
+Plutot que Signal(V), j'utiliserais un autre nom. Tu as deja un langage qui
+manipule des signaux, non ?
+
+Superdense(V) ? En reference aux travaux de Pnueli et al. Puis Edward Lee et al.
+Puis Caspi et al.
+
+Superdence(V) = R+ -> Nat -> V avec, a cote, une fonction qui pour chaque t in
+R+, donne le nombre de sauts (relire les articles; j'ai oublie le nom de cette
+fonction). step(t) in Nat.
+
+est-ce tres different de
+
+Nat -> R+ -> V avec a cote, une fonction qui indique le temps qui s'ecoule pour
+chaque n in Nat ?
+
+Dans un des articles de Edward Lee, il explique le type des signaux pour le
+temps continu. Regarde son livre (vert) sur mon etagere.
+
+Le texte de cette page (8 ?) est un peu confus.
+
+est-ce que les bouts de fonctions sont collees les uns aux autres ?
+
+none a le sens de "not yet" ? Comment recolles-tu les morceaux ?
+
+(h0, u0) (h1, u1) ... (hn, un)
+
+est-ce equivalent a :
+
+(h0, u0) none (h1, u1) none none (hn, un) ?
+
+A ce stade, tu n'as pas dit ce qu'etait $italic("CNode")(I, O, S_M, S'_M)$. A
+quoi ca ressemble ?
+
+Les elements du listing 4 doivent etre donnes dans le texte en maths. Le code
+caml doit etre mis aussi dans le listing 4. Je comprends bien mieux le code que
+les explications ! Sois plus precis dans le texte, ca aidera a comprendre.
+
+utiliser hmax plutot que h pour ivp. Eviter les surcharges de noms autant que
+possible.
+
+Plutot que "discrete node", utiliser "synchronous machine" ? Le terme "machine"
+est un terme technique des annees 60/70. Cf. articles historiques sur les
+automates et les "sequential machines". Tu peux dire, qui se caraterise par un
+etat initial, une fonction step lisant une entree et produisant une sortie,
+(plus une fonction reset) si tu veux.
+
+"Since time is logical in discrete nodes". Un peu trop pedant (et obscur) a mon
+avis. On manipule des suites de valeurs; il n'y a donc pas de temps reel ou
+physique. Le temps est juste la succession des valeurs.
+
+Le paragraphe "The question of discrete events..." est un peu confus. Donne un
+exemple. E.g., un timer: un controleur ouvre un robinet pendant 4s; une lumiere
+doit clignoter en alternant allume/eteint 1s/2s. ce sont des cas ou il y a des
+evenements en temps, cad en reference a la variable t (der t = 1 init 0).
+D'autres evenements ne sont pas des evenements en temps. E.g., une roue qui
+tourne (e.g. volant moteur). Compter le nombre de tours de roue par seconde avec
+un capteur. Le schema general est dit de "Zero-crossing". Etc. La balle qui
+rebondit, etc. Tu en trouveras d'autres !
+
+
+La suite est pas mal. Il faut expliquer (et possiblement donner le code, au
+moins en annexe) pour 2.5, 2.6, 2.7.
+
+Ca prend forme. Tu peux donner l'intuition de la compilation. Dans le source, on
+melange les der, present, up, etc. et le compilo., apres une analyse statique
+(qui va verifier des proprietes et donc rejeter certains programmes), va
+produire les differentes fonctions necessaires a la simulation qui alterne des
+phases d'integration (le temps ronronne) et des pas discrets (reactions
+instantanees).
+
+Super. Continue ! --Marc
+
+
diff --git a/doc/rep.typ b/doc/rep.typ
index 27f2899..38ca27c 100644
--- a/doc/rep.typ
+++ b/doc/rep.typ
@@ -1,24 +1,65 @@
#import "@preview/cetz:0.4.0"
+#import "@preview/cetz-plot:0.1.2"
#import "@preview/finite:0.4.1"
-#set page(paper: "a4")
-#set par(justify: true)
-#set raw(syntaxes: "zelus.sublime-syntax")
-#show raw: set text(font: "CaskaydiaCove NF")
+#set text(size: 11pt)
+#set page(paper: "a4", numbering: "1")
+#set par(
+ justify: true,
+ // first-line-indent: 10pt
+)
+#set raw(
+ syntaxes: "zelus.sublime-syntax",
+ // theme: none,
+)
+// #show raw: set text(font: "CaskaydiaCove NF")
+#set cite(style: "alphanumeric.csl")
+#set heading(numbering: "1.1")
-#let lustre = smallcaps[Lustre]
-#let esterel = smallcaps[Esterel]
-#let simulink = smallcaps[Simulink]
+#set figure(gap: 20pt, placement: top)
+
+#let zelus = smallcaps[Zélus]
+#let lustre = smallcaps[Lustre]
+#let ocaml = smallcaps[OCaml]
#let sundials = smallcaps[Sundials CVODE]
-#let zelus = smallcaps[Zélus]
-#let note(color, prefix, ..what) = {
- let msg = if what.pos().len() == 0 { "" } else { ": " + what.pos().join("") }
- block(fill: color, width: 100%, inset: 5pt, align(center, raw(prefix + msg)))
-}
-#let TODO(..what) = note(red, "TODO", ..what)
-#let MENTION(..what) = note(gray, "MENTION", ..what)
-#let adot(s) = $accent(#s, dot)$
-#let addot(s) = $accent(#s, dot.double)$
+#let simulink = smallcaps[Simulink]
+#let haskell = smallcaps[Haskell]
+
+#let zel(body) = raw(lang: "zelus", body)
+
+#let Time = math.italic[Time]
+#let stop = math.italic[stop]
+#let DNode = math.italic[DNode]
+#let DNodeC = math.italic[DNodeC]
+#let CNode = math.italic[CNode]
+#let HNode = math.italic[HNode]
+#let HNodeA = math.italic[HNodeA]
+#let Dense = math.italic[Dense]
+#let IVP = math.italic[IVP]
+#let ZCP = math.italic[ZCP]
+#let Stream = math.italic[Stream]
+#let Signal = math.italic[Signal]
+#let Solver = math.italic[Solver]
+#let CSolver = math.italic[CSolver]
+#let ZSolver = math.italic[ZSolver]
+#let Option = math.italic[Option]
+#let None = math.italic[None]
+#let Some = math.italic[Some]
+#let Unit = math.italic[Unit]
+#let List = math.italic[List]
+#let Zi = $Z_i$
+#let Zo = $Z_o$
+#let SimState = math.italic[State]
+#let DSim = math.italic[DSim]
+#let CSim = math.italic[CSim]
+#let HSim = math.italic[HSim]
+#let HASim = math.italic[HASim]
+
+#let show_notes = true
+#let note(body, prefix: "NOTE: ", color: rgb(0, 0, 0, 50)) = if show_notes {
+ rect(fill: color, inset: 3pt, width: 100%)[#prefix #body]
+} else []
+#let todo(body) = note(color: rgb(255, 0, 0, 50), prefix: "TODO: ")[#body]
#align(center)[
#text(16pt)[
@@ -29,215 +70,1305 @@
Henri Saudubray, supervised by Marc Pouzet, Inria PARKAS
]\
-#heading(outlined: false)[General Context]
-What is the report about? Where does the problem come from? What is the state of
-the art in that area?
+// #heading(outlined: false, numbering: none)[General Context]
+// What is the report about? Where does the problem come from? What is the state
+// of the art in that area?
-Hybrid systems modelers such as #simulink are essential tools in the development
-of embedded systems evolving in physical environments. They allow for precise
-descriptions of hybrid systems through both continuous-time behaviour defined
-through Ordinary Differential Equations (ODEs) and discrete-time reactive
-behaviour similar to what is found in the synchronous languages such as #lustre.
-The #zelus language aims to reconcile synchronous languages with hybrid systems,
-by taking a synchronous language kernel and extending it with continuous-time
-constructs similar to what is found in tools like #simulink. Continuous-time
-behaviour is computed through the use of an off-the-shelf ODE solver such as
-#sundials.
+*General Context* --- Hybrid systems modelers such as #simulink #footnote[
+ #link("https://www.mathworks.com/products/simulink")
+] are essential tools in the development of embedded systems which evolve in
+physical environments. They allow for precise descriptions of hybrid systems
+through both continuous-time behaviour defined with Ordinary Differential
+Equations (ODEs) and discrete-time reactive behaviour similar to what is found
+in the synchronous languages such as #lustre @cit:lustre. The #zelus language
+@cit:zelus_sync_lng_with_ode aims to reconcile synchronous languages with hybrid
+systems, by taking a synchronous language kernel and extending it with
+continuous-time constructs similar to what is found in tools like #simulink.
+Continuous-time behaviour is computed through the use of an off-the-shelf ODE
+solver such as #sundials @cit:sundials @cit:sundialsml.
-#heading(outlined: false)[Research Problem]
-What specific question(s) have you studied? What motivates the question? What
-are the applications/consequences? Is it a new problem? Why did you choose this
-problem?
+// #heading(outlined: false, numbering: none)[Research Problem]
+// What specific question(s) have you studied? What motivates the question? What
+// are the applications/consequences? Is it a new problem? Why did you choose
+// this problem?
-The simulation of hybrid system models, as done in #simulink and #zelus, uses a
-single ODE solver instance to simulate the entire model. This has various
-advantages: it is less computationally intensive, and simplifies the work of the
-compiler. Unfortunately, it also raises a difficult problem: sub-systems which
-seemingly should not interfere with each other end up affecting each other's
-results. This is due to the chosen integration method: an adaptive solver like
-#sundials will vary its step length throughout the integration process, and the
-addition of ODEs in the overall system can influence these step lengths,
-affecting the results obtained for pre-existing ODEs. This is particularly
-problematic in the case of runtime assertions, which are typically expected to
-be transparent: they should not affect the final result of the computation.
+*Research Problem* --- The simulation of hybrid system models, as in #simulink
+and #zelus, uses a single ODE solver instance to simulate the entire model. This
+has various advantages: it is less computationally intensive, and simplifies the
+work of the compiler. Unfortunately, it also raises a difficult problem:
+sub-systems which seemingly should not interfere with each other end up
+affecting each other's results. This is due to the chosen integration method. An
+adaptive solver like #sundials will vary its step length throughout the
+integration process, and the addition of new, unrelated ODEs in the system can
+influence these step lengths, affecting the results obtained for pre-existing
+ODEs. This is particularly problematic in the case of runtime assertions, which
+are typically expected to be transparent: they should not affect the final
+result of the computation.
-#heading(outlined: false)[Proposed Contributions]
-What is your answer to the research problem? Please remain at a high level;
-technical details should be given in the main body of the report. Pay special
-attention to the description of the _scientific_ approach.
+We therefore aim to define a new execution model for hybrid system models, which
+allows for clear separation between a program and its assertions, in such a way
+that the results obtained by executing a model with and without its assertions
+are the same.
-To solve this, we propose a new runtime for the #zelus language, which simulates
-assertions with their own solvers in order to maintain the separation between
-assertions and the model they operate on.
+// #heading(outlined: false, numbering: none)[Proposed Contributions]
+// What is your answer to the research problem? Please remain at a high level;
+// technical details should be given in the main body of the report. Pay special
+// attention to the description of the _scientific_ approach.
-#TODO()
+*Proposed Contributions* --- To solve this, we propose a new runtime for the
+#zelus language that simulates assertions with their own solvers in order to
+maintain the separation between assertions and the model they operate on. This
+is done by first describing the execution of a hybrid system model as a
+synchronous node akin to those found in languages such as #lustre, operating on
+streams of functions, and then implementing this vision in #ocaml. This
+interpretation allows us to _lift_ the runtime into the language, allowing for
+direct manipulation of hybrid simulations in the synchronous subset of #zelus.
+The addition of a few primitives then allows us to attain our final objective of
+transparent assertions.
-#heading(outlined: false)[Arguments Supporting Their Validity]
-What is the evidence that your solution is a good solution? (Experiments?
-Proofs?) Comment on the robustness of your solution: does it rely on
-assumptions, and if so, to what extent?
+// #heading(outlined: false, numbering: none)[Arguments Supporting Their Validity]
+// What is the evidence that your solution is a good solution? (Experiments?
+// Proofs?) Comment on the robustness of your solution: does it rely on
+// assumptions, and if so, to what extent?
-#TODO("Justify")
+*Arguments Supporting Their Validity* ---
+#todo[
+ Justify.
+]
-#heading(outlined: false)[Summary and Future Work]
-What did you contribute to the area? What comes next? What is a good _next_ step
-or question?
+// #heading(outlined: false, numbering: none)[Summary And Future Work]
+// What did you contribute to the area? What comes next? What is a good _next_
+// step or question?
-#TODO("Next steps")
+*Summary and Fugure Work* -- Future work includes the addition of assertions as
+a syntactic construct in the source language, with a compilation pass to
+translate them to their equivalent form in #ocaml. The direct manipulation of
+streams of functions which may depend on the solver's internal state also raises
+several questions, one of which is of particular interest: we cannot necessarily
+use `pre` on such values, since an integration step of the solver may modify its
+internal state, on which the closures depend. One could envision a static
+analysis to forbid such manipulations, and a simple answer could be typeclasses
+_à la_ #haskell, but the addition of typeclasses into the type system of #zelus
+is a problem of its own.
#pagebreak()
#outline()
+
+= Introduction
+
+#todo[
+ Write the introduction.
+ Justify the approach: we first present a low-level denotational semantics
+ of hybrid system models as a collection of functions operating on an inner
+ state, and pair this denotational semantics with a numerical ODE solver and a
+ zero-crossing detection mechanism to construct a synchronous operational
+ semantics of the simulation. We then use this operational semantics to
+ implement transparent observers, by separating the simulation of the model
+ from that of its observers.
+]
+
#pagebreak()
-= Background
+= Hybrid System Model Simulation
-== Hybrid systems
+The simulation of a hybrid system model is a function from signals to signals.
+Signals are functions from time to values, modelling the evolution of a value in
+time. The exact meaning of time depends on the nature of the model. Three
+possible situations may occur: discrete-time models, akin to those found in
+synchronous languages like #lustre @cit:lustre; continuous-time models with no
+discontinuities; and hybrid models, which involve both discrete and continuous
+behaviours.
-Hybrid systems are dynamical systems whose behaviour evolves through both
-continuous and discrete phases. They are used to model software interacting with
-physical systems. Continuous phases are described using ordinary differential
-equations (ODEs), while discrete phases can be represented as a reactive
-program in a synchronous language such as #lustre or #esterel.
+== Discrete-Time Models
-As an illustration, say we wished to model an extensively studied system: a
-bouncing ball. We could start by describing its distance from the ground $y$ as
-a function of time, with a second-order differential equation
-$ addot(y) = -9.81, $
-where $addot(y)$ denotes the second order derivative of $y$ with
-respect to time $(d^2y)/(d t^2)$ (the acceleration of the ball), and $9.81$ is
-the gravitational constant $g$: the acceleration of the ball is its negation. We
-now give the initial position and speed of the ball:
-$ y(0) = 50 space space space adot(y)(0) = 0 $
-We have just described an initial value problem: given an ODE and an initial
-value for its dependent variable, its solution is a function $y(t)$ returning
-the value of the variable at a given time $t$. We can find an approximation of
-this function using an ODE solver, such as #sundials.
+#zelus starts from a synchronous language kernel _à la_ #lustre, and extends it
+with continuous-time constructs @cit:zelus_sync_lng_with_ode. This synchronous
+language kernel allows for the description of variables evolving in time through
+streams of values. The notion of time is logical, and is represented as a series
+of discrete instants; time is then a value in $NN$, and streams are functions
+from $NN$ to their respective codomains:
-As of right now, our ball will fall until the end of time; we have not said
-anything about how it bounces when it hits the floor. To do so, we need a notion
-of _events_: we need to identify exactly when the ball hits the ground, so that
-we may take action to make it bounce. These events are modelled by
-zero-crossings: we monitor a certain value and stop when it goes from strictly
-negative to positive or null. For our purposes, we choose $-y$ as the monitored
-value, and call the zero-crossing event $z$. When $z$ occurs (i.e., when the
-ball touches the ground), we set the speed $adot(y)$ to
-$-k dot #raw(lang: "zelus", "last")\(adot(y))$, where
-$#raw(lang: "zelus", "last")\(y)$ denotes the left limit of $y$ (we cannot
-specify $adot(y)$ in terms of itself), and $k$ is a factor modelling the loss of
-inertia due to the collision (say, $0.8$). We can then resume the approximation
-of the solution.
+$ Stream(V) = NN -> V $
-@lst:ball.zls shows how such a model might be expressed in the concrete syntax
-of #zelus @cit:zelus_sync_lng_with_ode.
+Given a stream $s : Stream(V)$, we denote $s_n$ the $n$-th value $s(n)$ of the
+stream.
-#figure(placement: top, gap: 2em,
+Computation occurs in successive steps performed at each instant, and may depend
+on values computed at previous instants. This interpretation of time allows a
+program written in the synchronous kernel to be considered independently of its
+physical implementation. Nothing in this representation tells us anything about
+how much physical time passes between successive instants.
+
+#note[Should we give a precise definition of the kernel?]
+
+The programs expressed in this kernel, called _discrete nodes_, are functions on
+streams. At each time instant, given inputs $I$, they produce outputs $O$.
+Producing these outputs may also depend on previously computed values through
+operators like ```zelus pre(e)```, which returns the value of its sub-expression
+`e` at the previous instant, and ```zelus e1 -> e2```, which returns its
+left-hand side `e1` at the first instant and its right-hand side `e2`
+afterwards. This requires nodes to store some previously computed values. Nodes
+therefore operate on an inner state of type $S$: previously computed values must
+be stored inside this state in order to refer to them afterwards. The behaviour
+of a node is represented by a step function $f_"step" : S -> I -> O times S$ and
+an initial state $s_0 : S$, used at the first instant. Given a set of inputs and
+the current state, the $f_"step"$ function produces a set of outputs and a new,
+updated state. This function is then called at each instant, taking as input the
+current value of the input signal, and the state produced by the previous
+instant.
+
+#todo[
+ Clearly distinguish between stream semantics and Mealy machine implementation
+ (and mention Mealy machines explicitly).
+]
+
+Since programs may wish to reset the state of a node (for instance, when writing
+automata), nodes also define a reset function $f_"reset" : S -> R -> S$. Since
+nodes may be parameterized by a value, this reset function takes in an
+additional reset parameter $R$ and the previous state, and returns an updated
+state. A discrete model with input $I$ and output $O$ is then a triple of an
+initial state, and a step and reset function:
+
+#todo[Justify resets. Why not $R$ in $f_"step"$?]
+
+$ DNode(I, O, R, S) eq.def {
+ s_0 : S;
+ f_"step" : S -> I -> O times S;
+ f_"reset" : S -> R -> S;
+ } $
+
+The simulation of such a model then defines two streams: the inner state $s$ and
+the output $o$:
+
+$ DSim(M)(i_n) = o_n "where"
+ (o_n, s_(n+1)) = M.f_"step" (i_n, s_n), s_0 = M.s_0 $
+
+A possible implementation of this simulation in #ocaml, where streams are
+represented by lists of values, is given in @lst:discrete_sim.
+
+#todo[Make figures stand out more.]
+
+#figure(
+ ```ocaml
+ (** Discrete model representation. *)
+ type ('a, 'b, 'r, 's) dnode =
+ DNode of { s0 : 's; step : 's -> 'a -> 'b * 's; reset : 's -> 'r -> 's }
+
+ (** Run a model on a list of inputs. *)
+ let dsim (DNode model) input =
+ let rec run s = function
+ | [] -> []
+ | i :: is -> let (o, s) = model.step s i in (o :: run s is) in
+ run model.s0 input
+ ```,
+ caption: [Discrete simulation in #ocaml],
+)
+
+#figure(
```zelus
- let hybrid ball () = y where
- rec der y = y' init 50.0
- and der y' = -9.81 init 0.0 reset z -> -0.8 *. (last y')
- and z = up (-. y)
+ let h = 0.01 (* Integration time step. *)
+
+ (* Forward Euler integrator. *)
+ let node integr(x0: float, x': float) = (x: float) where
+ rec x = x0 -> pre(x +. x' *. h)
+
+ let node sincos_discrete() = (sin: float, cos: float) where
+ rec sin = integr(0.0, cos) (* (dsin/dt)(t) = cos(t), sin(0) = 0.0 *)
+ and cos = integr(1.0, -. sin) (* (dcos/dt)(t) = -sin(t), cos(0) = 1.0 *)
+ ```,
+ caption: [Sine and cosine approximation in discrete #zelus],
+)
+
+As an example, consider the program in @lst:sincos_discrete, written in the
+discrete subset of #zelus (it could have been written in #lustre in a similar
+way). This program computes approximations for the sine and cosine functions.
+The `integr` node implements a forward Euler integrator. It takes as input a
+stream `x0` representing the initial value of the signal to be integrated, and a
+stream `x'` representing the derivative of this same signal, sampled at a
+predefined integration step `h`. It then defines a new stream `x`, approximating
+the integral of `x'`, as follows:
+
+$ #zel("x")_0 = #raw("x0")_0 wide
+ #zel("x")_(n+1) = #zel("x")_n + #zel("x'")_n dot #zel("h") $
+
+The `integr` node is used to compute an approximation of the solution to a
+restricted form of an initial value problem: given a function $x'(t)$ computing
+the derivative of a variable $x$ with respect to time (that is,
+$(d x)/(d t)(t) = x'(t)$), and an initial value $x_0$ for this variable,
+its solution is a function of time $x(t)$ whose derivative is $x'$, and whose
+value at $t = 0$ is $x(0) = x_0$.
+
+The `sincos_discrete` node then uses this integrator to approximate solutions
+for the sine and cosine functions. Since $(d sin)/(d t)(t) = cos(t)$ and
+$(d cos)/(d t)(t) = -sin(t)$, and $sin(0) = 0$ and $cos(0) = 1$, we can
+formulate $sin$ and $cos$ through an initial value problem. We can then use two
+integrators to approximate solutions to $sin$ and $cos$. The output of `sincos`
+at instant $n$ is then the pair of $sin(n dot h)$ and $cos(n dot h)$.
+
+#todo[
+ Emphasize that the inner state of the subnodes is contained in the state of
+ the parent node, and that separate calls to subnodes are separate instances.
+]
+
+== Continuous-Time Models
+
+#figure(
+ cetz.canvas({
+ let data = csv("data/sincos_discrete.csv")
+ let dsin = data.map(t => (float(t.at(0)), float(t.at(1))))
+ let dcos = data.map(t => (float(t.at(0)), float(t.at(2))))
+ cetz-plot.plot.plot(size: (12, 2), axis-style: "left",
+ x-tick-step: 10, y-tick-step: 1, x-grid: "both", y-grid: "both",
+ legend: (11, 2.2), y-label: none, x-label: "Time", {
+ cetz-plot.plot.add(label: "sin", dsin, style: (stroke: (dash: "solid")))
+ cetz-plot.plot.add(label: "cos", dcos, style: (stroke: (dash: "dashed")))
+ })
+ }),
+ caption: [Simulation of @lst:sincos_discrete with `h = 0.01`],
+)
+
+
+While the model in @lst:sincos_discrete is simple to understand, it is somewhat
+rigid: the integration method is fixed, as well as the time step. The simulation
+results strongly depend on these parameters. Given a time step of $0.01$, for
+instance, the approximation quickly diverges from the analytical solution, as
+seen in @fig:sincos_discrete. This divergence is simple to understand: rather
+than using the derivative of $sin$ and $cos$ exactly, we sample it at specific
+instants (multiples of $h$), and consider it to be constant between samples.
+This causes the result of an integration "step" to be slightly imprecise. Over
+time, these imprecisions accumulate, leading to this divergence. This problem is
+not exclusive to the forward Euler integrator --- it is an inherent difficulty
+of numerical approximation of ODEs.
+
+#figure(
+ ```zelus
+ let hybrid sincos() = (sin: float, cos: float) where
+ rec der sin = cos init 0.0 (* (dsin/dt)(t) = cos(t), sin(0) = 0.0 *)
+ and der cos = -. sin init 1.0 (* (dcos/dt)(t) = -sin(t), cos(0) = 1.0 *)
+ ```,
+ caption: [Sine and cosine approximation in continuous #zelus]
+)
+
+Still, we can do better. Rather than remain in the discrete world, #zelus allows
+us to express a signal as a function of _continuous time_. Time is no longer
+logical and represented by a series of discrete instants, but rather physical
+and continuous. A model is now a function of signals on physical time. Given an
+input signal of type $I$, it defines a continuously evolving inner state of type
+$S$, and an output signal of type $O$. This is represented through an initial
+state $s_0 : S$ and two functions: the derivative function
+$f_"der" : I -> S -> S'$ computes the derivative $S'$ of the inner state $S$ at
+a given time using the value of the input signal and the inner state at that
+time; the output function $f_"out" : I -> S -> O$ computes the output of the
+model at a given time given the value of the input signal and the inner state at
+that time. A continuous model is then a tuple of an initial state and of these
+two functions:
+
+#todo[
+ The fact that the arguments of $f_"der"$ are at the same time is not visible
+ in the type.
+]
+
+$ CNode(I, O, S, S') eq.def
+ { s_0: S; f_"der": I -> S -> S'; f_"out": I -> S -> O; } $
+
+For instance, the model of @lst:sincos_discrete can be expressed in continuous
+time as seen in @lst:sincos_continuous. Here, `sin` and `cos` are expressed
+directly as initial value problems. The notation ```zelus der x = e init e0```
+expresses that the derivative of `x` with respect to time is `e`, and that the
+value of `x` at time `t = 0` is `e0`.
+
+A major difference between the discrete and continuous models is that the
+description of the continuous model is kept separate from the ODE solving
+machinery. Nothing in @lst:sincos_continuous expresses any constraints for how
+the two initial value problems of `sin` and `cos` are solved -- we leave this
+detail to the language implementation. This allows for greater flexibility in
+the simulation process, because independence from the solver means we can choose
+our approximation method.
+
+== Numerical ODE Solvers
+
+The simulation of a continuous model solves the initial value problem posed by
+the initial state $s_0$ and the derivative function $f_"der"$, and uses this
+solution in order to compute the output signal with $f_"out"$. This is done
+using a numerical solver which approximates the solution, such as #sundials ---
+the `integr` node from @lst:sincos_discrete is another example of a numerical
+solver (albeit not a very good one). In general, numerical ODE solvers can be
+considered through a simple interface: given an initial value problem for a
+signal $y : Time -> Y$, in the form of a maximum time $stop$, a derivative
+function $f : [0, stop] -> Y -> Y'$ such that for all $t in [0, stop]$,
+$(d y)/(d t)(t) = f(t, y(t))$), and an initial value $y_0 : Y$ such that
+$y(0) = y_0$, a numerical ODE solver provides a function
+
+$ italic("csolve")(f)(y_0) : (h : Time) ->
+ (h' : Time) times (italic("dky") : [0, h'] -> Y) $
+
+This function, given a requested horizon $h : Time$ (this represents the date up
+to which we wish to know the approximation of the solution), returns a new
+horizon $h' <= h$ and an approximation $italic("dky") : [0, h'] -> Y$ of the
+solution to the initial value problem, that is, $italic("dky")(t) approx y(t)$
+for all $t in [0, h']$. This function is called a _dense solution_. #footnote[
+ The notation $italic("dky")$ and the name _dense solution_ are taken from the
+ #sundials interface.
+]
+
+Of particular interest is the fact that numerical ODE solvers can be considered
+to compute approximations _sequentially_. Some solvers, such as #sundials,
+perform integration in successive steps. Each step produces a part of the
+approximation, and successive calls to the step function produce successive
+parts. More formally, a single call to the $italic("csolve")$ function provides
+us with an approximation of the solution up to the returned horizon $h'$, which
+may be less than the requested date $h$. To obtain an approximation of the
+solution at a later date, we must perform another call, this time with initial
+state $italic("dky")(h')$, which is the best approximation of the value of
+$y(h')$. This new call will provide us with a new horizon $h'' >= h'$, and a new
+approximation $italic("dky")' : [h', h''] -> Y$. This is then repeated as often
+as needed to build a larger approximation of the solution.
+
+This sequential process allows a synchronous interpretation of an ODE solver as
+a discrete node. Rather than producing a single function of continuous time, an
+ODE solver is a synchronous node that takes in a stream of requested horizons
+and produces a stream of dense functions and associated horizons.
+
+$ Dense(A) eq.def { h : Time; u : [0, h] -> A } $
+
+The ODE solver, given a stream of requested horizons, produces a stream of dense
+solutions, and operates on an internal state $S$, whose definition depends on
+the solver being used. Its reset parameter is an initial value problem:
+
+$ IVP(Y, Y') eq.def { y_0 : Y; stop : Time; f : [0, stop] -> Y -> Y' } $
+
+An ODE solver can thus be considered as a particular kind of discrete node:
+
+$ CSolver(Y, Y', S) eq.def DNode(Time, Dense(Y), IVP(Y, Y'), S) $
+
+A signal of type $V$ is now represented as a stream of interval-defined
+functions, that is, a function from $NN$ to $Dense(V)$. Successive values in the
+stream are interpreted as successive intervals on the time domain. Given a
+stream of dense functions $s : NN -> Dense(V)$, the corresponding signal
+$s' : Time -> V$ is defined as
+
+$ s'(t) = cases(
+ s_0.u(t) & "if" t in [0, e_0],
+ s_n.u(t - e_(n-1)) & "if" t in (e_(n-1), e_n] "for some" n > 0,
+ "undefined" & "otherwise"
+) wide wide e_n = sum^n_(i=0)s_i.h $
+
+#todo[Find better names than $s$ and $s'$.]
+
+where $e_n$ is the stream of instants at which the solver stops. We assume dense
+functions to be continuous on their domain. However, nothing prevents
+discontinuities from occurring at the joining points of the stream, that is, for
+the stream $s$ above, we might have that $s_n.u(s_n.h) != s_(n+1).u(0)$. The ODE
+solver does not itself introduce discontinuities; the only discontinuities in
+the system are those introduced by the input signal.
+
+It is important to note that the solver approximates solutions to the _entire_
+initial value problem at once. That is, if the initial value problem is composed
+of two or more unrelated ODEs (in the sense that they operate on distinct sets
+of variables), the solver does not consider these ODEs separately; rather, it
+computes approximations to the entire system at once. This can lead to some
+unexpected behaviour. Some solvers, such as #sundials, adapt their step length
+according to the system being approximated. If given a particularly "steep"
+curve (say, a sine wave with a high frequency), the step length is shortened to
+mitigate errors; instead, if given a "gentler" curve, the step length is
+lengthened to increase efficiency. Of course, the approximation obtained depends
+on the length of the integration steps the solver performs; integrating the same
+curve with different step lengths yields different results.
+
+The addition of a new, unrelated ODE to a pre-existing system can then alter the
+results obtained for this system. If the newly added ODE is "steep", the solver
+reduces its step length to mitigate error, and computes an approximation for the
+entire system using this new step length. This differs from the steps the solver
+would have taken had the new ODE not been included; and so, the results obtained
+for the rest of the system are different. This is particularly important: the
+simultaneous integration of two unrelated systems yields different results from
+their integration in isolation.
+
+Of course, one could consider a different approximation method, where each ODE
+is integrated independently with its own ODE solver, rather than all together as
+a whole system. This is called distributed simulation, and while it solves the
+issue of interference between unrelated systems, it raises other difficulties
+and performance concerns, and is more difficult to implement. #zelus chooses
+instead to live with the consequences of using a single solver for the entire
+system.
+
+The simulation of a continuous-time system with an ODE solver is now considered
+as a synchronous node. Rather than continuous-time signals, it operates on
+streams of interval-defined functions. At each step, it takes the value provided
+by its input signal, initializes the ODE solver with an appropriate initial
+state and derivative function, and performs a step of the solver to obtain an
+approximation of the solution to the initial value problem. It then uses this
+approximation and the model's output function to build an output value.
+
+Since the ODE solver does not necessarily reach the requested horizon in a
+single step, the simulation may need to execute several steps of the underlying
+solver for each dense function provided as input. That is, for an input value
+defined on the time interval $[0, h]$, the ODE solver produces a list of
+approximations such that their "concatenation" represents the solution over the
+full interval $[0, h]$, with each item in the list being the result of a step of
+the solver. Since the ODE solver does not introduce discontinuities, it is safe
+to consider this concatenation as a single continuous function.
+
+Unfortunately, some ODE solvers (such as #sundials) work in such a way that
+stepping the solver multiple times per step of the simulation is infeasible.
+Solvers operate on an internal state (in #sundials, this is called a session),
+and the approximation returned by a step of the solver depends on this internal
+state. Some solvers rewrite this internal state in-place during a step,
+invalidating previously produced approximations. We cannot then step the solver
+multiple times, as all but the last approximation produced will be unusable. But
+stepping the solver once per simulation step would mean that the solver would
+have to store its input values until the solver is ready to work on them, which
+we cannot do either: the input stream might be the output of another simulation,
+in which case all but the latest input value would be unusable as well.
+
+To solve this conundrum, we wrap the dense functions in our stream with an
+option type, representing the "readiness" of the simulation to accept a new
+input value.
+
+$ Signal(V) eq.def Option(Dense(V)) $
+
+Rather than stepping the solver multiple times per step of the simulation, we
+step the solver once, and return the output up to the horizon reached by the
+solver. If the solver has not reached the input's horizon, we perform another
+step, giving $None$ as input to the simulation, and do so until the solver
+reaches the input's horizon, after which the simulation simply returns $None$.
+Once this occurs, it is safe to provide a new dense function as input and begin
+integration again.
+
+The simulation can then take as input as many $None$ values as necessary for the
+solver to "have time" to reach the horizon of the input. The input stream must
+then contain as many successive $None$ as needed after a dense function for the
+solver to reach the horizon requested by this dense function. The simulation
+assumes the input signal takes this form; if a new input value is provided to
+the simulation before it is done integrating the previous one, no guarantee is
+made on the correctness of the results.
+
+The simulation of a continuous-time model with a solver is then a special case
+of a discrete node, with a complex internal state $S$:
+
+$ CSim : & CNode(I, O, S_M, S'_M) -> CSolver(S_M, S'_M, S_S) ->
+ \ & DNode(Signal(I), Signal(O), Unit, Signal(I) times S_M times S_S) $
+
+A step of the simulation can take three forms, depending on its input and on
+the state of the simulation. If the input is a new dense function, we assume we
+are done integrating the previous input. We reset the solver to take into
+account the new input value in the initial value problem (the model's derivative
+function uses the input, and so the derivative function given to the solver must
+change). If the input is $None$ and we are not done integrating, we call the
+solver again and use the approximation returned to build a dense value for the
+output stream. If the input is $None$ and we are done integrating, we do
+nothing: the simulation is waiting for the next input value. A possible
+implementation in #ocaml is given in @lst:continuous_sim.
+
+#figure(
+ ```ocaml
+ (** Continuous model. *)
+ type ('i, 'o, 's, 'sder) cnode =
+ CNode of { s0: 's; fder: 'i -> 's -> 'sder; fout: 'i -> 's -> 'o }
+ (** Dense values. *)
+ type 'a dense = { h : time; u : time -> 'a }
+ (** Initial value problem. *)
+ type ('y, 'yder) ivp = { y0 : 'y; f : time -> 'y -> 'yder; h : time }
+ (** ODE solver. *)
+ type ('y, 'yder, 's) csolver = (time, 'y dense, ('y, 'yder) ivp, 's) dnode
+
+ (** Simulation of a continuous model, as a discrete node. *)
+ let csim (CNode model) (DNode solver) =
+ let s0 = (None, model.s0, solver.s0) in
+ let step (current_input, mstate, sstate) new_input =
+ match (new_input, current_input) with
+ | (Some input, None) ->
+ let ivp_f t m = model.fder (input.u t) m in
+ let ivp = { y0=mstate; f=ivp_f; h=input.h } in
+ None, (Some i, mstate, solver.reset sstate ivp)
+ | (None, Some input) ->
+ let ({h; u=dky}, sstate) = solver.step sstate input.h in
+ let u t = model.fout (input.u t) (dky t) in
+ let current_input = if h >= input.h then None else current_input in
+ Some {h; u}, (current_input, dky h, sstate)
+ | (None, None) -> None, (None, ms, ss)
+ | (Some _, Some _) -> assert false in
+ let reset (_, ms, ss) () = (None, model.y0, solver.y0) in
+ DNode { s0; step; reset }
+ ```,
+ caption: [Continuous simulation in #ocaml]
+)
+
+== Hybrid Models
+
+Continuous-time models allow for precise descriptions of physical systems and
+continuous behaviours. However, they lack the ability to describe discrete
+_events_. For instance, consider the model of a bouncing ball. We can describe
+its behaviour in the air with two ODEs for the ball's position $y$ (the distance
+from the ground) and speed $y'$:
+
+$ (d y)/(d t)(t) = y'(t) wide (d y')/(d t)(t) = -g $
+
+where $g$ is the gravitational constant ($g approx 9.81$). Coupled with an
+initial position and speed, this gives us our initial value problem, which can
+be approximated as seen above. However, nothing here describes the ball's
+bouncing behaviour as it touches the ground: it will fall until the end of time.
+We would ideally like to identify the instant at which the ball touches the
+ground, stop the simulation at this instant, and perform some changes to the
+model state to represent the impact of the bounce (say, negate the speed and
+scale it by a constant), before resuming the simulation with the updated state.
+
+The question of discrete events comes up whenever we wish to include discrete
+behaviour in a continuous model. Since time is logical in discrete nodes,
+nothing tells us when, in continuous time, we should perform discrete steps.
+There are many possible choices. We could, for instance, pick a step length $p$
+and say that the discrete step is performed periodically at every $p$. In
+practice, hybrid system modelers like #simulink and #zelus use _zero-crossings_.
+They monitor a certain value during simulation, and perform a discrete step
+whenever this value changes from strictly negative to positive or null.
+
+More formally, a zero-crossing on a function $z : [0, h] -> RR$ occurs at time
+$t in [0, h]$ if any of the following conditions are met:
+$ & (z(t - epsilon) < 0) and (z(t) > 0) \
+ or & (z(t - epsilon) < 0) and (z(t) = 0) and (z(t + epsilon) >= 0) \
+ or & (z(t - epsilon) = 0) and (z(t) = 0) and (z(t + epsilon) > 0) $
+with $epsilon in Time$ a strictly positive, solver-dependent constant
+representing the maximum precision of the zero-crossing detection mechanism
+#footnote[
+ For instance, if the zero-crossing mechanism represented time with
+ floating-point numbers, a sensible choice for $epsilon$ could be
+ `epsilon_float`.
+] (see @sec:zero_crossing_solvers).
+
+An important point is that discrete events should take _no time_ to execute. The
+physical time of the model does not change during discrete steps. This is
+similar to the approach of the synchronous languages, where the execution of a
+step is considered to be instantaneous. Additionally, multiple discrete steps
+may occur directly after one another. The time basis should reflect this: if we
+use $Time = RR_+$, successive discrete steps would occur at the same time, and
+we have no way to distinguish the order of execution, or even represent as a
+function whose codomain is a single value. In the superdense semantics of
+@cit:op_sem_hyb_sys, the time basis is the set $RR_+ times NN$, ordered
+lexicographically ($(t_1, n_1) < (t_2, n_2)$ iff $t_1 < t_2$ or
+$t_1 = t_2 and n_1 < n_2$) #footnote[
+ This is not the only possible choice; for instance, in
+ @cit:nonstd_sem_hyb_sys_mod, the semantics of hybrid systems is expressed
+ using non-standard analysis, and the time basis is the set of hyperreals
+ $\ ^*RR$.
+]. At each physical instant $t : RR_+$, any number $n$ of discrete steps may
+occur in successive logical instants $(t, 0), (t, 1), ..., (t, n)$. In our
+stream representation of signals, discrete instants are instead represented by
+dense functions with horizon $h = 0$ (that is, defined on the interval $[0,0]$).
+The order of execution of successive discrete steps is simply the order given by
+the stream #footnote[
+ The interpretation of a stream of dense functions $s : NN -> Dense(V)$ as a
+ function of superdense time $s' : (RR_+ times NN) -> V$ is then defined
+ through the following recursive functions:
+ $ f(t, n, s, i) & = cases(
+ s_i.u(t) & "if" t < s_i.h,
+ g(s_i.u(0), n-1, s, i+1) & "if" t = s_i.h,
+ f(t - s_i.h, n, s, i+1) & "otherwise"
+ ) \
+ g(v, n, s, i) & = cases(
+ v & "if" n = 0,
+ s_i.u(0) & "if" s_i.h != 0,
+ g(s_i.u(0), n - 1, s, i + 1) & "otherwise"
+ ) \
+ s'(t, n) & = f(t, n, s, 0) $
+ The stream representation is quite similar to the hybrid sequences of
+ @cit:theory_timed_io_automata[sec. 3.4].
+].
+
+A hybrid model describes such systems whose behaviour goes through both
+continuous and discrete phases. Its state $S$ contains both discrete and
+continuous parts: the continuous part $Y$ is defined by ODEs, and evolves
+during continuous phases, while the discrete part is only modified during
+discrete steps, and must be constant during continuous phases #footnote[
+ This restriction is enforced by typing: see @cit:zelus_sync_lng_with_ode for
+ more details.
+]. The model defines functions $c_"get" : S -> Y$ and $c_"set" : S -> Y -> S$ to
+get and set the continuous state $Y$ from the whole state $S$. To handle
+zero-crossings, a model with input signal $I$ defines a zero-crossing function
+$f_"zero" : S -> I -> Y -> Zo$, where $Zo$ is a vector of values to be monitored
+for zero-crossings. The inner state $S$ also maintains a vector of boolean flags
+$Zi$, representing the events corresponding to the values in $Zo$ (a flag is set
+to true if its corresponding event has occured), and a function
+$z_"set" : S -> Zi -> S$ to update the state when an event has been detected.
+
+For implementation reasons, a hybrid model also defines two additional functions
+$f_"jump" : S -> BB$ and $f_"horizon" : S -> Time$. The $f_"horizon"$ function
+allows a model to provide a horizon after which no further integration must
+occur. This is used to indicate whether or not the simulation, after a discrete
+step, must perform another discrete step (the horizon is 0 in this case); and
+for the ```zelus period``` construct, which represents a recurring discrete
+event.
+// #footnote[
+// The ```zelus period``` construct could also be implemented using a sawtooth
+// ode:
+// ```zelus period(p)``` is roughly equivalent to \
+// ```zelus let rec der t = 1.0 init -. p reset z -> -. p and z = up(t) in z```
+// ].
+The $f_"jump"$ function is used to indicate whether or not a discrete step
+has introduced a discontinuity in the model's state: if so, the solver must be
+reset to take this change into account. These two functions exist mainly for
+implementation purposes: see @sec:sim_algorithm for more details.
+
+Finally, a hybrid model defines all functions required by discrete and
+continuous models:
+
+$ HNode(I, O, R, S, Y, Y', Zi, Zo) eq.def {
+ & s_0 : S; \
+ & c_"get" : S -> Y;
+ c_"set" : S -> Y -> S;
+ z_"set" : S -> Zi -> S; \
+ & f_"step" : S -> I -> Zi -> O times S;
+ f_"der" : S -> I -> Y -> Y'; \
+ & f_"out" : S -> I -> Y -> O;
+ f_"zero" : S -> I -> Y -> Zo; \
+ & f_"reset" : S -> R -> S;
+ f_"horizon" : S -> Time;
+ f_"jump" : S -> BB
+} $
+#todo[This definition is ugly.]
+
+#zelus provides several ways to specify zero-crossing events, of which the
+```zelus up(e)``` construct is the most common. It monitors its subexpression
+`e` for zero-crossings, and triggers an event whenever a zero-crossing occurs on
+`e`. Constructs like ```zelus present``` and ```zelus reset``` allow models to
+execute discrete behaviour when an event is triggered. These constructs are
+compiled down to an internal representation quite similar to $HNode$ (see
+@cit:sync_based_codegen_hyb_sys_lng for more details). Continuous-time models
+are a special case of hybrid ones, where the discrete step function does not
+do anything, and the zero-crossing function does not monitor any signals; and
+#zelus compiles them down to the same internal representation.
+
+Going back to our bouncing ball, we monitor the expression $-y$ for
+zero-crossings. Whenever this expression becomes positive, the ball touches the
+ground. When this zero-crossing event is triggered, we represent the effect of
+the bounce by negating the speed, so that the ball starts moving up again, and
+decreasing it by a small factor, to represent the loss of inertia from the
+collision. A possible implementation in #zelus is given in @lst:bouncing_ball.
+Whenever the zero-crossing event `z` occurs, `y'` is negated and scaled down;
+the notation ```zelus last y'``` represents the left-limit of `y'`.
+
+#figure(
+ ```zelus
+ let hybrid ball(y0, y'0) = y where
+ rec der y = y' init y0
+ and der y' = -9.81 init y'0 reset z -> -0.8 *. last y'
+ and z = up(-. y)
```,
caption: [The bouncing ball in #zelus]
-)
+)
+
+== Zero-crossing Detection
-More formally, and as done in @cit:alg_ana_hyb_sys, a hybrid system $cal(H)$ can
-be defined as a graph whose nodes represent continuous behaviour, and whose
-edges represent discrete transitions:
-$ cal(H) = (L o c, V a r, E d g, A c t, I n v, I n i) $
-where:
-- $L o c$ is a finite set of _locations_;
-- $V a r$ is a finite set of _variables_;
-- $E d g subset.eq L o c times F times L o c$ is a finite set of _transitions_
+The monitoring of the zero-crossing expressions $Zo$ requires a mechanism to
+detect zero-crossings, termed a zero-crossing solver. Multiple methods exist
+(#zelus uses the Illinois method @cit:illinois), but the zero-crossing solver
+can be summarized as providing a function
+$ italic("zsolve") : (Time -> Y -> Zo) -> Dense(Y) -> Time times Option(Zi) $
-== Executing models
+taking as input a zero-crossing function and a dense solution $v$, and returning
+a pair of a horizon $h in [0, v.h]$ and an optional vector of boolean flags $z$,
+such that if the zero-crossing solver detects one or more zero-crossing events,
+$h$ is the earliest instant at which a zero-crossing occurs, and $z$ is not
+null; otherwise, $h = v.h$, and $z$ is null.
-To execute such a model, we first compile it into a synchronous function, as
-described in @cit:sync_based_codegen_hyb_sys_lng. The details of this
-compilation step are not particularly relevant to our purposes, and can be
-ignored. What is more interesting is the output of this compilation step: a
-single synchronous function. The simulation loop is then itself described as a
-synchronous function operating on
+Similarly to the ODE solver, a zero-crossing solver can be seen as a synchronous
+node. Since the zero-crossing function does not change during continuous
+behaviour (it takes as argument the continuous part of the state, and the
+discrete part is considered constant during integration), it may be used as a
+reset parameter; the input is a stream of dense solutions, and the output is a
+stream of pairs of reached horizons and optional zero-crossings.
-#MENTION("Use of a single solver")
+$ ZSolver\(Y, Zi, Zo, S) eq.def
+ DNode\(Dense(Y), Time times Option(Zi), Time -> Y -> Zo, S) $
-#pagebreak()
+#todo[Finish this section.]
-// The compilation of a hybrid model into a synchronous function is described in
-// detail in @cit:sync_based_codegen_hyb_sys_lng, but can be summarized quite
-// succintly as follows. By pairing this synchronous function with an
-// off-the-shelf ODE solver like #sundials, we can then simulate the dynamics of
-// the system. This is done by repeatedly performing execution steps according to
-// two different modes: discrete and continuous.
+A zero-crossing solver may be combined with an ODE solver to obtain the full
+solver mechanism used by the simulation of a hybrid system. This full solver
+both performs approximation of the solution to the initial value problem of the
+model, as well as zero-crossing detection using this approximation. That is,
+given an ODE solver $c s : CSolver(Y, Y', S_C)$ and a zero-crossing solver
+$z s : ZSolver(Y, Zi, Zo, S_Z)$, their composition takes the form of a new
+synchronous node $s : Solver(Y, Y', Zi, Zo, S_C times S_Z)$, where
-// The continuous mode operates as follows: we first call the ODE solver in order
-// to approximate the dynamics of the model's continuous state.
+$ Solver(Y, Y', Zi, Zo, S) eq.def
+ DNode(Time, Dense(Y) times Option(Z_i), IVP(Y, Y') times ZCP(Y, Zo), S) $
-// Continuous steps first call the ODE solver to approximate the dynamics of the
-// model's continuous variables. The solver will return a function defined on a
-// time interval, which we then provide as input to the zero-crossing solver, which
-// will monitor the evolution of zero-crossing values along this interval. After
-// both solvers have been called, we choose what the next step's mode will be:
-// - if no zero-crossings have been detected, we output the entire solution
-// provided by the ODE solver, and the next step remains continuous;
-// - if a zero-crossing occurs, we return the solution provided by the ODE solver
-// up to the zero-crossing instant, and the next step becomes a discrete step.
+A step of the full solver $s$ takes in a horizon $h : Time$. It then performs a
+step of the underlying ODE solver $c s$ with input $h$, obtaining a dense
+function $v$ defined up to a horizon $h' <= h$ approximating the solution of the
+initial value problem, and uses this dense function as input to the
+zero-crossing solver $z s$, which returns a new horizon $h'' <= h'$ and an
+optional zero-crossing event $z$. The final horizon $h''$ is then used as the
+horizon of $v$ (since $v$ is defined on $[0, h']$ and $h'' <= h'$, then $v$ is
+defined on $[0, h'']$). Finally, it returns the dense function $v$ with horizon
+$h''$ and the optional zero-crossing event $z$.
-// Discrete steps perform state changes and side effects. We first call the model's
-// step function, which updates the state and outputs a value. We then decide what
-// the next step is. If a zero-crossing event occured due to the current step, the
-// next step is another discrete step. If no new event occured, we perform a
-// continuous step.
+$ s.f_"step" (h) eq.def (v', z) "where"
+ v = c s.f_"step" (h)", "
+ (h', z) = z s.f_"step" (v)
+ "and "v' = { h=h', u=v.u } $
-// Executing such a model is quite simple. There are two modes of execution:
-// discrete and continuous. In continuous mode, we call the ODE solver in order
-// to obtain an approximation of the variables defined through ODEs, and monitor
-// for zero-crossings. If a zero-crossing occurs, we enter the discrete mode, in
-// which we perform computation steps as needed, until no other zero-crossing
-// occurs, in which case we go back to the continuous mode, and repeat, as seen
-// in @automaton.
+The full solver mechanism is then paired with a hybrid model to construct a
+node representing the simulation of the model.
-// #figure(finite.automaton(
-// (D: (D: "cascade", C: "no cascade"),
-// C: (C: "no zero-crossing", D: "zero-crossing")),
-// initial: "D", final: (), layout: finite.layout.linear.with(spacing: 3)
-// ), caption: [High-level overview of the runtime], placement: top)
+== The Simulation Algorithm
-= Runtime
-To solve this issue, we need to redefine what the runtime of our hybrid system
-models is. The runtime is the piece of software that handles the pairing of a
-hybrid model with a solver and the different execution phases that need to take
-place.
+The simulation of a full hybrid model $m : HNode(I, O, R, S_M, Y, Y', Zi, Zo)$
+with a solver $s : Solver(Y, Y', Zi, Zo, S_S)$
+can also be seen as a synchronous node, operating on streams of dense functions.
+Simulation steps take two forms: discrete steps perform state changes and side
+effects, and continuous steps approximate the solution to the initial value
+problem of the model and monitors for zero-crossing events. The simulation
+alternates between these two modes as needed, switching from continuous to
+discrete steps if a zero-crossing event occurs, and from discrete to continuous
+steps if no additional discrete steps are necessary. A high-level overview of
+the simulation's behaviour is given in @fig:sim_automaton.
-To allow for independent simulation of the assertions, we need to simulate them
-with their own ODE solvers. Since assertions can themselves contain ODEs, they
-can be viewed as hybrid systems themselves. Assertions can also contain their
-own sub-assertions, and so recursively: our runtime needs to handle this
-accordingly.
+#figure(
+ finite.automaton(
+ (D: (D: "cascade", C: "no cascade"),
+ C: (C: "no zero-crossing", D: "zero-crossing")),
+ initial: "D", final: (), layout: finite.layout.linear.with(spacing: 3)
+ ),
+ caption: [
+ Overview of the simulation loop; $D$ and $C$ represent discrete
+ and continuous step modes
+ ]
+)
-== Hybrid models with assertions
-A hybrid model with assertions can be defined using a recursive data type:
-```ocaml
-type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a =
- HNodeA : {
- body : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec;
- assertions : ('p, 's, unit, 'y, 'yder, 'zin, 'zout) hnode_a list;
- } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a
+The simulation's internal state stores five things: the internal states
+$s_m : S_M$ and $s_s : S_S$ of the model and solver, respectively; the current
+simulation $m o d e$ (either idle, discrete or continuous); a boolean flag $r$
+indicating whether we should reset the solver before the next continous step
+(see @sec:full_sim); the current input $i : Option(Dense(I))$; and the current
+simulation time with respect to the input's domain $n o w in [0, i.h]$, used in
+discrete steps to obtain the correct input. We use the same trick as with
+continuous-time models to solve the problem of the ODE solver taking several
+steps to integrate a single input value: we ask that the input stream contains
+enough successive $None$ values for the ODE solver to finish integrating the
+current input value (as explained in @sec:ode_solvers).
+
+$ SimState(S_M, S_S, I) = {
+ s_m : S_M;
+ s_s : S_S;
+ m o d e : M o d e;
+ r : BB;
+ i : Option(Dense(I));
+ n o w: [0, i.h]
+ } $
+
+The simulation of a hybrid system is then a discrete node on streams:
+
+$ HSim : & HNode(I, O, R, S_M, Y, Y', Zi, Zo) ->
+ Solver(Y, Y', Zi, Zo, S_S) -> \
+ & DNode(Signal(I), Signal(O), R, SimState(S_M, S_S, I)) $
+
+Its step function's behaviour varies depending on the simulation mode; the
+following sections describe these behaviours in more details.
+
+=== Discrete steps
+
+A discrete step occurs whenever a zero-crossing event is triggered or a new
+value is obtained by the simulation. Zero-crossing events may be triggered in
+two different ways: either by detection using the zero-crossing solver during a
+previous continuous step, or resulting from the action of a previous discrete
+step as indicated by $f_"horizon"$. New input values require a discrete step to
+be performed in order to reset the underlying solver. Discrete steps may modify
+the entire model state, and perform side effects. The simulation's physical time
+does not advance during a discrete step.
+
+#note[Should we explain cascades in detail?]
+
+A discrete step of the simulation simply calls the model's $f_"step"$ function
+with the appropriate inputs, constructs a dense function defined on $[0,0]$
+using the output (as described in @sec:hybrid_mod, this represents a discrete
+step), and updates the simulation state as needed for the next step. Four
+possible situations arise. If the model's indicated horizon (obtained with the
+$f_"horizon"$ function defined by the model) requires us to perform another
+discrete step, we keep the simulation in discrete mode. If no other discrete
+step must be performed, but we are done integrating the current input (the
+current time $n o w$ is greater than or equal to the current input's horizon
+$i.h$), we cannot proceed further, and must wait for additional input; and so we
+switch to idle mode. If no other discrete step must be performed and we have not
+reached the input's horizon, we switch to continuous mode. In this case, we may
+have to reset the solver. This occurs if the current discrete step has caused a
+discontinuity (as indicated by the model's $f_"jump"$ function), or if the
+state's reset flag $r$ is set, in which case we build out a new initial value
+problem and zero-crossing problem using the current model state and reset the
+solver.
+
+A possible implementation of the discrete step in #ocaml is given in
+@lst:sim_step_discrete.
+
+#figure(
+ ```ocaml
+ (** Simulation mode. *)
+ type mode = Discrete | Continuous | Idle
+ (** Simulation state. *)
+ type ('i, 'sm, 'ss) sim_state =
+ { sm: 'sm; ss: 'ss; mode: mode; r : bool; i: 'i dense option; now: time }
+ (** Zero-crossing problem. *)
+ type ('y, 'zo) zcp = { f : time -> 'y -> 'zo; y0 : 'y }
+ (** Hybrid model. *)
+ type ('i, 'o, 'r, 's, 'y, 'yd, 'zi, 'zo) hnode = HNode of {
+ s0 : 's; cget : 's -> 'y; cset : 's -> 'y -> 's; zset : 's -> 'zi -> 's;
+ horizon : 's -> time; jump : 's -> bool; step : 's -> 'i -> 'o * 's;
+ fder : 's -> 'i -> 'y -> 'yd; fzer : 's -> 'i -> 'y -> 'zo;
+ fout : 's -> 'i -> 'y -> 'o; reset : 's -> 'r -> 's;
+ }
+
+ let dstep (HNode model) (DNode solver) state =
+ let i = Option.get state.i in
+ let (o, sm) = model.step state.sm (i.u state.now) in
+ let out = { u=(fun _ -> out); h=0.0 } in
+ let state =
+ let h = model.horizon sm in
+ if h <= 0 then { state with sm }
+ else if state.now >= i.h then { state with mode=Idle; i=None; sm }
+ else if model.jump sm || state.r then (* Reset solver. *)
+ let fder t y = model.fder sm (i.u t) y in
+ let fzer t y = model.fzer sm (i.u t) y in
+ let ivp = { f=fder; h=i.h; y0=model.cget sm } in
+ let zcp = { f=fzer; y0=model.cget sm } in
+ let ss = solver.reset (ivp, zcp) state.ss in
+ { state with mode=Continuous; sm; ss; r=false }
+ else { state with mode=Continuous; sm } in
+ (Some out, state)
+ ```,
+ caption: [Discrete simulation step in #ocaml]
+)
+
+=== Continuous steps
+
+Continuous steps advance time, approximate the solution to the model's initial
+value problem using the ODE solver, and monitor the model's zero-crossing
+function for zero-crossing events using the zero-crossing solver. They operate
+on a restricted part of the model's state; only the continuous part of the state
+may be modified, and the discrete part is constant. Furthermore, no side-effects
+may occur during continuous steps. These restrictions are enforced by a typing
+pass during the compilation process @cit:sync_based_codegen_hyb_sys_lng.
+
+A continuous step performs a call to the solver to obtain both an approximation
+of the solution to the model's initial value problem and an optional
+zero-crossing event. It builds a dense function representing the output on the
+approximation's domain using the model's $f_"out"$ function. Then, it updates
+the simulation state as needed for the next step. Once again, four possible
+situations arise. If a zero-crossing event has occured, we must perform a
+discrete step, and so we switch to discrete mode and update the model's state to
+take into account the zero-crossing event with the $z_"set"$ function. If no
+zero-crossing event has occurred, but we have reached the end of the current
+input (the horizon reached by the solver is greater than or equal to the current
+input's horizon), we must perform a discrete step as well, and so we switch to
+discrete mode. If no zero-crossing event has occured and we have not reached the
+current input's horizon, but we have reached the model's desired stopping point
+(as indicated by the model's $f_"horizon"$ function), we must again perform a
+discrete step, and so we switch to discrete mode. Otherwise, we can continue
+integrating; we keep the simulation mode as continuous.
+
+A possible implementation of the continuous step in #ocaml is given in
+@lst:sim_step_continuous.
+
+#figure(
+ ```ocaml
+ let cstep (HNode model) (DNode solver) state =
+ let i = Option.get state.i in
+ let stop = min (model.horizon state.sm) i.h in
+ let (({ h=now; u=dky }, z), ss) = solver.step state.ss stop in
+ let sm = model.cset state.sm (dky now) in
+ let out = { h=now; u=fun t -> model.fout sm (i.u t) (dky t) } in
+ let state = match z with
+ | Some z ->
+ let sm = model.zset sm z in
+ { state with mode=Discrete; sm; ss; now }
+ | None ->
+ if model.horizon sm <= 0.0 || now >= i.h
+ then { state with mode=Discrete; sm; ss; now }
+ else { state with mode=Continuous; sm; ss; now } in
+ (Some out, state)
+ ```,
+ caption: [Continuous simulation step in #ocaml]
+)
+
+=== Complete definition
+
+The full step function then performs the correct kind of step depending on the
+simulation mode; if the mode is Idle, it simply returns $None$ and the
+unmodified state. When a new dense function is provided as input, it updates the
+current input and time, sets the reset flag, and switches to discrete mode. Once
+again, it expects the input stream to contain as many successive $None$ values
+as needed for the solver to integrate the entire input, as explained in
+@sec:ode_solvers.
+
+The simulation's reset function simply resets the model using its $f_"reset"$
+function, sets the simulation mode to Idle, and sets the reset flag $r$ in the
+simulation state, so that the next discrete step resets the solver before
+integration.
+
+Finally, its initial state is simply the initial states of the model and solver,
+the mode set to idle, an empty current input ($None$) and a current time set at
+0\.
+
+Its implementation in #ocaml is given in @lst:sim_algorithm.
+
+#figure(
+ ```ocaml
+ let hsim (HNode model) (DNode solver) =
+ let s0 = { mode=Idle; i=None; now=0.0; ss=model.s0; ss=solver.s0; r=true } in
+ let step state i = match (i, state.mode) with
+ | Some _, Idle ->
+ let state = { state with mode=Discrete; i; now=0.0; r=true } in
+ dstep (HNode model) (DNode solver) state
+ | None, Discrete -> dstep (HNode model) (DNode solver) state
+ | None, Continuous -> cstep (HNode model) (DNode solver) state
+ | None, Idle -> (None, state) in
+ | Some _, _ -> assert false
+ let reset r state =
+ { state with mode=Idle; sm=model.reset r state.sm; r=true } in
+ DNode { s0; step; reset }
+ ```,
+ caption: [Simulation of a hybrid model in #ocaml]
+)
+
+== Implementation Details
+
+While the above algorithm works, it suffers from a few flaws which limit its
+efficiency. Indeed, resetting the solver at every new input value is
+counterproductive. ODE solvers with adaptive step lengths (such as #sundials)
+begin integration by performing very small steps in time, and increase their
+step length later, as they obtain more information on the function they are
+currently integrating. Resetting the solver slows down the progress of the
+simulation, as such ODE solvers will perform shorter steps than if they had not
+been reset.
+
+If two successive input values can be considered to be continuous (that is, the
+second one only extends the first, with no discontinuity at the joining point),
+there is no particular reason why we should reset the solver. This occurs for
+instance between successive continuous steps. The ODE solver does not itself
+introduce discontinuities, and so the simulation should not reset its solver if
+taking as input two successive continuous steps of an ODE solver. As a first
+solution, we can equip our dense values with an additional bit of information
+$c$ representing whether the next value in the stream is simply an extension of
+themselves:
+
+$ Dense(V) eq.def { h : Time; u : [0, h] -> V; c : BB } $
+
+When building the output value, the simulation knows how this output value
+behaves compared to the next one: if we just performed a continuous step and the
+next step is also continuous, we know that the next output value will simply be
+an extension of the current one, and we can include this information in the
+current output value. In all other cases, it is safe to consider that successive
+values are discontinuous. When a simulation receives an input value, it can then
+reset the solver only if necessary.
+
+Another issue comes from the impossibility of stepping the solver more than once
+per step of the simulation, as seen in @sec:ode_solvers. Since adaptive solver
+begin with small integration steps, output values will be defined on small
+intervals. If we compose simulations together, the resulting output will be
+defined on smaller and smaller intervals, even though this is not always
+necessary, as the ODE solvers do not introduce discontinuities between their
+steps.
+
+We can impose another restriction on our ODE solvers to mitigate this issue. If
+the solver provides a function to copy its internal state, allowing us to
+preserve the validity of the previous approximations, we can safely step the
+solver multiple times per step of the simulation and concatenate the results. A
+discrete node with state copies operating on a state $S$ defines an additional
+function $f_"copy" : S -> S$ returning a copy of the inner state, which may be
+used for the rest of the simulation. Previously computed approximations then
+depend on the original copy of the state, which remains untouched by later
+steps.
+
+$ DNodeC(I, O, R, S) eq.def {
+ s_0 : S;
+ f_"step" : S -> I -> O times S;
+ f_"reset" : S -> R -> S;
+ f_"copy" : S -> S
+ } $
+
+Given a solver with state copies, the simulation can then perform multiple steps
+of the solver, performing a state copy in between each step and concatenating
+the approximations returned by the solver. This concatenation
+$j o i n : Dense(V) times Dense(V)-> Dense(V)$ is defined as expected:
+
+$ j o i n(l, r) = {
+ h = l.h + r.h;
+ u = lambda t. "if" t <= l.h "then" l.u(t) "else" r.u(t + l.h)
+ } $
+
+This does not free us from the option type in our signals, however; the
+simulation may still produce more than one dense function as output per dense
+function as input (for instance, if a zero-crossing event occurs).
+
+== Lifting the Runtime
+
+An interesting consequence of interpreting simulations as discrete nodes is that
+we can reasonably consider manipulating them directly inside the language. This
+takes the form of a module in #zelus' standard library, providing several
+primitives and utility functions to create and manipulate simulations and
+signals. For instance, the function
+
+```zelus
+ val solve : ('i -C-> 'o) -S-> ('i signal -D-> 'o signal)
```
-Notice that the list of assertions takes as input the inner state of the parent
-model: assertions are checked on the model state, and any variable which is
-required by the assertion becomes a state variable.
+takes as input a continuous-time model (this is indicated by the `-C->` arrow)
+from `'i` to `'o` and producing a discrete-time node (indicated by the `-D->`
+arrow) from `'i signal` to `'o signal` #footnote[
+ Due to some restrictions in the higher-order facilities of #zelus, the
+ argument to `solve` must be statically known (that is, we must be able to
+ allocate the necessary memory before starting the execution of the program);
+ this is represented by the `-S->` arrow.
+]. It represents the simulation _with a dedicated instance_ of an ODE solver of
+its argument; that is, the ODE solver used to simulate the argument of `solve`
+is separate from the one used for the rest of the program. We can now choose
+which parts of our program we wish to simulate in isolation from the rest. The
+primitives `compose` and `synchr`, whose signatures are
+#grid(columns: (1fr, 1fr),
+ ```zelus
+ val compose :
+ ('a signal -D-> 'b signal) -S->
+ ('b signal -D-> 'c signal) -S->
+ ('a signal -D-> 'c signal)
+ ```,
+ ```zelus
+ val synchr :
+ ('a signal -D-> 'b signal) -S->
+ ('a signal -D-> 'c signal) -S->
+ ('a signal -D-> ('b * 'c) signal)
+ ```
+)
+allow for the composition and synchronization of independent simulations.
+Indeed, standard composition of simulations does not work: the output of the
+first would not take into account that the second simulation requires $None$
+values as input until it is done integrating its current input. We need to
+ensure that this requirement is met. This is simple: to simulate
+$f circle.small g$, step $g$ and then $f$ when $f$ is done integrating (i.e.
+when $f$ returns $None$), otherwise only step $f$. Furthermore, two simulations
+will not necessarily advance at the same speed, and a synchronization primitive
+is required to simulate two systems in parallel. Its behaviour is simple: step
+only the simulation that has not progressed as far as the other one, and return
+the solution only on the interval on which both are defined.
-== Solvers as synchronous nodes
-== Simulations as synchronous nodes
-#MENTION("the new runtime")
+A possible implementation of both of these in #ocaml is given in
+@sec:additional_code, but it is interesting to note that both of these could
+just as well be implemented in discrete #zelus, given the right tools to
+manipulate dense functions (in particular, getting their horizon and splitting
+them into two smaller, successive dense functions). One could even imagine that
+simulations themselves are implemented in #zelus\; they are only discrete nodes
+implementing a particular automaton, and discrete #zelus provides all of the
+necessary constructs to implement automata. Given an interface allowing us to
+instanciate and call a solver, this seems entirely feasible.
-= Assertions
-#MENTION("how assertions are done")
+= Hybrid Observers and Assertions
+
+#todo[
+ Describe assertions as a run-time defensive construct (as in #ocaml), and
+ their equivalent in #lustre: synchronous observers. Emphasize that the
+ synchronous subset does not cause any problem, since the observer has no
+ impact on its underlying model. Describe their naïve extension to continuous
+ time, and their unexpected impact on the behaviour of the continuous model:
+ recall the examples of Section 1.1. Conclusion: we need to simulate continuous
+ observers with their own solver.
+]
+
+During the design and implementation of software, programmers often use
+assertions in order to check certain properties, both before and during
+execution. In #ocaml, for instance, the ```ocaml assert``` instruction checks
+that a certain expression evaluates to ```ocaml true``` at runtime, and raises
+an exception otherwise. This is an example of a defensive use of assertions as a
+way to avoid unwanted behaviour: the programmer chooses to suspend execution
+rather than proceed with a state which they consider invalid.
+
+An important property of assertions is that they are _transparent_: their
+presence does not affect the result of the computation. Of course, this property
+requires the expression evaluated in the assertion to not cause any visible
+side-effects: in #ocaml, if the assertion modifies mutable data or performs I/O
+operations, executing the assertion is not transparent. Still, if the
+subexpression respects certain criteria, we can safely assume that the presence
+of the assertion will not affect the final result.
+
+In synchronous languages like #lustre, the classical way to represent such an
+assertion is an _observer_: a node whose sole purpose is to monitor its input
+stream to check a certain property. @lst:observer_discrete gives an example of
+an assertion in discrete #zelus, and an idealized translation of this assertion
+into a separate observer node. The `assertion_f` node monitors its input stream,
+and warns the user if a property is not respected (here, that a certain value,
+obtained by integrating another stream with the `integr` node of
+@lst:sincos_discrete, is always positive); otherwise it has no effect. In
+particular, if its input stream respects the property, the result of the
+simulation does not change whether the call to `assertion_f` in `f` is included
+or not.
+
+#figure(
+ placement: top,
+ grid(
+ columns: (100fr, 100fr),
+ align: horizon,
+ ```zelus
+ let node f (*...*) =
+ let v = (*...*) in
+ assert (
+ let p = integr(0.0, v) in
+ p >= 0.0
+ ); (*...*)
+ ```,
+ ```zelus
+ let node assertion_f(v) =
+ let p = integr(0.0, v) in
+ p >= 0.0
+
+ let node f (*...*) =
+ let v = (*...*) in
+ (match assertion_f(v) with
+ | true -> print_endline "ERROR"
+ | false -> ());
+ (*...*)
+ ```,
+ ),
+ caption: [A discrete assertion and its idealized translation as an observer. ]
+)
+
+In general, a discrete assertion on a model $M : DNode(I, O, R, S)$ can be seen
+as another node whose input stream is the state $S$ of the parent model, and
+whose output stream is a Boolean value; it defines its own internal state, with
+its own local variables, subnodes, and so on and so forth. The parent model then
+calls the assertion with its internal state during each step.
+
+In continuous-time, one could imagine a similar way of encoding such behaviour.
+@lst:observer_continuous presents a possible implementation of the behaviour of
+@lst:observer_discrete in continuous time. Rather than using inequalities, which
+are considered discontinuous and are not allowed in continuous code, we use a
+zero-crossing event to monitor the signal, and perform the appropriate side
+effect if the event occurs. The integration is once again handled by the
+simulation, as described in @sec:sim_algorithm.
+
+#figure(
+ placement: top,
+ ```zelus
+ let hybrid assertion_f(v) =
+ let der p = v init 0.0 in
+ up(-. p)
+
+ let hybrid f (*...*) =
+ let v = (*...*) in
+ (present(assertion_f(v)) ->
+ print_endline "ERROR"
+ else ());
+ (*...*)
+ ```,
+ caption: [A naïve implementation of a continuous observer.]
+)
+
+Unfortunately, this implementation does not meet our criteria for assertions.
+Indeed, adding ODEs to a model changes the approximation returned by the ODE
+solver, as explained in @sec:ode_solvers, even if the new ODEs are entirely
+unrelated to the existing ones. The implementation in @lst:observer_continuous
+does not separate the body of the assertion from its parent model, and the
+simulation runs both the model and its assertion at the same time. Therefore,
+the assertion may impact the results of its parent model, and is not
+transparent.
+
+We wish instead to simulate the assertion independently from its parent model,
+with its own ODE solver.
+
+== Models With Assertions
+
+#todo[
+ Since we are able to transform a hybrid model into a discrete one (albeit with
+ a complicated type), we can isolate parts of a model in order to simulate them
+ with a dedicated solver, and provide the obtained stream as input to another
+ hybrid model simulation, which performs the assertion on this precomputed
+ solution. This gives rise to a new datatype for a model together with its
+ assertions.
+]
+
+Since an assertion can be considered as a separate model operating on the inner
+state of its parent model, we can represent a model with assertions as a pair
+of the parent model $m$, operating on its inner state $S_M$, and a list of
+assertion models. All assertions operate on the same input datatype $S_M$ (the
+state of the parent model) and return Boolean output signals $BB$.
+
+$ HNodeA(I, O, R, S_M, Y, Y', Zi, Zo) eq.def {
+ & m : HNode(I, O, R, S_M, Y, Y', Zi, Zo); \
+ & a : List(HNodeA(S_M, BB, R_A, S_A, Y_A, Y'_A, Zi_A, Zo_A)) } $
+
+Note that the output signal is Boolean even in continuous time. There are
+two possible interpretations of this: either assertions benefit from relaxed
+typing rules for their output, and are allowed a limited subset of discrete
+behaviours in continuous time; or assertions in continuous time are defined in
+terms of zero-crossing events, and as such the output of the assertion will be
+constant during continuous phases (in fact, the output will be constantly true,
+as the simulation would not have entered a continuous step otherwise). Both
+interpretations lead to the same updated simulation algorithm, as we will see in
+@sec:updated_sim.
+
+The $HNodeA$ datatype is recursive: indeed, nothing prevents assertions from
+containing their own assertions, and so on and so forth.
+
+== Updated Simulation
+
+#todo[
+ Present a version of the simulation with assertions (probably the simple
+ version with a single assertion, rather than the complex one with polymorphic
+ recursion and such).
+]
+
+#todo[
+ Maybe present a high-level overview of the compilation passes? This seems a
+ little out of scope.
+]
+
+= Future Work
+
+#todo[
+ What do we do next? Depends on how much the #zelus compiler is already able to
+ do: compiling the source language into the internal representation with an
+ additional model is a good next step, compiling the source language into a
+ complex, recursive datatype, etc...
+]
+
+= Conclusion
+
+#todo[
+ Conclude.
+]
#pagebreak()
-= Annex
-#set heading(level: 2)
-#bibliography("sources.bib", full: true)
-#set heading(level: 1)
+
+#bibliography(
+ "sources.bib",
+ // full: true,
+ style: "alphanumeric.csl",
+ title: [Bibliography]
+)
+
+= Appendix
+
+#todo[Find a proper bibliographic style.]
+
+== Additional Code
+
+#set figure(placement: none)
+#figure(
+ ```ocaml
+ let compose (DNode f) (DNode g) =
+ let s0 = (f.s0, g.s0) in
+ let step (sf, sg) = function
+ | Some i ->
+ let (v, sg) = g.step sg (Some i) in
+ let (o, sf) = f.step sf v in
+ (o, (sf, sg))
+ | None ->
+ let (o, sf) = f.step sf None in
+ match o with Some _ -> (o, (sf, sg)) | None ->
+ let (o, sg) = g.step sg None in
+ match o with None -> (None, (sf, sg)) | Some _ ->
+ let (o, sf) = f.step sf o in
+ (o, (sf, sg)) in
+ let reset (sf, sg) (rf, rg) = (f.reset sf rf, g.reset sg rg) in
+ DNode { s0; step; reset }
+ ```,
+ caption: [Composition of simulations in #ocaml]
+)
+
diff --git a/doc/sources.bib b/doc/sources.bib
index 4c8f0a4..3b10a08 100644
--- a/doc/sources.bib
+++ b/doc/sources.bib
@@ -41,7 +41,6 @@
modelers).},
}
@inbook{cit:op_sem_hyb_sys,
- address = {Berlin, Heidelberg},
series = {Lecture Notes in Computer Science},
title = {Operational Semantics of Hybrid Systems},
volume = {3414},
@@ -128,10 +127,88 @@
author = {Alur, Rajeev and Courcoubetis, Costas and Halbwachs, Nicolas and
Henzinger, Thomas A and Ho, P-H and Nicollin, Xavier and Olivero,
Alfredo and Sifakis, Joseph and Yovine, Sergio},
- journal = {Theoretical computer science},
+ journal = {Theoretical Computer Science},
volume = {138},
number = {1},
pages = {3--34},
year = {1995},
publisher = {Elsevier},
}
+@inproceedings{cit:lustre,
+ title = {LUSTRE: A declarative language for programming synchronous systems},
+ author = {Pilaud, Daniel and Halbwachs, N and Plaice, J.A.},
+ booktitle = {Proceedings of the 14th Annual ACM Symposium on Principles of
+ Programming Languages (14th POPL 1987). ACM, New York, NY},
+ volume = {178},
+ pages = {188},
+ year = {1987},
+ organization = {Citeseer},
+}
+@article{cit:sundials,
+ title = {SUNDIALS: Suite of nonlinear and differential/algebraic equation
+ solvers},
+ author = {Hindmarsh, Alan C and Brown, Peter N and Grant, Keith E and Lee,
+ Steven L and Serban, Radu and Shumaker, Dan E and Woodward, Carol S
+ },
+ journal = {ACM Transactions on Mathematical Software (TOMS)},
+ volume = {31},
+ number = {3},
+ pages = {363--396},
+ year = {2005},
+ publisher = {ACM New York, NY, USA},
+}
+@inproceedings{cit:sundialsml,
+ title = {Sundials/ML: interfacing with numerical solvers},
+ author = {Bourke, Timothy and Inoue, Jun and Pouzet, Marc},
+ booktitle = {ACM Workshop on ML},
+ year = {2016},
+}
+@book{cit:theory_timed_io_automata,
+ address = {Cham},
+ series = {Synthesis Lectures on Distributed Computing Theory},
+ title = {The Theory of Timed I/O Automata},
+ rights = {https://www.springernature.com/gp/researchers/text-and-data-mining
+ },
+ ISBN = {978-3-031-00875-7},
+ url = {https://link.springer.com/10.1007/978-3-031-02003-2},
+ DOI = {10.1007/978-3-031-02003-2},
+ abstractNote = {This monograph presents the Timed Input/Output Automaton
+ (TIOA) modeling framework, a basic mathematical framework to
+ support description and analysis of timed (computing)
+ systems. Timed systems are systems in which desirable
+ correctness or performance properties of the system depend on
+ the timing of events, not just on the order of their
+ occurrence. Timed systems are employed in a wide range of
+ domains including communications, embedded systems, real-time
+ operating systems, and automated control. Many applications
+ involving timed systems have strong safety, reliability and
+ predictability requirements, which makes it important to have
+ methods for systematic design of systems and rigorous
+ analysis of timing-dependent behavior. An important feature
+ of the TIOA framework is its support for decomposing timed
+ system descriptions. In particular, the framework includes a
+ notion of external behavior for a timed I/O automaton, which
+ captures its discrete interactions with its environment. The
+ framework also defines what it means for one TIOA to implement
+ another, based on an inclusion relationship between their
+ external behavior sets, and defines notions of simulations,
+ which provide sufficient conditions for demonstrating
+ implementation relationships. The framework includes a
+ composition operation for TIOAs, which respects external
+ behavior, and a notion of receptiveness, which implies that a
+ TIOA does not block the passage of time.},
+ publisher = {Springer International Publishing},
+ author = {Kaynar, Dilsun K. and Lynch, Nancy and Segala, Roberto and
+ Vaandrager, Frits},
+ year = {2011},
+ collection = {Synthesis Lectures on Distributed Computing Theory},
+ language = {en},
+}
+@article{cit:illinois,
+ title = {Inverse interpolation, a real root of f (x)= 0},
+ author = {Snyder, J.N.},
+ journal = {University of Illinois Digital Computer Laboratory, ILLIAC I
+ Library Routine H1-71},
+ volume = {4},
+ year = {1953},
+}
diff --git a/doc/zelus.sublime-syntax b/doc/zelus.sublime-syntax
index f67a028..44022a5 100644
--- a/doc/zelus.sublime-syntax
+++ b/doc/zelus.sublime-syntax
@@ -6,20 +6,30 @@ scope: source
contexts:
main:
- - match: \b(let|in|where|rec|and|local)\b
+ - match: \b(val|let|in|where|rec|and|local)\b
scope: keyword.control
- - match: \b(if|then|else)\b
+ - match: \b(if|then|else|match|with)\b
scope: keyword.control.conditional
- match: \b(hybrid|node)\b
scope: keyword.control
- - match: \b(up|assert|der|init|reset|last)\b
- scope: entity.name.constant
+ - match: \b(period|up|assert|der|init|reset|present|last|fby|pre)\b
+ scope: keyword.control
+
+ - match: \b(float|int)\b
+ scope: entity.name.type
+
+ - match: \b[0-9]*(\.[0-9]+)?\b
+ scope: constant
+
+ - match: \b(true|false)\b
+ scope: constant
- match: '"'
push: string
+
- match: \(\*
push: comment