
Notifications
You must be signed in to change notification settings  Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
 Loading branch information
Showing
1 changed file
with
132 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number  Diff line number  Diff line change 

@@ 0,0 +1,132 @@  
#!/usr/bin/env python3  
# * coding: utf8 *  
"""  
SteinsGate2D.py  
Created on Sat March 6, 2021  
@author: David Nolte  
Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019)  
2D simulation of Stein's Gate Divergence Meter  
"""  
import numpy as np  
from scipy import integrate  
from matplotlib import pyplot as plt  


plt.close('all')  


def solve_flow(param,lim = [6,6,6,6],max_time=20.0):  




def flow_deriv(x_y, t0, alpha, beta, gamma):  
#"""Compute the timederivative ."""  
x, y = x_y  


w = 1  
R2 = x**2 + y**2  
R = np.sqrt(R2)  
arg = (R2)/0.1  
env1 = 1/(1+np.exp(arg))  
env2 = 1  env1  


f = env2*(x*(1/(R1.99)**2 + 1e2)  x) + env1*(w*y + w*x*(1  R))  
g = env2*(y*(1/(R1.99)**2 + 1e2) + y) + env1*(w*x + w*y*(1  R))  


return [f,g]  
model_title = 'Steins Gate'  


plt.figure()  
xmin = lim[0]  
xmax = lim[1]  
ymin = lim[2]  
ymax = lim[3]  
plt.axis([xmin, xmax, ymin, ymax])  


N = 24*4 + 47  
x0 = np.zeros(shape=(N,2))  
ind = 1  
for i in range(0,24):  
ind = ind + 1  
x0[ind,0] = xmin + (xmaxxmin)*i/23  
x0[ind,1] = ymin  
ind = ind + 1  
x0[ind,0] = xmin + (xmaxxmin)*i/23  
x0[ind,1] = ymax  
ind = ind + 1  
x0[ind,0] = xmin  
x0[ind,1] = ymin + (ymaxymin)*i/23  
ind = ind + 1  
x0[ind,0] = xmax  
x0[ind,1] = ymin + (ymaxymin)*i/23  
ind = ind + 1  
x0[ind,0] = 0.05  
x0[ind,1] = 0.05  


for thetloop in range(0,10):  
ind = ind + 1  
theta = 2*np.pi*(thetloop)/10  
ys = 0.125*np.sin(theta)  
xs = 0.125*np.cos(theta)  
x0[ind,0] = xs  
x0[ind,1] = ys  


for thetloop in range(0,10):  
ind = ind + 1  
theta = 2*np.pi*(thetloop)/10  
ys = 1.7*np.sin(theta)  
xs = 1.7*np.cos(theta)  
x0[ind,0] = xs  
x0[ind,1] = ys  


for thetloop in range(0,20):  
ind = ind + 1  
theta = 2*np.pi*(thetloop)/20  
ys = 2*np.sin(theta)  
xs = 2*np.cos(theta)  
x0[ind,0] = xs  
x0[ind,1] = ys  


ind = ind + 1  
x0[ind,0] = 3  
x0[ind,1] = 0.05  
ind = ind + 1  
x0[ind,0] = 3  
x0[ind,1] = 0.05  
ind = ind + 1  
x0[ind,0] = 3  
x0[ind,1] = 0.05  
ind = ind + 1  
x0[ind,0] = 3  
x0[ind,1] = 0.05  
ind = ind + 1  
x0[ind,0] = 6  
x0[ind,1] = 0.00  
ind = ind + 1  
x0[ind,0] = 6  
x0[ind,1] = 0.00  


colors = plt.cm.prism(np.linspace(0, 1, N))  


# Solve for the trajectories  
t = np.linspace(0, max_time, int(250*max_time))  
x_t = np.asarray([integrate.odeint(flow_deriv, x0i, t, param)  
for x0i in x0])  


for i in range(N):  
x, y = x_t[i,:,:].T  
lines = plt.plot(x, y, '', c=colors[i])  
plt.setp(lines, linewidth=1)  




plt.show()  
plt.title(model_title)  




return t, x_t  


param = (0.02,0.5,0.2) # Steins Gate  
lim = (6,6,6,6)  


t, x_t = solve_flow(param,lim)  


plt.savefig('Steins Gate') 