# -*- coding: utf-8 -*-
from scipy import *
from scipy.integrate import ode, odeint
from math import *
from pylab import *
from mpl_toolkits.mplot3d import Axes3D
def dipole(x,y,z) :#cartesian dipole magnetic field
xyz=[x,y,z]
len_r = vec_len(xyz)
dir_r = vec_dir(xyz)
mi0=1e-7
A = mi0/(len_r**3)*m_velikost
return A * (3*((np.dot(m_smer,dir_r))*dir_r) - m_smer)
#direction and length of a vector
def vec_len(x) :
return np.sqrt(x[0]**2+x[1]**2+x[2]**2)
def vec_dir(x) :
return np.array([x[0]/vec_len(x), x[1]/vec_len(x), x[2]/vec_len(x)])
'''initial velocity -> from energy'''
def velocity_rel(E,mass,c):
return np.sqrt(c**2-mass**2*c**6/E**2)
def velocity_class(E,masa):
return np.sqrt(2*E/masa)
'''spherical - cartesian'''
def cart_to_spher(x,y,z):
x += 10**(-10)
y += 10**(-10)
z += 10**(-10)
r = np.sqrt(x**2+y**2+z**2)
fi = (np.arctan(y/x))
fi = (fi - 0.5*(np.abs(x)/x - 1) * np.pi) % (2*np.pi)#if x<0
theta = (np.arctan(np.sqrt(x**2+y**2)/z))
theta = (theta - 0.5*(np.abs(z)/z - 1) * np.pi)#if z<0
return r,fi,theta
def spher_to_cart(r,fi,theta):
return [r*np.cos(fi)*np.sin(theta), r*np.sin(fi)*np.sin(theta), r*np.cos(theta)]
def vec_len(x) :
return np.sqrt(x[0]**2+x[1]**2+x[2]**2)
'''gives trajectory of a particle'''
def SolveNewtonLorenz_2(x0, y0, z0, vx0, vy0, vz0, tend, dt, m, q,field=dipole):
vsq = vx0**2 + vy0**2 + vz0**2
gamma = 1.0/sqrt(1 - vsq/c**2)
def f(Y, t=0):
x,y,z,vx,vy,vz = Y
Bx, By, Bz = field(x,y,z)
fac = q/(m*gamma)
return [ vx, vy, vz,
fac*(vy*Bz - vz*By),
fac*(vz*Bx - vx*Bz),
fac*(vx*By - vy*Bx) ]
t = arange(0, tend, dt)
y = odeint(f, [x0, y0, z0, vx0, vy0, vz0], t)
coll_mat=[]
for i in range(len(y)):
vec=[y[i,0],y[i,1],y[i,2]]
if vec_len(vec)<=Re:
print("Aaay papi! It's a collision!")
coll_mat=np.append(coll_mat,vec)
break
return y, coll_mat'''gives distribution of particles on the surface (theta, phi)'''
def distrib(N,r0, E0, e, mass, tend, dt):
theta_tab=[]
fi_tab=[]
hits=0
#I generate particles at distance r0 (theta, phi are random), with certain direction
for i in range(N):
#Gaussian velocity distribution
E=normal(E0,scale=E0/10)
#incoming direction
x=rand(1)
theta=acos(2*x-1)
y=rand(1)
fi=(2*np.pi*y)
#position of particle
fi0=uniform(0,2*np.pi)
x1=rand(1)
theta0=acos(2*x1-1)
v0=velocity_class(E,mass)
v_cart=spher_to_cart(v0,fi,theta)
#if velocity too large -> general relativity
if v0 >=0.1*c:
v0=hitrost(E,masa,c)
v_kart=spher_to_cart(v0,fi,theta)
vx=v_kart[0]
vy=v_kart[1]
vz=v_kart[2]
#I adjust direction of particles (all are flying towards or next to the planet - not from it!)
if theta0>=0 and theta0<np.pi/2 and fi0 >=0 and fi0<np.pi/2:#1.octant
vx=-np.abs(vx)
vy=-np.abs(vy)
vz=-np.abs(vz)
elif theta0>=0 and theta0<np.pi/2 and fi0 >np.pi/2 and fi0<np.pi:#2. octant
vx=np.abs(vx)
vy=-np.abs(vy)
vz=-np.abs(vz)
elif theta0>=0 and theta0<np.pi/2 and fi0 >np.pi and fi0<np.pi*3/2:#3. octant
vx=np.abs(vx)
vy=np.abs(vy)
vz=-np.abs(vz)
elif theta0>=0 and theta0<np.pi/2 and fi0 >np.pi*3/2 and fi0<2*np.pi:#4. octant
vx=-np.abs(vx)
vy=np.abs(vy)
vz=-np.abs(vz)
elif theta0>=np.pi/2 and theta0<np.pi and fi0 >=0 and fi0<np.pi/2:#5. octant
vx=-np.abs(vx)
vy=-np.abs(vy)
vz=np.abs(vz)
elif theta0>=np.pi/2 and theta0<np.pi and fi0 >np.pi/2 and fi0<np.pi:#6. octant
vx=np.abs(vx)
vy=-np.abs(vy)
vz=np.abs(vz)
elif theta0>=np.pi/2 and theta0<np.pi and fi0 >np.pi and fi0<np.pi*3/2:#7. octant
vx=np.abs(vx)
vy=np.abs(vy)
vz=np.abs(vz)
elif theta0>=np.pi/2 and theta0<np.pi and fi0 >np.pi*3/2 and fi0<2*np.pi:#8. octant
vx=-np.abs(vx)
vy=+np.abs(vy)
vz=np.abs(vz)
elif theta0==0:
vz=-np.abs(vz)
if theta0==np.pi:
vz=np.abs(vz)
if theta0==np.pi/2 and fi0==0 or fi0==2*np.pi:
vx=-np.abs(vx)
if theta0==np.pi/2 and fi0==np.pi:
vx=np.abs(vx)
if theta0==np.pi/2 and fi0==np.pi/2:
vy=-np.abs(vx)
if theta0==np.pi/2 and fi0==2*np.pi:
vy=np.abs(vy)
startVel=[float(vx),float(vy),float(vz)]
startPos=spher_to_cart(r0,fi0,theta0)
posMat,cPos= SolveNewtonLorenz_2(startPos[0],startPos[1],startPos[2], startVel[0],startVel[1],startVel[2], tend,dt, mass,e, field=dipole)
if len(cPos)>1:#if collision occurred
r_coll,fi_coll,theta_coll=cart_to_spher(cPos[0],cPos[1],cPos[2])
theta_tab=np.append(theta_tab,theta_coll)
fi_tab=np.append(fi_tab,fi_coll)
hits+=1
print(hits)
return fi_tab,theta_tab, posMat,cPos