# A Distrubution of cosmic rays

1. Apr 29, 2017

### RyuuJin

Hey guys!
I am trying to simulate distribution of cosmic ray particles, which change their trajectory due to planet's magnetic field (no atmosphere; dipole approximation). I tried considering protons with non-relativistic velocities, falling on an Earth-like planet, which has rotational axis aligned with the dipole axis. I was expecting a distribution maximum at the poles, but instead I got a result with maxima located at the latitudes of roughly 70 degrees (north and south). Do any of you know of an explanation for this occurrence? Or is just something wrong with my code?
Thanks!

2. Apr 29, 2017

Can you please post the code you have written??

3. Apr 29, 2017

### RyuuJin

Here's the code (Python). I added comments for easier interpretation. I honestly have no idea, what could I be doing wrong...

Code (Python):
# -*- 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 occured

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

4. Apr 30, 2017

### RyuuJin

5. Apr 30, 2017

### stefan r

May be the polar cusp.

6. Apr 30, 2017

### RyuuJin

Are polar cusps formed in simple dipole approximation? I didn't classify regions like magnetosheath or magnetopause in my simulation.

7. Apr 30, 2017

### stefan r

The magnetic field lines stretch further outward at higher latitude. Your density would stop increasing because the magnetic field is no longer deflecting. At 0 degrees (equator) the proton feels a force at 90 degrees. A proton at high latitudes does not feel much force so it continues closer to forward. If a proton is coming in from above northern Siberia moving toward Canada there is not much force to change the direction. It would cross the north pole.

There is also more kilometers at 70 degrees relative to 80 degrees. The 90 degree latitude line is 0 kilometers long. Is your "maximum" a count of protons hitting between degree x and degree y? If so you should expect 0 impacts in the 90+ category. Should get increasing impacts with increasing area until the magnetic defection becomes strong and the map geometry effect becomes weak.