What Causes a Shift in Cosmic Ray Distribution on Earth?

In summary, the conversation discusses the simulation of cosmic ray particles and their trajectory when affected by a planet's magnetic field. The code used for the simulation, written in Python, is shared and explained. The results obtained show a distribution of maximums at latitudes of approximately 70 degrees instead of at the poles as expected. The speaker is seeking an explanation for this occurrence or if there is an error in the code. A code snippet is also provided for easier interpretation.
  • #1
RyuuJin
7
0
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!
 
Astronomy news on Phys.org
  • #2
Can you please post the code you have written??
 
  • #3
Here's the code (Python). I added comments for easier interpretation. I honestly have no idea, what could I be doing wrong...
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
RyuuJin said:
Here's the code (Python). I added comments for easier interpretation. I honestly have no idea, what could I be doing wrong...

EDIT: fixed typing/naming error in line 113

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=velocity(E,mass,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
 
  • #5
May be the polar cusp.
 
  • #6
stefan r said:
May be the polar cusp.
Are polar cusps formed in simple dipole approximation? I didn't classify regions like magnetosheath or magnetopause in my simulation.
 
  • #7
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.
 

What are cosmic rays?

Cosmic rays are high-energy particles that originate from outer space. They can range in size from protons and electrons to larger atomic nuclei and can travel at near the speed of light.

How are cosmic rays distributed throughout the universe?

Cosmic rays are distributed throughout the universe in a relatively uniform manner. However, their distribution can be affected by various factors such as the magnetic fields of stars and galaxies, as well as the presence of supernovae and other cosmic events.

What is the primary source of cosmic rays?

The primary source of cosmic rays is still not fully understood, but it is believed that the majority come from supernova explosions in our galaxy. Other possible sources include active galactic nuclei, gamma ray bursts, and black holes.

What is the energy range of cosmic rays?

Cosmic rays can have a wide range of energies, from very low energy particles to extremely high energy particles. The highest energy cosmic rays can have energies up to 10^20 electron volts (eV), which is millions of times more energetic than the particles produced in the Large Hadron Collider.

How do cosmic rays interact with Earth's atmosphere?

When cosmic rays enter Earth's atmosphere, they collide with atoms and molecules, producing secondary particles such as pions, muons, and electrons. These secondary particles can then collide with other particles, creating a cascade of particles that eventually reach the Earth's surface.

Similar threads

I Are Cosmic Rays Intensifying Across the Globe?
    • Astronomy and Astrophysics
Replies
4
Views
8K
B Where do the electrons go? (in stellar nuclear fusion)
    • Astronomy and Astrophysics
Replies
4
Views
1K
B "Without Earth's atmosphere, cosmic rays would be lethal"
    • Astronomy and Astrophysics
Replies
9
Views
2K
Has a change in the cosmic rays hitting Earth been causing Global Warming?
    • Earth Sciences
Replies
8
Views
1K
I What´s the cosmic ray's energy dependence?
    • High Energy, Nuclear, Particle Physics
Replies
3
Views
1K
B What would it be like on Earth if it were not a sphere?
    • Astronomy and Astrophysics
2
Replies
42
Views
9K
Why aren't all Cosmic Ray energies 'effective'?
    • High Energy, Nuclear, Particle Physics
2
Replies
43
Views
4K
I E- & Ions Deflection by Earth's Magnetic Field
    • Astronomy and Astrophysics
Replies
30
Views
5K
LHC vs. Cosmic Rays: What Role Does Density Play in Particle Research?
    • High Energy, Nuclear, Particle Physics
Replies
25
Views
3K
I Orientation of the Earth, Sun and Solar System in the Milky Way
    • Astronomy and Astrophysics
5
Replies
142
Views
113K

Hot Threads

Recent Insights

Back
Top