Skip to content
Permalink
Browse files

Updates 4-7-21

  • Loading branch information
nolte committed Apr 7, 2021
1 parent 44fa415 commit a273f632291fa2cfb34600f4ca2473643686e262
Showing with 487 additions and 4 deletions.
  1. +11 −3 Kuramoto.py
  2. BIN Medio.png
  3. +1 −1 NetDynamics.py
  4. +116 −0 caustic.py
  5. +134 −0 gravcaustic.py
  6. +74 −0 gravfront.py
  7. +151 −0 raycaustic.py
@@ -27,6 +27,17 @@
N = 50 # 50
width = 0.2

# model_case 1 = complete graph (Kuramoto transition)
# model_case 2 = Erdos-Renyi
model_case = int(input('Input Model Case (1-2)'))
if model_case == 1:
facoef = 0.2
nodecouple = nx.complete_graph(N)
elif model_case == 2:
facoef = 5
nodecouple = nx.erdos_renyi_graph(N,0.1)


# function: omegout, yout = coupleN(G)
def coupleN(G):

@@ -110,8 +121,6 @@ def flow_deriv(y,t0):
omega = omegatemp - meanomega
sto = np.std(omega)

#nodecouple = nx.complete_graph(N)
nodecouple = nx.erdos_renyi_graph(N,0.1)

lnk = np.zeros(shape = (N,), dtype=int)
for loop in range(N):
@@ -128,7 +137,6 @@ def flow_deriv(y,t0):
xx = np.zeros(shape=(Nfac,))
for facloop in range(Nfac):
print(facloop)
facoef = 0.2

fac = facoef*(16*facloop/(Nfac))*(1/(N-1))*sto/mnomega
for nodeloop in range(N):
BIN +101 KB Medio.png
Binary file not shown.
@@ -26,7 +26,7 @@


Nfac = 25 # 25
N = 36 # 50
N = 100 # 50
width = 0.2

# model_case 1 = Complete Graph
@@ -0,0 +1,116 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Feb 16 19:50:54 2021
caustic.py
@author: nolte
D. D. Nolte, Optical Interferometry for Biology and Medicine (Springer,2011)
"""

import numpy as np
from matplotlib import pyplot as plt
from numpy import random as rnd
from scipy import signal as signal

plt.close('all')

N = 256

def gauss2(sy,sx,wy,wx):

x = np.arange(-sx/2,sy/2,1)
y = np.arange(-sy/2,sy/2,1)
y = y[..., None]

ex = np.ones(shape=(sy,1))
x2 = np.kron(ex,x**2/(2*wx**2));

ey = np.ones(shape=(1,sx));
y2 = np.kron(y**2/(2*wy**2),ey);

rad2 = (x2+y2);

A = np.exp(-rad2);

return A

def cauchy(sy,sx,wy,wx):

x = np.arange(-sx/2,sy/2,1)
y = np.arange(-sy/2,sy/2,1)
y = y[..., None]

ex = np.ones(shape=(sy,1))
x2 = np.kron(ex,x**2/(2*wx**2))

ey = np.ones(shape=(1,sx))
y2 = np.kron(y**2/(2*wy**2),ey)

rad2 = (x2+y2)

A = 1/(1+rad2)

return A

def speckle2(sy,sx,wy,wx):

Btemp = 2*np.pi*rnd.rand(sy,sx);
B = np.exp(complex(0,1)*Btemp);

C = gauss2(sy,sx,wy,wx);

Atemp = signal.convolve2d(B,C,'same');

Intens = np.mean(np.mean(np.abs(Atemp)**2));

D = np.real(Atemp/np.sqrt(Intens));

Dphs = np.arctan2(np.imag(D),np.real(D));

return D, Dphs


Sp, Sphs = speckle2(N,N,N/16,N/16)

#Sp = cauchy(N,N,N/16,N/16)

plt.figure(2)
plt.matshow(Sp,2,cmap=plt.cm.get_cmap('seismic')) # hsv, seismic, bwr
plt.show()

fx, fy = np.gradient(Sp);

fxx,fxy = np.gradient(fx);
fyx,fyy = np.gradient(fy);

J = fxx*fyy - fxy*fyx;

D = np.abs(1/J)

plt.figure(3)
plt.matshow(D,3,cmap=plt.cm.get_cmap('gray')) # hsv, seismic, bwr
plt.clim(0,0.5e7)
plt.show()

eps = 1e-7
cnt = 0
E = np.zeros(shape=(N,N))
for yloop in range(0,N-1):
for xloop in range(0,N-1):

d = N/2

indx = int(N/2 + (d*(fx[yloop,xloop])+(xloop-N/2)/2))
indy = int(N/2 + (d*(fy[yloop,xloop])+(yloop-N/2)/2))

if ((indx > 0) and (indx < N)) and ((indy > 0) and (indy < N)):
E[indy,indx] = E[indy,indx] + 1

plt.figure(4)
plt.imshow(E,interpolation='bicubic',cmap=plt.cm.get_cmap('gray'))
plt.clim(0,20)
plt.xlim(N/4, 3*N/4)
plt.ylim(N/4,3*N/4)
@@ -0,0 +1,134 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Feb 16 19:50:54 2021
caustic.py
@author: nolte
D. D. Nolte, Optical Interferometry for Biology and Medicine (Springer,2011)
"""

import numpy as np
from matplotlib import pyplot as plt
from numpy import random as rnd
from scipy import signal as signal

plt.close('all')

N = 256

def grav(sy,sx,w):

x = np.arange(-sx/2,sy/2,1)
y = np.arange(-sy/2,sy/2,1)
y = y[..., None]

ex = np.ones(shape=(sy,1))
x2 = np.kron(ex,x**2);

ey = np.ones(shape=(1,sx));
y2 = np.kron(y**2,ey);

rad2 = (x2+y2);

A = np.positive(1.0/np.sqrt((1.0 - w/np.sqrt(rad2)))**2 + 1e-10) ;

return A

def gauss2(sy,sx,wy,wx):

x = np.arange(-sx/2,sy/2,1)
y = np.arange(-sy/2,sy/2,1)
y = y[..., None]

ex = np.ones(shape=(sy,1))
x2 = np.kron(ex,x**2/(2*wx**2));

ey = np.ones(shape=(1,sx));
y2 = np.kron(y**2/(2*wy**2),ey);

rad2 = (x2+y2);

A = np.exp(-rad2);

return A

def cauchy(sy,sx,wy,wx):

x = np.arange(-sx/2,sy/2,1)
y = np.arange(-sy/2,sy/2,1)
y = y[..., None]

ex = np.ones(shape=(sy,1))
x2 = np.kron(ex,x**2/(2*wx**2))

ey = np.ones(shape=(1,sx))
y2 = np.kron(y**2/(2*wy**2),ey)

rad2 = (x2+y2)

A = 1/(1+rad2)

return A

def speckle2(sy,sx,wy,wx):

Btemp = 2*np.pi*rnd.rand(sy,sx);
B = np.exp(complex(0,1)*Btemp);

C = gauss2(sy,sx,wy,wx);

Atemp = signal.convolve2d(B,C,'same');

Intens = np.mean(np.mean(np.abs(Atemp)**2));

D = np.real(Atemp/np.sqrt(Intens));

Dphs = np.arctan2(np.imag(D),np.real(D));

return D, Dphs


#Sp = grav(N,N,N/8)

Sp = cauchy(N,N,N/8,N/8)

plt.figure(2)
plt.matshow(Sp,2,cmap=plt.cm.get_cmap('seismic')) # hsv, seismic, bwr
plt.show()

fx, fy = np.gradient(Sp);

fxx,fxy = np.gradient(fx);
fyx,fyy = np.gradient(fy);

J = fxx*fyy - fxy*fyx;

D = np.abs(1/J)

plt.figure(3)
plt.matshow(D,3,cmap=plt.cm.get_cmap('gray')) # hsv, seismic, bwr
plt.clim(0,0.5e7)
plt.show()

eps = 1e-7
cnt = 0
E = np.zeros(shape=(N,N))
for yloop in range(0,N-1):
for xloop in range(0,N-1):

d = 5*N/2

indx = int(N/2 + (d*(fx[yloop,xloop])+(xloop-N/2)/2))
indy = int(N/2 + (d*(fy[yloop,xloop])+(yloop-N/2)/2))

if ((indx > 0) and (indx < N)) and ((indy > 0) and (indy < N)):
E[indy,indx] = E[indy,indx] + 1

plt.figure(4)
plt.imshow(E,interpolation='bicubic',cmap=plt.cm.get_cmap('gray'))
plt.clim(0,20)
plt.xlim(N/4, 3*N/4)
plt.ylim(N/4,3*N/4)
@@ -0,0 +1,74 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Mar 30 19:47:31 2021
gravfront.py
@author: David Nolte
Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019)
Gravitational Lensing
"""

import numpy as np
from matplotlib import pyplot as plt

plt.close('all')

def refindex(x):
n = n0/(1 + abs(x)**expon)**(1/expon);
return n


delt = 0.001
Ly = 10
Lx = 5
n0 = 1
expon = 2 # adjust this from 1 to 10


delx = 0.01
rng = np.int(Lx/delx)
x = delx*np.linspace(-rng,rng)

n = refindex(x)

dndx = np.diff(n)/np.diff(x)

plt.figure(1)
lines = plt.plot(x,n)

plt.figure(2)
lines2 = plt.plot(dndx)

plt.figure(3)
plt.xlim(-Lx, Lx)
plt.ylim(-Ly, 2)
Nloop = 160;
xd = np.zeros((Nloop,3))
yd = np.zeros((Nloop,3))
for loop in range(0,Nloop):
xp = -Lx + 2*Lx*(loop/Nloop)
plt.plot([xp, xp],[2, 0],'b',linewidth = 0.25)

thet = (refindex(xp+delt) - refindex(xp-delt))/(2*delt)
xb = xp + np.tan(thet)*Ly
plt.plot([xp, xb],[0, -Ly],'b',linewidth = 0.25)

for sloop in range(0,3):
delay = n0/(1 + abs(xp)**expon)**(1/expon) - n0
dis = 0.75*(sloop+1)**2 - delay
xfront = xp + np.sin(thet)*dis
yfront = -dis*np.cos(thet)

xd[loop,sloop] = xfront
yd[loop,sloop] = yfront

for sloop in range(0,3):
plt.plot(xd[:,sloop],yd[:,sloop],'r',linewidth = 0.5)





0 comments on commit a273f63

Please sign in to comment.
You can’t perform that action at this time.