Variable changes value without changing it at all in the program

Question

The numpy array pos_0 changes its value without anything happening to it in the program.

Relevant steps:

1. I assign a value to pos_0
2. I set pos=pos_0
3. I change pos (in a while loop)
4. I print both pos and pos_0 and pos_0 is now equal to the value that pos has after the while loop.

Nowhere after the assignment of pos=pos_0 does the name of the variable even come up.

Also, the program takes FOREVER to run. I know it's a long loop so it's not a surprise, but any advice on how to speed it up would be so great.

Thanks a ton!

Here is the code:

import math
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import sys
import random

#starting conditions:
q_b=5.19230065e-39                                                      #the magnetic charge of the particle in Ampere*kpc by dirac quantization cond.
m=1.78266184e-10                                                        #mass (kg). An estimate from Wick (2002).
pos_0=np.array([-8.33937979, 0.00, 0.00])                               #starting position (kpc)
vel_0=np.array([0, 0, 0])*(3.24077929e-20)                       #starting velocity (kpc/s), but enter the numbers in m, conversion is there.
dt=1e8                                                                 #the timestep (smaller = more accuracy, more computing time) in seconds
distance_to_track= .01                                                    #how far to track the particle (kpc)

#disk parameters. B is in tesla.
b1=0.1e-10
b2=3.0e-10
b3=-0.9e-10
b4=-0.8e-10
b5=-2.0e-10
b6=-4.2e-10
b7=0.0e-10
b8=2.7e-10
b_ring=0.1e-10
h_disk=0.4
w_disk=0.27

#halo parameters
B_n=1.4e-10
B_s=-1.1e-10
r_n=9.22
r_s=16.7 #NOTE: parameter has a large error.
w_h=0.2
z_0=5.3

#X-field parameters
B_X=4.6e-10
Theta0_X=49.0
rc_X=4.8
r_X=2.9

#striation parameters (the striated field is currently unused)
#gamma= 2.92
#alpha= 2.65 #taken from page 9, top right paragraph
#beta= gamma/alpha-1

#other preliminary measures:
c=9.7156e-12    #speed of light in kpc/s

def mag(V):
    return math.sqrt(V[0]**2+V[1]**2+V[2]**2)

initialposition=pos_0

pos=pos_0
vel=vel_0

trailx=(pos[0],)                    #trailx,y,z is used to save the coordinates of the particle at each step, to plot the path afterwards
traily=(pos[1],)
trailz=(pos[2],)

gam=1/math.sqrt(1-(mag(vel)/c)**2)

KE=m*c**2*(gam-1)*5.942795e48       #KE, converted to GeV
KEhistory=(KE,)

distance_tracked=0                  #set the distance travelled so far to 0

time=0

#boundary function (between disk and halo fields)
def L(Z,H,W):
    return (1+math.e**(-2*(abs(Z)-H)/W))**-1

#halo boundary spirals:
i=11.5                                                        #this is the opening "pitch" angle of the logarithmic spiral boundaries.
r_negx=np.array([5.1, 6.3, 7.1, 8.3, 9.8, 11.4, 12.7, 15.5])  #these are the radii at which the spirals cross the x-axis

def r1(T):
    return r_negx[0]*math.e**(T/math.tan(math.radians(90-i)))
def r2(T):
    return r_negx[1]*math.e**(T/math.tan(math.radians(90-i)))
def r3(T):
    return r_negx[2]*math.e**(T/math.tan(math.radians(90-i)))
def r4(T):
    return r_negx[3]*math.e**(T/math.tan(math.radians(90-i)))
def r5(T):
    return r_negx[4]*math.e**(T/math.tan(math.radians(90-i)))
def r6(T):
    return r_negx[5]*math.e**(T/math.tan(math.radians(90-i)))
def r7(T):
    return r_negx[6]*math.e**(T/math.tan(math.radians(90-i)))
def r8(T):
    return r_negx[7]*math.e**(T/math.tan(math.radians(90-i)))

#X-field definitions:

def r_p(R,Z):
    if abs(Z)>=math.tan(Theta0_X)*math.sqrt(R)-math.tan(Theta0_X)*rc_X:
        return R-abs(Z)/math.tan(Theta0_X)
    else:
        return R*rc_X/(rc_X+abs(Z)/math.tan(Theta0_X))

def b_X(R_P):
    return B_X*math.e**(-R_P/r_X)

def Theta_X(R,Z):
    if abs(Z)>=math.tan(Theta0_X)*math.sqrt(R)-math.tan(Theta0_X)*rc_X:
        return math.atan(abs(Z)/(R-r_p(R,Z)))
    else:
        return Theta0_X

#preliminary check:
if mag(vel) >= c:
    print("Error: Velocity cannot exceed the speed of light. Currently, it is",mag(vel)/c,"times c.")
    sys.exit("Stopped program.")


#print initial conditions:
print()
print()
print("=========================PARTICLE TRAIL INFO=========================")
print("Your Initial Parameters: \nq_b =",q_b,"A*kpc     m =",m,"kg     dt =",dt,"s     distance to track =",distance_to_track,"kpc","KE =",KE,"GeV")
print("initial position (kpc)   =",pos_0,"\ninitial velocity (kpc/s) =",vel_0)
print()



#ok, let's start tracking the monopole. Most of the calculations in the loop are to find the bfield.

while distance_tracked<distance_to_track:
    #definitions:
    r=math.sqrt(pos[0]**2+pos[1]**2)
    theta=math.atan(pos[1]/pos[0])
    gam=1/math.sqrt(1-(mag(vel)/c)**2)

    #now for bfield calculation, component by component: halo, then X-field, then disk (striated not currently used)

    #halo component:
    if pos[2]>=0:
        bfield = math.e**(-abs(pos[2])/z_0)*L(pos[2],h_disk,w_disk)* B_n*(1-L(r,r_n,w_h)) * np.array([-1*pos[1]/r, pos[0]/r,0])
    else:
        bfield = math.e**(-abs(pos[2])/z_0)*L(pos[2],h_disk,w_disk)* B_s*(1-L(r,r_s,w_h)) * np.array([-1*pos[1]/r, pos[0]/r,0])

    #X-field component:
    if abs(pos[2])>=math.tan(Theta0_X)*math.sqrt(r)-math.tan(Theta0_X)*rc_X:
        bfield += b_X(r_p(r,pos[2]))*(r_p(r,pos[2])/r)**2*np.array([math.cos(Theta_X(r,pos[2]))**2,
                                                                    math.sin(Theta_X(r,pos[2]))*math.cos(Theta_X(r,pos[2])),
                                                                    math.sin(Theta_X(r,pos[2]))])
    else:
        bfield += b_X(r_p(r,pos[2]))*(r_p(r,pos[2])/r)*np.array([math.cos(Theta_X(r,pos[2]))**2,
                                                                math.sin(Theta_X(r,pos[2]))*math.cos(Theta_X(r,pos[2])),
                                                                math.sin(Theta_X(r,pos[2]) ) ])

    #disk component:
    if r>=3.0 and r< 5.0 :
        bfield+=b_ring*(1-L(pos[2],h_disk,w_disk))
    elif r>=5.0 and r<=20.0:
        if r>=r1(theta) and r<r2(theta): #region 6
            bfield+=(b6/r)*(1-L(pos[2],h_disk,w_disk))*np.array([math.sin(math.radians(11.5))*pos[0]/r - math.cos(math.radians(11.5))*pos[1]/r,
                                                                 math.sin(math.radians(11.5))*pos[1]/r + math.cos(math.radians(11.5))*pos[0]/r,
                                                                 1])
        elif r>=r2(theta) and r<r3(theta): #region 7
            bfield+=(b7/r)*(1-L(pos[2],h_disk,w_disk))*np.array([math.sin(math.radians(11.5))*pos[0]/r - math.cos(math.radians(11.5))*pos[1]/r,
                                                                 math.sin(math.radians(11.5))*pos[1]/r + math.cos(math.radians(11.5))*pos[0]/r,
                                                                 1])
        elif r>=r3(theta) and r<r4(theta): #region 8
            bfield+=(b8/r)*(1-L(pos[2],h_disk,w_disk))*np.array([math.sin(math.radians(11.5))*pos[0]/r - math.cos(math.radians(11.5))*pos[1]/r,
                                                                 math.sin(math.radians(11.5))*pos[1]/r + math.cos(math.radians(11.5))*pos[0]/r,
                                                                 1])
        elif r>=r4(theta) and r<r5(theta): #region 1
            bfield+=(b1/r)*(1-L(pos[2],h_disk,w_disk))*np.array([math.sin(math.radians(11.5))*pos[0]/r - math.cos(math.radians(11.5))*pos[1]/r,
                                                                 math.sin(math.radians(11.5))*pos[1]/r + math.cos(math.radians(11.5))*pos[0]/r,
                                                                 1])
        elif r>=r5(theta) and r<r6(theta): #region 2
            bfield+=(b2/r)*(1-L(pos[2],h_disk,w_disk))*np.array([math.sin(math.radians(11.5))*pos[0]/r - math.cos(math.radians(11.5))*pos[1]/r,
                                                                 math.sin(math.radians(11.5))*pos[1]/r + math.cos(math.radians(11.5))*pos[0]/r,
                                                                 1])
        elif r>=r6(theta) and r<r7(theta): #region 3
            bfield+=(b3/r)*(1-L(pos[2],h_disk,w_disk))*np.array([math.sin(math.radians(11.5))*pos[0]/r - math.cos(math.radians(11.5))*pos[1]/r,
                                                                 math.sin(math.radians(11.5))*pos[1]/r + math.cos(math.radians(11.5))*pos[0]/r,
                                                                 1])
        elif r>=r7(theta) and r<r8(theta): #region 4
            bfield+=(b4/r)*(1-L(pos[2],h_disk,w_disk))*np.array([math.sin(math.radians(11.5))*pos[0]/r - math.cos(math.radians(11.5))*pos[1]/r,
                                                                 math.sin(math.radians(11.5))*pos[1]/r + math.cos(math.radians(11.5))*pos[0]/r,
                                                                 1])
        elif r>=r8(theta) and r<r1(theta): #region 5
            bfield+=(b5/r)*(1-L(pos[2],h_disk,w_disk))*np.array([math.sin(math.radians(11.5))*pos[0]/r - math.cos(math.radians(11.5))*pos[1]/r,
                                                                 math.sin(math.radians(11.5))*pos[1]/r + math.cos(math.radians(11.5))*pos[0]/r,
                                                                 1])
    #striated fields (unfinished, unused):
    #zeroorone=randrange(2)
    #if zeroorone==0:
    #    bfield-= math.sqrt(beta)*bfield
    #if zeroorone==1:
    #    bfield+= math.sqrt(beta)*bfield

    #CALCULATION OF THE CHANGE IN POSITION:
    #nonrelativistic:
    #acc=bfield*(q_b/m)

    #pos=np.array([pos[0]-vel[0]*dt-0.5*acc[0]*(dt**2),
    #              pos[1]-vel[1]*dt-0.5*acc[1]*(dt**2),
    #              pos[2]-vel[2]*dt-0.5*acc[2]*(dt**2)])
    #distance_tracked+=math.sqrt((vel[0]*dt+0.5*acc[0]*(dt**2))**2+(vel[1]*dt+0.5*acc[1]*(dt**2))**2+(vel[2]*dt+0.5*acc[2]*(dt**2))**2)

    #vel=np.array([vel[0]-acc[0]*dt,
    #              vel[1]-acc[1]*dt,
    #              vel[2]-acc[2]*dt])

    #trailx=trailx+(pos[0],)
    #traily=traily+(pos[1],)
    #trailz=trailz+(pos[2],)



    #KE=9.521406e38*6.24150934e15*gam*m*c**2 #calculate KE, and convert from kg*kpc^2/s^2 to J to MeV
    #KEhistory=KEhistory+(KE,)

    #RELATIVISTIC:

    force=q_b*bfield #In kg*kpc/s^2

    acc_prefactor=(gam*m*(c**2+gam**2*mag(vel)**2))**-1
    acc=np.array([acc_prefactor*(force[0]*(c**2+gam**2*vel[1]**2+gam**2*vel[2]**2) - gam**2*vel[0]*(force[1]*vel[1]+force[2]*vel[2])),
                  acc_prefactor*(force[1]*(c**2+gam**2*vel[0]**2+gam**2*vel[2]**2) - gam**2*vel[1]*(force[0]*vel[0]+force[2]*vel[2])),
                  acc_prefactor*(force[2]*(c**2+gam**2*vel[0]**2+gam**2*vel[1]**2) - gam**2*vel[2]*(force[0]*vel[0]+force[1]*vel[1]))])

    vel_i=vel

    vel+= -acc*dt

    pos+= -vel_i*dt-0.5*acc*dt**2

    KE=m*c**2*(gam-1)*5.942795e48 #converted to GeV from kg*kpc^2/s^2
    KEhistory=KEhistory+(KE,)

    time+=dt

    if random.randint(1,100000)==1:
        print("distance traveled:",distance_tracked)
        print("pos",pos)
        print("vel",vel)
        print("KE",KE)
        print("time",time)

    distance_tracked+=math.sqrt((vel_i[0]*dt+0.5*acc[0]*dt**2)**2+(vel_i[1]*dt+0.5*acc[1]*dt**2)**2+(vel_i[2]*dt+0.5*acc[2]*dt**2)**2)

    trailx=trailx+(pos[0],)
    traily=traily+(pos[1],)
    trailz=trailz+(pos[2],)

#print(trailx)
#print()
#print(traily)
#print()
#print(trailz)

print(pos_0)
distance_from_start=math.sqrt( (pos[0]-pos_0[0])**2 +(pos[1]-pos_0[1])**2 +(pos[2]-pos_0[2])**2)

print("The final position (kpc) is ( " + str(pos[0]) + ", " + str(pos[1]) + ", " + str(pos[2]) + ")." )
print("The final velocity (kpc/s) is ( " + str(vel[0]) + ", " + str(vel[1]) + ", " + str(vel[2]) + ")." )
print("The final Kinetic Energy (GeV) is",KE)
print()
print("Distance from initial position is", distance_from_start,"kpc")
print("The journey took",time,"seconds.")
print()
print("The galactic center is plotted as a blue dot, and the sun is plotted as a yellow dot.")
print()
print()

fig = plt.figure(figsize=(5,8))
fig.canvas.set_window_title('Monopole Trail')

ax = fig.add_subplot(211, projection='3d')
ax.plot(trailx,traily,trailz)
ax.plot([-8.33937979],[0],[0],'yo')
#ax.plot([0],[0],[0],'bo')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.suptitle("Particle Trail (click and drag to change view)",fontsize=12,position=(0.5,0.93),weight='bold')
plt.title('Particle Kinetic Energy vs Time', position=(0.5,-0.2),fontsize=12,weight='bold')

t_array=np.arange(0,dt*len(KEhistory),dt)
ax2=fig.add_subplot(212)
ax2.plot(t_array,KEhistory)
ax2.set_xlabel("Time (s)")
ax2.set_ylabel("Particle Kinetic Energy (GeV)")
plt.grid(True)

plt.show()

Answer 1

Once you do pos = pos_0, changing pos changes pos_0 because numpy arrays much like lists are mutable.

Here's a simplified example:

>>> a = [1, 2, 3]
>>> b = a
>>> b.append(4)
>>> b
[1, 2, 3, 4]
>>> a
[1, 2, 3, 4]

In python, when you create a list by saying a = [1, 2, 3], what's happening is that a list object is created and an a tag is put on it. When you assign b to a, the same object now gets a b tag on it too. So, a and b refer to the same object and changing one changes the other.

Answer 2

You are defining pos as an alias to pos_0. Later, you are making in in-place modification of pos with this line:

pos+= -vel_i*dt-0.5*acc*dt**2

Answer 3

Python uses objects so that after you say pos = pos0, any change to either name changes the object to which they both point. Thus, you are changing one when you refer to the other. I had this same problem with dictionary entries and had to use the copy method to handle it.

Answer 4

OP here.

These solutions don't actually work for numpy arrays, only for regular lists. For numpyarrays, the solution is to use numpy.copyto(destination,source).