In [1]:
import numpy as np
In [2]:
from matplotlib import pyplot as plt
In [3]:
%matplotlib inline
In [4]:
# %load Multiple_Coll_01.py
import numpy as np
from matplotlib import pyplot as plt
def orto(x):
if np.dot(x,x) == 0:
return 'No se puede: ese es el vector cero!'
else:
if 0 not in x:
v1 = 1
v2 = -(x[0]/x[1])
v3 = 0
#return np.array([v1,v2,v3])
else:
if x[0] == 0:
if x[1] == 0:
v1 = 1
v2 = 0
v3 = 0
else:
v1 = 0
v2 = 0
v3 = 1
elif x[1] == 0:
v1 = 0
v2 = 1
v3 = 0
else:
v1 = 0
v2 = 0
v3 = 1
return np.array([v1,v2,v3])
#Funcion que regresa dos vectores; numpy arrays de 3D, ortogonales al vector de input x.
#Esto es, devuelven la base al espacio ortogonal definido por el vector x de entrada.
#@jit
def base_ort_nor(x):
y = orto(x)
v1 = y/np.linalg.norm(y)
z = np.cross(x,v1)
v2 = z/np.linalg.norm(z)
return v1, v2
#Esta funcion genera un vector con distrubucion uniforme en las direcciones sobre un plano tangente a la esfera de radio R.
#@jit
def vector_des(v1,v2):
na = 2*np.pi*np.random.rand()
vn = v1*np.cos(na) + v2*np.sin(na)
return vn/np.linalg.norm(vn)
R = 1
#Normalizamos al vector de desplazamiento para que intersecte al vector de la nueva posicion de acuerdo con que el
#desplazamiento (s) sobre la esfera, sobre este arco de circulo maximo, sea el determinado por el movimiento browniano particular.
#@jit
def vector_q(x,s):
q = (R)*np.tan(s/(R))
return q*x
#Dados todos los datos anteriores, esta funcion actualiza la posicion de la particula.
#Lo que hace es que la mueve sobre el plano tangente a la esfera en la direccion predeterminada de tal suerte que el desplazamiento efectivo
#s sobre una geodesica de la esfera, se el deseado, y posteriormente la proyecta sobre la superficie de la esfera.
#@jit
def nuevo_r(r, vector_q):
y = r + vector_q
y = y/np.linalg.norm(y)
return (R)*y
#Esta funcion ensambla todo lo anterior: como imput necesita una posicion inicial y un arco de desplazamiento
#Como output da un vector de posicion nuevo dada un tipo de desplazamiento.
#@jit
def actualiza(r,s):
v1, v2 = base_ort_nor(r)
pre_q = vector_des(v1,v2)
q = vector_q(pre_q, s)
return nuevo_r(r, q)
#Esta funcion actualiza la posicion de todos los elementos de una lista; particula brownianas.
#@jit
def act_n(lista,s):
l = []
for v in lista:
l.append(actualiza(v,s))
return l
#Huella de la trayectoria
#La siguiente funcion hace una particion de la trayectoria sobre s en n pedazos y regresa
#una lista de los vectores de esas posiciones sobre la esfera.
#Usa al operador de rotacion.
#@jit
def b_steps_(ri,rf,n):
l = [ri]
r0 = ri
lamb = (np.dot(ri,rf))/((np.linalg.norm(ri))*(np.linalg.norm(rf)))
if abs(lamb) > 1:
#print 'Is fucked up: there was a rounding '
if lamb < 0:
lamb = -1
else:
lamb = 1
theta = np.arccos(lamb)
#if theta < 1e17:
#return l
if theta == 0:
return [ri,rf]
else:
normal = np.cross(ri, rf)/ np.linalg.norm(np.cross(ri,rf))
for i in range(1,n + 1):
#vi = rot_theta(r0, theta/n, normal)
vi = rot_finita(r0, -normal, theta/n)
l.append(vi)
r0 = vi
return l
#Operador de Rotacion
#Depende de los parametros r, el vector o punto que queremos rotar; theta el angulo de rotacion; n el vector que define el eje de rotacion y el signo de rotacion.
#@jit
def rot_theta(r, theta, u):
x = np.array([np.cos(theta) + (u[0]*u[0])*(1 - np.cos(theta)), u[0]*u[1]*(1 - np.cos(theta)) - u[2]*np.sin(theta), u[0]*u[2]*(1 - np.cos(theta)) + u[1]*np.sin(theta)])
y = np.array([u[1]*u[0]*(1 - np.cos(theta)) + u[2]*np.sin(theta), np.cos(theta) + u[1]*u[1]*(1 - np.cos(theta)), u[1]*u[2]*(1 - np.cos(theta)) - u[0]*np.sin(theta)])
z = np.array([u[2]*u[0]*(1 - np.cos(theta)) - u[1]*np.sin(theta), u[2]*u[1]*(1 - np.cos(theta)) + u[0]*np.sin(theta), np.cos(theta) + u[2]*u[2]*(1 - np.cos(theta))])
R = np.array([x,y,z])
return np.dot(R, r)
#Transformacion de coordenada de esfericas a cartesianas.
#@jit
def trans_s_c(r,theta, phi):
x = r*np.sin(theta)*np.cos(phi)
y = r*np.sin(theta)* np.sin(phi)
z = r*np.cos(theta)
return x, y, z
#Transformacion de coordenadas de cartesianas a esfericas.
#@jit
def trans_c_s(x,y,z):
r = np.sqrt(x**2 + y**2 + z**2)
#print r
cociente = z/r
if abs(cociente) > 1:
if cociente < 0:
theta = np.arccos(-1)
else:
theta = np.arccos(1)
else:
theta = np.arccos(z/r)
phi = np.arctan(y/x)
return r, theta, phi
#@jit
def r_uni(theta, phi):
x = np.sin(theta)*np.cos(phi)
y = np.cos(theta)*np.cos(phi)
z = np.cos(theta)
return np.array([x,y,z])
#@jit
def theta_uni(theta, phi):
x = np.cos(theta)*np.cos(phi)
y = np.cos(theta)*np.sin(phi)
z = -np.sin(theta)
return np.array([x,y,z])
#@jit
def phi_uni(theta, phi):
x = -np.sin(phi)
y = np.cos(phi)
z = 0
return np.array([x,y,z])
#@jit
def nombre(s):
diferencia = 4 - len(str(s))
ceros = ''
for i in range(diferencia):
ceros = ceros + '0'
variable = ceros + str(s)
return variable
#Varianza para una distribucion bigaussiana; difusion en 2D
#@jit
def var(D, delta_t):
return 4*D*delta_t
#Arco de circulo maximo con distribucion normal alrededor de cero y una varianza dada por
#@jit
def ese(D,delta_t):
return abs(np.random.normal(loc = 0., scale = np.sqrt(var(D,delta_t)),size = None))
#Funcion de rotacion finita
#@jit
def rot_finita(r_ini, N, Phi):
n = N/np.linalg.norm(N)
r_fin = np.cos(Phi)*r_ini + (np.dot(n,r_ini))*(1 - np.cos(Phi))*n + (np.sin(Phi))*(np.cross(r_ini,n))
return r_fin
#Funcion que regresa una lista de n numpy arrays que son l
def Trayectoria(ri,rf,n):
l = [ri]
r0 = ri
theta = np.arccos((np.dot(ri,rf))/((np.linalg.norm(ri))*(np.linalg.norm(rf))))
N = np.cross(ri, rf)
for i in range(1,n + 1):
vi = rot_finita(r0, N, theta/n)
l.append(vi)
r0 = vi
return l
#Collision_check es una función que, dada una trayectoria: una lista de vectores que
#pasan por puntos sucesivos de la trayectoria, verifica si alguna de estas posiciones
#interesecto a alguno de los obstáculos. En caso de que así sea, actualiza conforme una
#colision elastica. En caso de no intersectar a ningun obstaculo regresa una lista
#con dos vectores: posicion inicial y posicion final en ese orden.
#@jit
def penetrate_obs(lista_vect, lista_obs, size):
metiches = []
for obs in lista_obs:
theta_omega = size
r_omega = obs
frontera = .2
#metiches = []
for v in lista_vect:
tamanho = np.cos(theta_omega - frontera)
if np.dot(v,r_omega) > tamanho:
print 'Penetro el mother fucker obstacle'
metiches.append(v)
else:
continue
#print 'no choco el mother fucker'
#valor = False
return metiches
#Esta funcion cuando es llamada grafia la posicion de las partoculas brownianas.
#sobre la superficie de una esfera sobre la que se esta difundiendo.
#@jit
def plot_particles(lista, vpolar, vazim, numero):
from mpl_toolkits.mplot3d import axes3d
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
#import matplotlib.pyplot as plt
#import numpy as np
from itertools import product, combinations
fig = plt.figure(figsize=(20,10))
ax = fig.gca(projection='3d')
ax.set_aspect("equal")
ax._axis3don = False
#draw sphere
R = 1
u, v = np.mgrid[0:2*np.pi:50j, 0:np.pi:50j]
x=R*np.cos(u)*np.sin(v)
y=R*np.sin(u)*np.sin(v)
z=R*np.cos(v)
#ax.plot_surface(x, y, z, color="r", alpha = 0.15)
ax.plot_surface(x, y, z, cmap=cm.YlGnBu_r,rstride=1, cstride=1, alpha = 0.10, linewidth = 0.15)
ax.view_init(vpolar, vazim)
#draw an arrow or a set of arrow
ax.quiver(0,0,1.5,0,0,1, length=0.5, arrow_length_ratio = .5)
#draw patch
#u, v = np.mgrid[0:2*np.pi:50j, 0:(np.pi/7):50j]
#x=R*np.cos(u)*np.sin(v)
#y=R*np.sin(u)*np.sin(v)
#z=R*np.cos(v)
#ax.plot_surface(x, y, z, color="r", alpha = 0.25)
#draw points
for p in lista:
ax.scatter([p[0]],[p[1]],[p[2]],color="b",s=50, alpha = 0.6)
#fig.savefig('Obs_Dist_Images{}.png'.format(nombre(numero)))
#ax.view_init(80, 30)
plt.show()
#plt.close()
#@jit
def polo_n(n, R):
l = []
for i in range(n):
l.append(np.array([0,0,R]))
return l
#@jit
def particion_esfera(ccero, Nphi):
Ntheta = int(4*np.pi/(ccero*Nphi))
print 'Ntheta', Ntheta, 'Nphi', Nphi, 'Ntheta*Nphi', Ntheta*Nphi
sigmaPhi = 2*np.pi/Nphi
deltaphi = 2*np.pi/Nphi
thetas = []
phis = [0]
cociente = ccero/sigmaPhi
for i in range(Ntheta + 1):
theta = np.arccos(1 - (i)*cociente)
thetas.append(theta)
for j in range(Nphi):
phis.append(phis[j] + deltaphi)
return thetas, phis
#@jit
def secuencia_part(tamini, Nfi, numero):
l1, l2 = particion_esfera(4*np.pi/tamini, Nfi)
particion = []
for i in range(len(l2)):
for j in range(len(l1)):
x, y, z = trans_s_c(1, l1[j], l2[i])
particion.append(np.array([x, y, z]))
return plot_particles(particion, 45, 45, numero)
#Funcion que regresa las coordenadas del centro de dos arreglos para
#las coordenadas theta y phi
#@jit
def coordenadas_centro(l1,l2):
thetas_centro = []
phis_centro = []
for i in range(len(l1) - 1):
theta_media = l1[i] + (l1[i + 1] - l1[i])/2.
thetas_centro.append(theta_media)
for j in range(len(l2) - 1):
phi_media = l2[j] + (l2[j + 1] - l2[j])/2.
phis_centro.append(phi_media)
return thetas_centro, phis_centro
#@jit(nopython=True)
def secuencia_obs(N, Nfi, numero):
l1_prima, l2_prima = particion_esfera(4*np.pi/N, Nfi)
l1, l2 = coordenadas_centro(l1_prima, l2_prima)
particion = []
for i in range(len(l2)):
for j in range(len(l1)):
x, y, z = trans_s_c(1, l1[j], l2[i])
particion.append(np.array([x, y, z]))
print len(particion)
#return plot_particles(particion, 0, 0, numero)
return particion
def plot_particle_traj_obs(lista_obstaculos, trayectoria, vpolar, vazim, numero):
from mpl_toolkits.mplot3d import axes3d
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
#import matplotlib.pyplot as plt
#import numpy as np
from itertools import product, combinations
fig = plt.figure(figsize=(20,10))
ax = fig.gca(projection='3d')
ax.set_aspect("equal")
ax._axis3don = False
#draw sphere
R = 1
u, v = np.mgrid[0:2*np.pi:50j, 0:np.pi:50j]
x=R*np.cos(u)*np.sin(v)
y=R*np.sin(u)*np.sin(v)
z=R*np.cos(v)
#ax.plot_surface(x, y, z, color="r", alpha = 0.15)
ax.plot_surface(x, y, z, cmap=cm.YlGnBu_r,rstride=1, cstride=1, alpha = 0.10, linewidth = 0.15)
ax.view_init(vpolar, vazim)
#draw patch
#u, v = np.mgrid[0:2*np.pi:50j, 0:(np.pi/7):50j]
#x=R*np.cos(u)*np.sin(v)
#y=R*np.sin(u)*np.sin(v)
#z=R*np.cos(v)
#ax.plot_surface(x, y, z, color="r", alpha = 0.25)
#draw obstacles
for p in lista_obstaculos:
ax.scatter([p[0]],[p[1]],[p[2]], color="b", s=10, alpha = 0.2)
#draw trajectory
for p in trayectoria:
ax.scatter([p[0]],[p[1]],[p[2]], color="k",s=20, alpha = 0.7)
#Plot the x positive direction
ax.quiver(1.5,0,0,1,0,0, length=0.5, arrow_length_ratio = .5)
ax.quiver(0,1.5,0,0,1,0, length=0.5, arrow_length_ratio = .5)
ax.quiver(0,0,1.5,0,0,1, length=0.5, arrow_length_ratio = .5)
#fig.savefig('BS_24_Obs_test_01{}.png'.format(nombre(numero + 1)))
#ax.view_init(80, 30)
#plt.close()
plt.show()
def obs_uniforme(N, R, size):
list_obs = []
omega = np.cos(size)
while len(list_obs) < N:
x, y, z = np.random.uniform(-1,1), np.random.uniform(-1,1), np.random.uniform(-1,1)
v = np.array([x, y, z])
norma = np.linalg.norm(v)
if norma <= R:
n = v/norma
if not np.dot(n, np.array([0.,0.,1.]))/R > omega:
list_obs.append(R*n)
return list_obs
def puntos_obs_j(r_omega, theta_omega, n):
r , theta, phi = trans_c_s(r_omega[0],r_omega[1],r_omega[2])
rp = rot_finita(r_omega, phi_uni(theta, phi), theta_omega)
puntos_obs_j = [rp]
for i in range(1,n):
x = rot_finita(rp, r_omega, 2*np.pi/n)
puntos_obs_j.append(x)
rp = x
return puntos_obs_j
def puntos_obs(lista_obstaculos, size):
mis_obs = []
for i in range(len(lista_obstaculos)):
a = lista_obstaculos[i]
b = size
mis_obs = mis_obs + puntos_obs_j(a,b,100)
return mis_obs
#Collision_check es una función que, dada una trayectoria: una lista de vectores que
#pasan por puntos sucesivos de la trayectoria, verifica si alguna de estas posiciones
#interesecto a alguno de los obstáculos. En caso de que así sea, actualiza conforme una
#colision elastica. En caso de no intersectar a ningun obstaculo regresa una lista
#con dos vectores: posicion inicial y posicion final en ese orden.
def penetrate_obs(lista_vect, lista_obs, size):
metiches = []
for obs in lista_obs:
r_omega, theta_omega = obs, size
#frontera = .2
#metiches = []
for v in lista_vect:
tamanho = np.cos(theta_omega)
if np.dot(v,r_omega) >= tamanho:
#print 'Penetro el mother fucker obstacle'
metiches.append(v)
else:
continue
#print 'no choco el mother fucker'
#valor = False
return metiches
def check_collision(lista_vect, lista_obs, size):
for obs in lista_obs:
r_omega, theta_omega = obs, size
for v in lista_vect:
tamanho = np.cos(theta_omega)
if np.dot(v,r_omega) > tamanho:
return True
else:
continue
return False
In [5]:
def tangent_space(x1,x2,xo):
np_prima = np.cross(np.cross(xo, x1), x1)
nor_p = np_prima/np.linalg.norm(np_prima)
up_prima = np.cross(np.cross(x1, x2), x1)
up = up_prima/np.linalg.norm(up_prima)
tp_prima = np.cross(x1, nor_p)
tp = tp_prima/np.linalg.norm(tp_prima)
y = (np.dot(up,tp))*tp - (np.dot(up, nor_p))*nor_p
v_rot_prima = np.cross(x1, y)
v_rot = v_rot_prima/np.linalg.norm(v_rot_prima)
return v_rot
In [6]:
def collision_check(Lv, Lo, size):
if len(Lv) == 1:
return Lv[0], Lv[0], 0
elif len(Lv) == 2:
sproduct = np.dot(Lv[0], Lv[1])
if abs(sproduct) > 1:
if sproduct < 0:
sproduct = -1.
else:
sproduct = 1.
alpha = np.arccos(sproduct)
for obs in Lo:
r_o = obs
if np.dot(Lv[1],r_o) > np.cos(size):
v_rot = tangent_space(Lv[0], Lv[1], r_o)
x_final = rot_finita(Lv[1], -v_rot, alpha)
return Lv[1], x_final, 0
else:
for k in range(1,len(Lv)):
for obs in Lo:
r_o = obs
if np.dot(Lv[k],r_o) > np.cos(size):
v_rot = tangent_space(Lv[0], Lv[k], r_o)
if k == len(Lv) - 1:
sproduct = np.dot(Lv[k-1], Lv[k])
if abs(sproduct) > 1:
if sproduct < 0:
sproduct = -1.
else:
sproduct = 1.
else:
sproduct = np.dot(Lv[k], Lv[-1])
if abs(sproduct) > 1:
if sproduct < 0:
sproduct = -1.
else:
sproduct = 1.
alpha = np.arccos(sproduct)
x_final = rot_finita(Lv[k], -v_rot, alpha)
return Lv[k], x_final, k
return Lv[0], Lv[-1], 0
In [ ]:
In [ ]:
In [7]:
def uniform_inside_omega(N, R, size):
list_obs = []
omega = np.cos(size)
while len(list_obs) < N:
x, y, z = np.random.uniform(-1,1), np.random.uniform(-1,1), np.random.uniform(-1,1)
v = np.array([x, y, z])
norma = np.linalg.norm(v)
if norma <= R:
n = v/norma
if not np.dot(n, np.array([0.,0.,1.]))/R < omega:
list_obs.append(R*n)
return list_obs
In [ ]:
In [49]:
D = 1e-2
dt = (2.3**2/4.)*1e-5
theta1 = 0.4
theta2 = 0.1
In [ ]:
In [50]:
pos_ini = uniform_inside_omega(100,1, theta1)
In [51]:
plot_particles(pos_ini,90,0,1)
In [9]:
In [21]:
Out[21]:
In [ ]:
In [ ]:
In [ ]:
In [60]:
#Simulacion del Primer tiempo de salida para una region de tamaño.
pos_ini = uniform_inside_omega(10000,1, theta1-0.05)
tiempo = 0
tiempos_salida = []
tamanho = np.cos(theta1)
r_omega = np.array([0,0,1])
#plot_particles(pos_ini,90,0,1)
frozen = []
#plot_particle_Region_frozen(estructura_obs, pos_ini, frozen, 90, 0, 1)
for i in range(10000):
tiempo = tiempo + dt
nuevas_pos = act_n(pos_ini, ese(D,dt))
#plot_particles(nuevas_pos, 90,0,1)
#plot_particle_Region_frozen(estructura_obs, nuevas_pos, frozen, 90, 0, 1)
survey = []
for v in nuevas_pos:
if np.dot(r_omega, v) <= tamanho:
tiempos_salida.append(tiempo)
frozen.append(v)
else:
survey.append(v)
if len(survey) == 0:
break
else:
pos_ini = survey