Skip to content

Commit 8f2c1cc

Browse files
mekisemekise
committed
Write dynamics, ode solver and files to run steady state
1 parent c2d9196 commit 8f2c1cc

File tree

6 files changed

+153
-38
lines changed

6 files changed

+153
-38
lines changed

‎.gitattributes

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
notebooks/* linguist-vendored
1+
*.ipynb linguist-vendored

‎plot.ipynb

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": 1,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"import numpy as np\n",
10+
"import matplotlib.pyplot as plt"
11+
]
12+
},
13+
{
14+
"cell_type": "code",
15+
"execution_count": 14,
16+
"metadata": {},
17+
"outputs": [],
18+
"source": [
19+
"data = np.load(\"src/run.npz\")"
20+
]
21+
},
22+
{
23+
"cell_type": "code",
24+
"execution_count": 16,
25+
"metadata": {},
26+
"outputs": [
27+
{
28+
"data": {
29+
"text/plain": [
30+
"[<matplotlib.lines.Line2D at 0x7fe0907570d0>]"
31+
]
32+
},
33+
"execution_count": 16,
34+
"metadata": {},
35+
"output_type": "execute_result"
36+
},
37+
{
38+
"data": {
39+
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfHUlEQVR4nO3deXBc5Znv8e/T3VLL1mLJtiRbixfwKrzJFkuYBAgxYAzBYjEFZKWSopwJqak792Yge2ZIZpKbydQkEwiXm5Dt3hmGzcYQA2FzuAMBLG94N8LGtrxIMl4k2Wh/7x/dsmW5ZbXkbp1efp+qLuss6n5eGv10dPqc5zXnHCIikvx8XhcgIiKxoUAXEUkRCnQRkRShQBcRSREKdBGRFKFAFxFJEQGvXnjs2LFu0qRJXr28iEhSWrt27WHnXGGkbQMGupk9CtwINDjnZkXYbsDPgMXASeCLzrl1Az3vpEmTqKmpGWg3ERHpxcz29LctmlMuvwUWnWP79cDU8OMe4JeDKU5ERGJjwEB3zr0OHDnHLkuA37uQt4B8MxsfqwJFRCQ6sfhQtBTY12u5LrzuLGZ2j5nVmFlNY2NjDF5aRER6xCLQLcK6iA1inHOPOOeqnHNVhYURz+mLiMgQxSLQ64DyXstlwIEYPK+IiAxCLAJ9JfB5C7kMOO6cOxiD5xURkUGI5rLF/wCuAsaaWR3wPSADwDn3MLCK0CWLtYQuW7w7XsXGUmdXN+1d3bR3hh5tnaeX2zpPr2/v6mLm+DzGjxrhdckiIuc0YKA75+4cYLsDvhqzigaw9UATKzbsPzN4u7pp6+g6I6AjhXXvfboH0QZ+1IgMXv7bKynMDcZvYCIi58mzO0WHau+RE/zhL3vIDPhCD7+PYM/X4eXMgI+crEB4vf/UumCfffp+Hex5vgwfmX4/mQEfLW0dLPvDOv7+2S384q75Xg9fRKRfSRfoi2aNZ9Gs4b3M/d6rp/AvL+2kel49CyuKh/W1RUSipeZcUVh25YVML87l2ys209za4XU5IiIRKdCjkBnw8aNbZ1Pf3MqPX9judTkiIhEp0KNUOaGAuy+fzP95ay9rPjhXJwQREW8o0Afhv187jdL8Edz31Lu0dnR5XY6IyBkU6IOQHQzwj7fMZlfjCR58rdbrckREzqBAH6QrpxVyS2Upv1z9PtsONnldjojIKQr0IfjOjRWMGpHB/U+9S9dg7lASEYkjBfoQFGRn8t1PV7Cx7ji/eWO31+WIiAAK9CG7aW4JV88o4qd/2sm+Iye9LkdERIE+VGbGA9Wz8Bl8c/kmQi1tRES8o0A/D6X5I7jv+hn8v/cO8/S6/V6XIyJpToF+nj576UQWTCzggT9u5XBLm9fliEgaU6CfJ5/P+NEtsznZ1sXfP7vV63JEJI0p0GNganEuX/3kFJ7deIBXttV7XY6IpCkFeox85aoLmVaco46MIuIZBXqMhDoyzuFQUyv/84UdXpcjImlIgR5D8ycU8MXLJ/GHt/ZQo46MIjLMFOgx9j+una6OjCLiCQV6jGUHA/zw5lm833iCh9SRUUSGkQI9Dq6aXsTNlaU8tPp9th9SR0YRGR4K9Dj5zo0V5I3I4L6nNqkjo4gMCwV6nIzOzuR7n65g475j/O7ND7wuR0TSgAI9jm6aW8Inpxfykxd3qCOjiMSdAj2OzIwf3Dwbn8G3VmxWR0YRiSsFepyV5o/g7xbN4PWdjSxfr46MIhI/CvRh8LnLJjJ/Qj7/8Jw6MopI/CjQh4HPZ/z41jmcbOviH9SRUUTiRIE+THo6Mq7ceIBXt6sjo4jEXlSBbmaLzGyHmdWa2f0Rto8ys2fNbKOZbTGzu2NfavI71ZFx+WZa2jq9LkdEUsyAgW5mfuBB4HqgArjTzCr67PZVYKtzbi5wFfBTM8uMca1Jr6cj48GmVn7ywnavyxGRFBPNEfolQK1zbpdzrh14DFjSZx8H5JqZATnAEUCHoBHMn1DAFz42id+/tYe1e9SRUURiJ5pALwX29VquC6/r7RfATOAAsAn4G+dcd98nMrN7zKzGzGoaGxuHWHLy+/p10ykZNYL7ntpEW6c6MopIbEQT6BZhXd87ZK4DNgAlwDzgF2aWd9Y3OfeIc67KOVdVWFg4yFJTR09HxtqGFh587X2vyxGRFBFNoNcB5b2Wywgdifd2N/C0C6kFdgMzYlNiaurpyPjL1bXsONTsdTkikgKiCfQ1wFQzmxz+oPMOYGWfffYCnwIws2JgOrArloWmou/cWEFuVgb3PfWuOjKKyHkbMNCdc53AvcCLwDbgcefcFjNbZmbLwrs9AFxuZpuAV4D7nHOH41V0qujpyLhh3zF+/5cPvC5HRJKcedUwqqqqytXU1Hjy2onEOcfdv13DO7uP8Kf/dgVlBSO9LklEEpiZrXXOVUXapjtFPWZm/PDm2QB8c7k6MorI0CnQE0Bp/gj+7rrpvL6zkRUb1JFRRIZGgZ4gPvexSaGOjM9u5UN1ZBSRIVCgJwh/uCNjc2snD63WtekiMngK9AQytTiXay8qZvn6/bR3nnWjrYjIOSnQE8zSBeUcOdGuFrsiMmgK9ATzialjKc4L8kRNndeliEiSUaAnmIDfxy3zy3htRwMNTa1elyMiSUSBnoCWLiij28HTmlRaRAZBgZ6ALijMoWpiAY/X7NONRiISNQV6grq9qpxdjSdYt/eY16WISJJQoCeoxXPGMyLDzxM1+wbeWUQEBXrCygkGuGHOeJ579yAn2zWbn4gMTIGewG6vKqelrZPnNx3yuhQRSQIK9AR28aQCJo0ZyeM67SIiUVCgJzAzY2lVOW/vPsKeD094XY6IJDgFeoK7ZX4pPoMn1+rOURE5NwV6ghs/agSfmFrIU2vrNO+oiJyTAj0JLK0q48DxVt6o1TStItI/BXoSuKaimPyRGTyh0y4icg4K9CQQDPhZMreEF7cc4vjJDq/LEZEEpUBPEkurymnv7GblRjXsEpHIFOhJYlbpKGaOz+Nx9UkXkX4o0JPI7VVlbNp/nG0Hm7wuRUQSkAI9iVTPKyXT79NsRiISkQI9iRRkZ7KwoogVGzSJtIicTYGeZJZWaRJpEYlMgZ5krphaSHFeUB+OishZFOhJxu8zbp1fxuodDdRrEmkR6SWqQDezRWa2w8xqzez+fva5ysw2mNkWM/tzbMuU3m7rmUR6na5JF5HTBgx0M/MDDwLXAxXAnWZW0WeffOAh4Cbn3EXA0tiXKj0uKMzh4kkFPLFWk0iLyGnRHKFfAtQ653Y559qBx4Alffa5C3jaObcXwDnXENsypa+lC3omkT7qdSkikiCiCfRSoPeUOXXhdb1NAwrMbLWZrTWzz8eqQIls8ZzxjMz065p0ETklmkC3COv6/p0fABYANwDXAd8xs2lnPZHZPWZWY2Y1jY2Ngy5WTssJBlg8ezzPbjygSaRFBIgu0OuA8l7LZcCBCPu84Jw74Zw7DLwOzO37RM65R5xzVc65qsLCwqHWLGG3V5Vzor1Lk0iLCBBdoK8BpprZZDPLBO4AVvbZ5xngE2YWMLORwKXAttiWKn1pEmkR6W3AQHfOdQL3Ai8SCunHnXNbzGyZmS0L77MNeAF4F3gH+JVzbnP8yhbQJNIicqaorkN3zq1yzk1zzl3onPtheN3DzrmHe+3zE+dchXNulnPuX+NUr/ShSaRFpIfuFE1ymkRaRHoo0FPA7VXlmkRaRBToqWBhRRH5IzP04ahImlOgp4BgwE/1vFL+tLWeYyfbvS5HRDyiQE8Rty0oC08i3fcWARFJFwr0FDGrdBQV4/PUCkAkjSnQU8hSTSItktYU6ClEk0iLpDcFegopyM7kmopiTSItkqYU6CnmtqoyTSItkqYU6CnmiqmFjMvL0iTSImlIgZ5i/D7jlvmlmkRaJA0p0FPQ0qpyTSItkoYU6Clo8thsTSItkoYU6ClqaZUmkRZJNwr0FHXDbE0iLZJuFOgpKjsY4AZNIi2SVhToKWxpeBLpVZpEWiQtKNBTWM8k0k+oT7pIWlCgpzBNIi2SXhToKe7W+WWaRFokTSjQU9y4UVlcMa2QJzWJtEjKU6CngaULyjmoSaRFUp4CPQ1oEmmR9KBATwOaRFokPSjQ08TSKk0iLZLqFOhp4qISTSItkuoU6Gnkdk0iLZLSFOhpZIkmkRZJaQr0NNIzifTy9XWaRFokBUUV6Ga2yMx2mFmtmd1/jv0uNrMuM7stdiVKLC2tKuPoyQ5e2aZJpEVSzYCBbmZ+4EHgeqACuNPMKvrZ78fAi7EuUmLnE+FJpJ9QKwCRlBPNEfolQK1zbpdzrh14DFgSYb+vAU8BDTGsT2LM7zNuXaBJpEVSUTSBXgr0vsWwLrzuFDMrBW4GHj7XE5nZPWZWY2Y1jY2Ng61VYuS2BZpEWiQVRRPoFmFd3y5P/wrc55zrOtcTOececc5VOeeqCgsLoyxRYm3y2GwumTSaJ2o0ibRIKokm0OuA8l7LZUDf2w2rgMfM7APgNuAhM6uORYESH7dVlbHrsCaRFkkl0QT6GmCqmU02s0zgDmBl7x2cc5Odc5Occ5OAJ4G/ds6tiHWxEjs9k0g/vkYfjoqkisBAOzjnOs3sXkJXr/iBR51zW8xsWXj7Oc+bS2LqmUT6uXcPMLkwm6lFOUwtyqWsYAQ+X6SzbCKS6Myrc6hVVVWupqbGk9eWkJ31zXzpd2vYd+SjU+uyMnxMKcphWlEuU4pDIT+tOIeygpH4FfQinjOztc65qojbFOhy/KMOahuaea++hfcaWthZ30xtQwsHj5++rDEYCAX91KIcphbnnvp3wmgFvchwOlegD3jKRVLfqBEZLJg4mgUTR5+xvqm1g9qGFmrrQyH/XkML7+w+wooNpz8Tzwz4uLAwFPTTinOYEj6inzB6JAG/OkuIDCcFuvQrLyuD+RMKmD+h4Iz1za0dvN944tSR/Hv1zazdc/SMXuuZfh8XFGafPpoPH9FPGqOgF4kXBboMWm5WBvPK85lXnn/G+hNtnbzf2MLO+hbeC5/C2bDvKM/2Cvq/mjKG//vly4a5YpH0oECXmMkOBphTls+csvwz1p9s7+T9hhM8ta6O3775AdsPNTFjXJ43RYqkMP3tK3E3MjPA7LJRfO3qKQR8xor1mgZPJB4U6DJsxuQEuWJaISs37Ke7Wy0HRGJNgS7DqrqylAPHW3nngyNelyKSchToMqyumVlMdqafFevV6VEk1hToMqxGZPq5btY4/rjpIK0d52zOKSKDpECXYVc9r5Tm1k5W79BcKCKxpECXYXf5hWMozA3qaheRGFOgy7AL+H18ek4Jr25v4PjJDq/LEUkZCnTxxM2VpbR3dfP85oNelyKSMhTo4olZpXlcUJjNcl3tIhIzCnTxhJlx87xS3t59hP3HPhr4G0RkQAp08cySeaUArNygD0dFYkGBLp6ZMGYkCyYW8MwGnXYRiQUFuniqel4J2w81s+1gk9eliCQ9Bbp46oY5JaEOjDpKFzlvCnTx1OjsTK6cVsjKDQfUgVHkPCnQxXPVlaUcPN7K27vVgVHkfCjQxXML1YFRJCYU6OK5EZl+Fs0az6rN6sAocj4U6JIQqitLaG7t5LXt6sAoMlQKdEkIl184NtSBUVe7iAyZAl0Sgt9n3DS3hNe2N6oDo8gQKdAlYfR0YFylDowiQ6JAl4RxUUkeF6oDo8iQKdAlYZgZ1fNKeUcdGEWGJKpAN7NFZrbDzGrN7P4I2z9jZu+GH2+a2dzYlyrpoKcDoxp2iQzegIFuZn7gQeB6oAK408wq+uy2G7jSOTcHeAB4JNaFSnro6cC4Yv1+nFMrAJHBiOYI/RKg1jm3yznXDjwGLOm9g3PuTefc0fDiW0BZbMuUdFJdWcrO+ha2HWz2uhSRpBJNoJcC+3ot14XX9edLwPORNpjZPWZWY2Y1jY2N0VcpaeWG2eMJ+EynXUQGKZpAtwjrIv4tbGafJBTo90Xa7px7xDlX5ZyrKiwsjL5KSSujszO5anohz2w4QJc6MIpELZpArwPKey2XAWfNGWZmc4BfAUuccx/GpjxJV0vmlXKoqZW3d+t/JZFoRRPoa4CpZjbZzDKBO4CVvXcwswnA08DnnHM7Y1+mpJuFM4vJCQbUgVFkEAYMdOdcJ3Av8CKwDXjcObfFzJaZ2bLwbt8FxgAPmdkGM6uJW8WSFkZk+rnuonE8v+mQOjCKRCkQzU7OuVXAqj7rHu719ZeBL8e2NEl3N1eW8tS6Ol7d3sDi2eO9Lkck4elOUUlYH7twDEW5QZ12EYmSAl0S1qkOjDsaOHay3etyRBKeAl0SWnVlKR1djlWbDnldikjCU6BLQruoJI8pRTk67SISBQW6JLRQB8YS3vngCHVHT3pdjkhCU6BLwjvdgfGs+9lEpBcFuiS88tEjqVIHRpEBKdAlKVRXlvJeQwtbDzZ5XYpIwlKgS1I43YFRp11E+qNAl6RQkJ3JVdOLeGbDfnVgFOmHAl2SRnVlCfVNbby9Sx0YRSJRoEvS6OnAuFzXpItEpECXpJGV4WfRrHG8sFkdGEUiUaBLUrm5spTmtk5e2dbgdSkiCUeBLknlsgvCHRg136jIWRToklT8PmPJvBJWqwOjyFkU6JJ0lswLdWD846aDXpciklAU6JJ0LirJY6o6MIqcRYEuScfMqK4sZc0HR9l3RB0YRXoo0CUp3TS3BICVG9UKQKSHAl2SUvnokVw8qYDl6sAocooCXZJWdWUptQ0tbDmgDozx1N3teG17A2v3HFEfnQQX8LoAkaG6YfZ4vr9yC89s2M+s0lFel5OS1u89yvdXbmFj3XEAxuZk8qkZxSysKObjU8YyItPvcYXSmwJdklb+yJ4OjAe4//qZ+H3mdUkpo6GplR+9sJ2n1+2nKDfIPy+dS2bAx8tb61m1+SD/WbOPrAwfH59SyLUVxVw9s4ixOUGvy057CnRJatXzSnlpaz1v7fqQv5oy1utykl5bZxe/eeMD/u2V9+jocnzlqgv56ienkBMMRcVNc0to7+zmnd1HeHlbPS9treflbfWYwfwJBVxTUczCmcVMKcrxeCTpybz6QKmqqsrV1NR48tqSOlo7urj4By9z3axx/PPSuV6Xk9Re3V7PPzy7lQ8+PMnCmUV8+4YKJo3NPuf3OOfYdrCZl7bW89K2Q2zeH/o844Kx2aFwryhm/oQC/fUUQ2a21jlXFXGbAl2S3def2Mjzmw9R8+2FZGXonO5gvd/YwgPPbWX1jkYuKMzmuzdWcNX0oiE914FjH/HKtnr+FP6rqaPLMSY7k6tnFLGwophPTB3LyEydGDgfCnRJaW/WHuauX73NL+6q5MY5JV6XkzSaWzv4t1dr+c0bu8kK+PmbhVP5/McmkRmIzcVvza0d/HlnIy9tree17Q00tXYSDPj4+JSxXFNRzKdmFlOYq/Pug3WuQNevSkl6l14whuK8ICvWH1CgR6G72/HUujp+/MIOPjzRxtIFZXz9uhkxD9fcrAxunFPCjXNK6OjqZs3uI7wUPu/+yvYGzDYxrzyfayqKuSZ83t1Mp2bOR1RH6Ga2CPgZ4Ad+5Zz7UZ/tFt6+GDgJfNE5t+5cz6kjdImlf1y1jUf/azdrvrWQguxMr8tJWOv3HuX7z25l475jVE7I5/ufvoi55fnDWoNzju2Hmnl5az0vbavn3fAlkZPGjDz1oeqCiQUE/KG/FDq7umnt7Kato4u2zm5aw/+2hdf1t+301120dZz+t/XUcmhda0c3Ab9RnJtFcV6Q4lFZFOdmMW5UeDkvi5xgIGF+2ZzXKRcz8wM7gWuAOmANcKdzbmuvfRYDXyMU6JcCP3POXXqu51WgSyxtOXCcG37+X/ygehafvWyi1+UknIamVn78wg6eWldHUW6Q+6+fQfW8UnwJ8GHloeOtp66Y+cv7H9Le1c2IDD9m0NbZfd43MwUDPoIBH1kZfoIZPoIB/+nl8LaOLkd9UyuHmlppbu086zlGZvoZl5dFUV6QcXlZFJ/xCIV+UV6QYCD+n+Gc7ymXS4Ba59yu8JM9BiwBtvbaZwnwexf67fCWmeWb2XjnnPqbyrCoGJ/HtOJQB0YF+mntnd385o3d/Dx8GeKyKy/k3qtPX4aYCMaNyuKzl03ks5dNpKWtk9d3NrLmgyP4zQhm+MgKnA7irEiB3Hdbr32CAd+gj6xPtndS39RGfVNrr0cbh5paaWhqZe3eo9Q3tdHe2X3W947OzqQoNxg6us/NCh3t9/klMCY7M26/SKN5V0uBfb2W6wgdhQ+0TymgQJdhYWYsmVfKT17cwb4jJykfPdLrkjz36vZ6HnhuG7sPn4j6MkSv5QQDLJ49nsWzx3tWw8jMAJPHBph8jv9WzjmOneygvjkU9vXHW08d4ff8Mth6oInDLW30/QMj4DP++qoL+dtrp8e89mgCPdKvkr5/A0WzD2Z2D3APwIQJE6J4aZHoLZlXwk9e3MEzG/Zz79VTvS7HM7vClyG+Fr4M8bd3XzzkyxAlMjOjIDuTguxMZozrf7/Orm4Ot7SfCvuG8L+VEwriUlc0gV4HlPdaLgP69iyNZh+cc48Aj0DoHPqgKhUZQFnBSC6ZNJrl6/fz1U9OSZgPsYZLc2sHv3i1lkff2E0w4Odbi2fyhctjdxmiDF7A72PcqNAHrMNx21s0gb4GmGpmk4H9wB3AXX32WQncGz6/filwXOfPxQvVlaV8c/kmthxoinvDrpPtnWzcd5x1e4+yds9R3mtoJn9EJoW5QQpzghTmBinK6/V1bhaFucGYN7TqfRni4ZY2bq+Kz2WIkvgGDHTnXKeZ3Qu8SOiyxUedc1vMbFl4+8PAKkJXuNQSumzx7viVLNK/xbPH8b2Vm1mxPrYdGJ1z7D/2EWv3HGXdnqOs23uMrQebTl2BMaUoh3nlBbS0dlDf1Mrm/ccjnj+F0Hniwtzg6UfOmcHfE/6jszMHvGW+72WIv/5C1bBfhiiJQ3eKSsq55/c1bNh3jL9841ND7iHS3tnNlgPHQwEePgKvb2oDQpewzSvPZ8HEAuZPLKCyPJ/8kWdf+97V7Th6sp2GpjYaW9pobG6jobmVxua2sx7NbWdfKuczGJNzduAX5QYZmxtk9Y5GnlxbR2FukPsXzeDmysS4DFHiS3eKSlqprizlT+Frmj8+NboOjIdb2li35yhr94aOwDfWHT91WVpZwQguu2BMKMAnFDBjXO6pm17Oxe8zxuYEo2or+1F7VyjcW04HfkOvwG9obmP7wWYOt7TRGT7sz/BbQl6GKN7R/wWScq6eUURuMMDy9fsjBnpXt2NnffOp0ydr9x5lz4ehyaYz/T5mlebx+csmnjoCL87LinvNIzL9TBgzkgljzn25ZXe349hHHTQ0t1IwMnNYapPkoUCXlJOV4ef62eNYtekQP2ifRUd3Nxv2Hjt1+mT93mO0hE9xjM0JsmBiPp+5dAILJhZwUcmohO7Y6PMZo7MzGa32BhKBAl1SUnVlKY/X1LHwX/7MgeMf4VzonPT0cXlUV5awYGIBCyaMpnz0iLS7vFFSlwJdUtJlk8fwqRlFdHQ7bq8qZ8HEAuaWjyI3K8Pr0kTiRoEuKcnnM379xYu9LkNkWOkWMhGRFKFAFxFJEQp0EZEUoUAXEUkRCnQRkRShQBcRSREKdBGRFKFAFxFJEZ61zzWzRmAPMAo43mvTuZZ7vh4LHI5RKX1fb6j79bc90vpoxth3W7qMuffXsRpztOONZl+Nuf/1Q/lZhuQZ82Df477LsRrzROdcYcQtzjlPH8Aj0S73fA3UxOv1h7pff9sjrY9mjOk65j5fx2TM0Y5XYz6/MQ/lZzmZxjzY93g4xtz3kQinXJ4dxHLfbfF4/aHu19/2SOsHM8Z0G7OX441mX425//XJ8rMczb7RvJ+R1g33mM/g2SmX82FmNa6fGTtSlcacHjTm9BCvMSfCEfpQPOJ1AR7QmNODxpwe4jLmpDxCFxGRsyXrEbqIiPShQBcRSREKdBGRFJFygW5mM83sYTN70sy+4nU9w8HMqs3sf5vZM2Z2rdf1DAczu8DMfm1mT3pdSzyZWbaZ/S78/n7G63qGQ7q8t73F7Gc4Hhe3n8eNAY8CDcDmPusXATuAWuD+KJ/LB/za6zEN85gL0nDMT3o9nniOH/gc8Onw1//pde3D+Z4n43sbgzGf18+w54PuM5grgPm9/wMAfuB94AIgE9gIVACzgef6PIrC33MT8CZwl9djGq4xh7/vp8B8r8c0zGNOuh/6QY7/G8C88D7/7nXtwzHmZH5vYzDm8/oZTqhJop1zr5vZpD6rLwFqnXO7AMzsMWCJc+6fgBv7eZ6VwEoz+yPw73Es+bzFYsxmZsCPgOedc+viXPJ5i9X7nKwGM36gDigDNpDEp0gHOeatw1xeXAxmzGa2jRj8DCfD/yClwL5ey3XhdRGZ2VVm9nMz+1/AqngXFyeDGjPwNWAhcJuZLYtnYXE02Pd5jJk9DFSa2TfiXdww6G/8TwO3mtkvGYZbx4dZxDGn4HvbW3/vc0x+hhPqCL0fFmFdv3dDOedWA6vjVcwwGeyYfw78PH7lDIvBjvlDIFl/eUUScfzOuRPA3cNdzDDpb8yp9t721t+YY/IznAxH6HVAea/lMuCAR7UMF405PcbcWzqOX2OO8ZiTIdDXAFPNbLKZZQJ3ACs9rineNOb0GHNv6Th+jTnWY/b6k+A+nwr/B3AQ6CD0m+xL4fWLgZ2EPh3+ltd1aswas8avMSfimNWcS0QkRSTDKRcREYmCAl1EJEUo0EVEUoQCXUQkRSjQRURShAJdRCRFKNBFRFKEAl1EJEUo0EVEUsT/B5wJNLMXx7wFAAAAAElFTkSuQmCC",
40+
"text/plain": [
41+
"<Figure size 432x288 with 1 Axes>"
42+
]
43+
},
44+
"metadata": {
45+
"needs_background": "light"
46+
},
47+
"output_type": "display_data"
48+
}
49+
],
50+
"source": [
51+
"plt.figure()\n",
52+
"plt.semilogx(data[\"T\"], data[\"S\"][:,2])"
53+
]
54+
},
55+
{
56+
"cell_type": "code",
57+
"execution_count": null,
58+
"metadata": {},
59+
"outputs": [],
60+
"source": []
61+
}
62+
],
63+
"metadata": {
64+
"kernelspec": {
65+
"display_name": "Python 3.10.5 64-bit",
66+
"language": "python",
67+
"name": "python3"
68+
},
69+
"language_info": {
70+
"codemirror_mode": {
71+
"name": "ipython",
72+
"version": 3
73+
},
74+
"file_extension": ".py",
75+
"mimetype": "text/x-python",
76+
"name": "python",
77+
"nbconvert_exporter": "python",
78+
"pygments_lexer": "ipython3",
79+
"version": "3.10.5"
80+
},
81+
"orig_nbformat": 4,
82+
"vscode": {
83+
"interpreter": {
84+
"hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
85+
}
86+
}
87+
},
88+
"nbformat": 4,
89+
"nbformat_minor": 2
90+
}

‎src/Dynamics.jl

Lines changed: 21 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -1,35 +1,21 @@
1-
function diffeqsolver(N, Δt, J::GenericSD, noise::Noise, Cω::Coupling, distro=Normal(0., 1/sqrt(Δt)))
2-
3-
end
4-
5-
6-
7-
# def M_sim(nsx, nsy, nsz, matrix):
8-
# w0, Gamma, A = prm
9-
# matrix_c_w_2 = matrix @ np.transpose(matrix)
10-
# bx_int = intpl.interp1d(tva_extended,
11-
# b(barTv, S0, prm, Nsam, whatNoise, nsx)[:, 0],
12-
# bounds_error='False',
13-
# fill_value=0.0)
14-
# by_int = intpl.interp1d(tva_extended,
15-
# b(barTv, S0, prm, Nsam, whatNoise, nsy)[:, 0],
16-
# bounds_error='False',
17-
# fill_value=0.0)
18-
# bz_int = intpl.interp1d(tva_extended,
19-
# b(barTv, S0, prm, Nsam, whatNoise, nsz)[:, 0],
20-
# bounds_error='False',
21-
# fill_value=0.0)
22-
# bn = lambda t: matrix @ np.array([bx_int(t), by_int(t), bz_int(t)])
23-
# def system(t, V, res):
24-
# s = V[0:3]
25-
# p = V[3:6]
26-
# x = V[6:9]
27-
# Beff = np.sign(gam) * (Bn + matrix_c_w_2 @ x + bn(t))
28-
# res[0] = s[1] * Beff[2] - s[2] * Beff[1]
29-
# res[1] = s[2] * Beff[0] - s[0] * Beff[2]
30-
# res[2] = s[0] * Beff[1] - s[1] * Beff[0]
31-
# res[3:6] = -(w0**2) * x - Gamma * p + np.sign(gam) * A * s
32-
# res[6:9] = p
33-
# resa = odeint(system, tva, inspi, rtol=ode_tol, atol=ode_tol)
34-
# M = np.sign(gam) * resa.values.y[:, 0:3]
35-
# return M
1+
function diffeqsolver(s0, tspan, J::LorentzianSD, bfields, matrix::Coupling; S0=1/2, Bext=[0, 0, 1])
2+
u0 = [s0[1], s0[2], s0[3], 0, 0, 0, 0, 0, 0]
3+
Cω2 = matrix.C*transpose(matrix.C)
4+
bn = t -> matrix.C*[bfields[1](t), bfields[2](t), bfields[3](t)];
5+
function f(du, u, p, t)
6+
s = @view u[1:3] # not allocating values. no hard copy.
7+
v = @view u[4:6]
8+
w = @view u[7:9]
9+
du[1:3] = cross(s, Bext + bn(t)/S0 + Cω2*v)
10+
du[4:6] = w
11+
du[7:9] = -J.ω0^2*v -J.Γ*w +J.α*s
12+
end
13+
prob = ODEProblem(f, u0, tspan)
14+
sol = solve(prob, Vern7(), abstol=1e-8, reltol=1e-8, maxiters=Int(1e7))
15+
s = zeros(length(sol.t), 3)
16+
for n in 1:length(sol.t)
17+
s[n,:] = sol.u[n][1:3]
18+
end
19+
sinterp = t -> sol(t)[1:3]
20+
return sinterp, sol.t, s
21+
end

‎src/SpiDy.jl

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ using DifferentialEquations
55
using FFTW
66
using Random
77
using Distributions
8+
using Interpolations
89

910
include("Noise.jl");
1011
include("SpectralDensity.jl");
@@ -15,8 +16,8 @@ include("Dynamics.jl");
1516
# export external files to build documentation
1617
export Noise, SpectralDensity, StochasticField, CouplingTensor, Dynamics
1718
# export structures to build documentation
18-
export ClassicalNoise, QuantumNoise, NoZeroQuantumNoise, GenericSD, LorentzianSD
19+
export ClassicalNoise, QuantumNoise, NoZeroQuantumNoise, GenericSD, LorentzianSD, AnisoCoupling, IsoCoupling
1920
# export modules and functions to build documentation
2021
export psd, bfield, spectrum, sd, sdoverω, reorgenergy, kernel, diffeqsolver
2122

22-
end
23+
end

‎src/movert.jl

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
include("./SpiDy.jl")
2+
using .SpiDy
3+
4+
using Statistics
5+
6+
using ProgressBars
7+
using NPZ
8+
9+
Δt = 0.15
10+
N = 72_000
11+
tspan = (0., N*Δt)
12+
13+
J = LorentzianSD(1., 7., 5.);
14+
15+
matrix = AnisoCoupling([-sin(π/4) 0. 0.
16+
0. 0. 0.
17+
cos(π/4) 0. 0.]);
18+
19+
T = 10 .^ LinRange(-3, 2, 12)
20+
Sss = zeros(length(T), 3)
21+
22+
navg = 12
23+
24+
Threads.@threads for n in ProgressBar(1:length(T))
25+
noise = ClassicalNoise(T[n]);
26+
s = zeros(navg, 3)
27+
for k in 1:navg
28+
bfields = [bfield(N, Δt, J, noise),
29+
bfield(N, Δt, J, noise),
30+
bfield(N, Δt, J, noise)];
31+
sol = diffeqsolver([0.8, 0., -0.6], tspan, J, bfields, matrix);
32+
# s[k,:] = mean(sol[3][end-(length(sol[2])÷4):end,:], dims=1)
33+
s[k,:] = mean([sol[1](t)[n] for t in (3*N÷4)*Δt:Δt:N*Δt, n in 1:3], dims=1)
34+
end
35+
Sss[n,:] = mean(s, dims=1)
36+
end
37+
38+
npzwrite("run.npz", Dict("T" => T, "S" => Sss))

‎src/run.npz

738 Bytes
Binary file not shown.

0 commit comments

Comments
 (0)