Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

A Distrubution of cosmic rays

  1. Apr 29, 2017 #1
    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. jcsd
  3. Apr 29, 2017 #2
    Can you please post the code you have written??
     
  4. Apr 29, 2017 #3
    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




     
     
  5. Apr 30, 2017 #4
     
  6. Apr 30, 2017 #5

    stefan r

    User Avatar
    Gold Member

    May be the polar cusp.
     
  7. Apr 30, 2017 #6
    Are polar cusps formed in simple dipole approximation? I didn't classify regions like magnetosheath or magnetopause in my simulation.
     
  8. Apr 30, 2017 #7

    stefan r

    User Avatar
    Gold Member

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Distrubution of cosmic rays
  1. Cosmic Rays (Replies: 1)

  2. Cosmic rays (Replies: 2)

  3. Cosmic ray data (Replies: 3)

Loading...