May seem like a silly question is there any way to create variables depending on the number of objects in a list. The context to this is I am trying to create an n body simulation. One in which the user can manually chose the planets they want in the solar system and then run the script with only those planets. So prior to them running the chosing the planets and running the script I will not know how many variables to create. This problem of how many variables to create is show by: The mass is created by a class objects: class Objects():
position = np.full((1, 3), 0)
velocity = np.full((1, 3), 0)
acceleration = np.full((1, 3), 0)
name = ""
mass = np.full((1, 1), 0)
planets_init = np.full((1, 3), 0)
def __init__(self, Name, Mass, initPosition, initVelocity, initAcceleration):
au = 149597870.700e3
v_factor = 1731460
self.name = Name
self.mass = Mass
The function is solved by using scipy.integrate.solve_ivp by:
three_body_sol = sci.integrate.solve_ivp(fun=Objects.ThreeBodyEquations,t_span=domain,y0=init_params,args=(G,planets_mass,N), max_step=max_step)
Where the function is:
def ThreeBodyEquations(t,w,G,mass,N):
# N being the number of objects in the system so in this case 2
m1, m2 = mass
#Unpack all the variables from the array "w"
r1=w[:3]
r2=w[3:6]
v1=w[6:9]
v2=w[9:12]
# Harry's attempt
G = G
planets_pos = np.vstack((r1, r2))
planets_mass = mass # np.vstack((m1, m2))
# positions r = [x,y,z] for all particles
x = planets_pos[:,0:1]
y = planets_pos[:,1:2]
z = planets_pos[:,2:3]
# matrix that stores all pairwise particle separations: r_j - r_i
dx = x.T - x
dy = y.T - y
dz = z.T - z
# matrix that stores 1/r^3 for all particle pairwise particle separations
inv_r3 = (dx**2 + dy**2 + dz**2)
inv_r3[inv_r3>0] = inv_r3[inv_r3>0]**(-1.5)
ax = G * (dx * inv_r3) @ planets_mass
ay = G * (dy * inv_r3) @ planets_mass
az = G * (dz * inv_r3) @ planets_mass
# planets_acceleration = np.sqrt(ax**2 + ay**2 + az**2)
planets_acceleration = np.vstack((ax,ay,az))
planets_acceleration = planets_acceleration.flatten()
dv1bydt=planets_acceleration[0::N]
dv2bydt=planets_acceleration[1::N]
# print(planets_acceleration)
dr1bydt=v1
dr2bydt=v2
#Package the derivatives into one final size-18 array
r12_derivs=np.concatenate((dr1bydt,dr2bydt))
r_derivs = r12derivs
v12_derivs=np.concatenate((dv1bydt,dv2bydt))
v_derivs= v12_derivs
derivs=np.concatenate((r_derivs,v_derivs))
return derivs
My main question centres around this function. When the user defines what planets they want to use I have no idea what the number of planets will be. I know the range which is might be but that’s all. Is there a feasible way to do this or is this a lost cause?
I hope this clarifys the question with some additional code. Sorry about the previous lack of code its, I didn’t realise it was valuable at first and didn’t want to burden with too much code.