开发者

translate this block of C physics code to python

I have only learnt python for for few months and totally a newb in C, I got a C code from the web, and I am dying to study it. But i only understand python language, so would someone can help to translate the following code to python would be great. Thanks in advance!

for(i=0; i<n; i++) { /* Foreach particle "i" ... */
  ax=0.0;
  ay=0.0;
  az=0.0;
  for(j=0; j<n; j++) { /* Loop over all particles "j" */
    dx=x[j]-x[i];
    dy=y[j]-y[i];
    dz=yz[j]-z[i];
    invr = 1.0/sqrt(dx*dx + dy*dy + dz*dz + eps);
    invr3 = invr*invr*invr;
    f=m[j]*invr3;
    ax += f*dx; /* accumulate the acceleration from gravitational attraction */
    ay += f*dy;
    az += f*dx;
  }
  xnew[i] = x[i] + dt*vx[i] + 0.5*dt*dt*ax; /* update position of particle "i" */
  ynew[i] = y[i] + dt*vy[i] + 0.5*dt*dt*ay;
  znew[i] = z[i] + dt*vz[i] + 0.5*dt*dt*az;
  vx[i] += dt*ax; /* update velocity of particle "i" */
开发者_运维百科  vy[i] += dt*ay;
  vz[i] += dt*az;
}

Thanks again!


Treat this as an ideal opportunity to learn some C, the syntax is not so unlike Python.

It will serve you well


Here's a literal translation

(Edit: was a literal translation but there are some anomalies in the original C code which @hughdbrown pointed out in comments -- looks like you're trying to study this subject from code so utterly broken it wouldn't even compile, which seems very unwise on your part, so, anyway, I'm taking the liberty of fixing the apparent errors and absurdities of the original C code in this Python "almost"-literal translation...):

import math

for i in range(n):
  ax = ay = az = 0.0
  for j in range(n):
    dx=x[j]-x[i]
    dy=y[j]-y[i]
    dz=z[j]-z[i]
    invr = 1.0/math.sqrt(dx*dx + dy*dy + dz*dz + eps)
    f=m[j]*invr**3
    ax += f*dx  # accumulate the acceleration from gravitational attraction
    ay += f*dy
    az += f*dz
  xnew[i] = x[i] + dt*vx[i] + 0.5*dt*dt*ax
  ynew[i] = y[i] + dt*vy[i] + 0.5*dt*dt*ay
  znew[i] = z[i] + dt*vz[i] + 0.5*dt*dt*az
  vx[i] += dt*ax  # update velocity of particle "i"
  vy[i] += dt*ay
  vz[i] += dt*az

Of course it assumes lists x, y, z, ... vx, vy, vz are all pre-existing and of the same length n, just like the C code snippet assumes for the same-named arrays.


Honestly, except for the loops this is already very close to python ..

 for(j=0; j < n; ++j) 
simply maintains the loop counter j, which is incremented after each loop iteration. The loop runs for all
0 <= j < n
. Good luck.


Incidentally this code is extremely buggy, assuming the intent is to numerically integrate the n-body problem.

  • The loop to calculate the force doesn't correctly handle the case i == j, resulting in a very large (depending on the size of eps) spurious error being added to each force.

  • Speaking of eps, for well-posed problems it serves no purpose except to add error to your force calculation.

  • The update rules for position and velocity are a strange mix of Velocity Verlet and Forward Euler, and the resulting integrator likely has all of the bad properties of the latter, at best. Velocity Verlet, being second-order, explicit, and symplectic, will give much better results for negligible cost and should be used here instead.

0

上一篇:

下一篇:

精彩评论

暂无评论...
验证码 换一张
取 消

最新问答

问答排行榜