Python. How to constrain the particles to only move on the surface of the sphere?-CodePudding

I am using vpython to simulate the movement of the particles on the sphere when repelled by other particles located on a plate below the sphere. I calculate each force using the following code.

  for i in range(len(charOS)): #COS
    for j in range(len(charges_on_plate)): #COP
      d = sqrt((charOS[i].pos.x - charges_on_plate[j].pos.x)**2   (charOS[i].pos.y - charges_on_plate[j].pos.y)**2   (charOS[i].pos.z - charges_on_plate[j].pos.z)**2)
      newFx = (con/(d**2)) * np.cos(np.arctan((charOS[i].pos.y - charges_on_plate[j].pos.y)/(charOS[i].pos.x - charges_on_plate[j].pos.x)))
      newFy = (con/(d**2)) * np.sin(np.arctan((charOS[i].pos.y - charges_on_plate[j].pos.y)/(charOS[i].pos.x - charges_on_plate[j].pos.x)))
      newFz = (con/(d**2)) * np.sin(np.arcsin((charOS[i].pos.z - charges_on_plate[j].pos.z)/(d)))
      Fx = Fx   newFx
      Fy = Fy   newFy
      Fz = Fz   newFz
    list_Fx.append(Fx)
    list_Fy.append(Fy)
    list_Fz.append(Fz)

Then I want to move each particle on the sphere based on the x, y, z forces. But they have to be necessarily constrained to the surface of the sphere. I have developed this code but it doesn't seem to work well.

My code for their constraint:

  for i in range(len(charOS)):
    if charOS[i].pos.x   radius/60 >= radius - dr*1.01 or charOS[i].pos.x   radius/60 <= -1*radius   dr*1.01:
      charOS[i].velocity = S_C_V*vector(list_Fx[i]*0, list_Fy[i], list_Fz[i])
      if charOS[i].pos.y   radius/60 >= posvec1.y - dr*1.01 or charOS[i].pos.y   radius/60 <= dr*1.01   posvec2.y:
        charOS[i].velocity = S_C_V*vector(list_Fx[i]*0, list_Fy[i]*0, list_Fz[i])
        if charOS[i].pos.z   radius/60 >= radius - dr*1.01 or charOS[i].pos.z   radius/60 <= dr*1.01   -1*radius:
          charOS[i].velocity = S_C_V*vector(list_Fx[i]*0, list_Fy[i]*0, list_Fz[i]*0)
    else:
      charOS[i].velocity = S_C_V*vector(list_Fx[i], list_Fy[i], list_Fz[i])
      if charOS[i].pos.y   radius/60 >= posvec1.y - dr*1.01 or charOS[i].pos.y   radius/60 <= dr*1.01   posvec2.y:
        charOS[i].velocity = S_C_V*vector(list_Fx[i], list_Fy[i]*0, list_Fz[i])
        if charOS[i].pos.z   radius/60 >= radius - dr*1.01 or charOS[i].pos.z   radius/60 <= dr*1.01   -1*radius:
          charOS[i].velocity = S_C_V*vector(list_Fx[i], list_Fy[i]*0, list_Fz[i]*0)
      else:
        charOS[i].velocity = S_C_V*vector(list_Fx[i], list_Fy[i], list_Fz[i])
        if charOS[i].pos.z   radius/60 >= radius - dr*1.01 or charOS[i].pos.z   radius/60 <= dr*1.01   -1*radius:
          charOS[i].velocity = S_C_V*vector(list_Fx[i], list_Fy[i], list_Fz[i]*0)
        else:
          charOS[i].velocity = S_C_V*vector(list_Fx[i], list_Fy[i], list_Fz[i])
    charOS[i].pos = charOS[i].pos   charOS[i].velocity*deltat

I know it could be fairly messy but I genuinely have no clue how to optimize this and get the desired result, thanks for the suggestions!!

CodePudding user response：

I have finally found an answer. In fact, it is fairly easy, at first I allow the new position to be calculated outside the sphere radius and then I can perform a vector normalization by dividing each vector x, y, z by the magnitude of the vector. This gives a unit vector and then finally the unit vector is multiplied by the desired sphere's radius. The code for this:

  for i in range(len(charOS)):
    charOS[i].velocity = S_C_V*vector(list_Fx[i], list_Fy[i], list_Fz[i])
    charOS[i].pos = charOS[i].pos   charOS[i].velocity*deltat
    point_vector = charOS[i].pos/mag(charOS[i].pos)
    charOS[i].pos = point_vector * radius

The calculation is significantly shorter than my previous attempt and fairly simple. :)