tags:

views:

177

answers:

4

Hi all, 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!

+6  A: 

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

It will serve you well

pavium
+1, about the only thing in there that might trip him up is the for loop syntax (http://cprogramminglanguage.net/c-for-loop-statement.aspx). And they're basically the C equivalent to `for i in range(n)`. The rest is might even pass as valid python unmodified.
Adam Luchjenbroers
As Alex Martelli's translation shows...
pavium
+1  A: 

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.

Alexander Gessler
+2  A: 

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.

Alex Martelli
A literal translation -- right down to the typo in `dz=yz[j]-z[i]`.
hughdbrown
It's a typo? I thought it was referring to another list named `yz` (I even listed that among the lists that needed to be predefined in the text right after the code!). Surely the OP wouldn't be copying and pasting **broken** code that wouldn't even _compile_, would they?!
Alex Martelli
Okay, what do you think of this line of code: `az += f*dx`. The code he took it from (http://www.browndeertechnology.com/docs/BDT_OpenCL_Tutorial_NBody.html) seems dubious to me.
hughdbrown
@hughdbrown, yep, really weird (though at least that wouldn't break at compile time as the other line would if `yz` was not an existing array, so it's conceivable as a "silent bug") -- still, you've convinced me to fix the probable errors in the code, just in case (tx).
Alex Martelli
@Alex: If you're going to do that, you may as well go back to the wikipedia page (http://en.wikipedia.org/wiki/N-body_problem) and implement it properly. This code is buggy enough that I don't think it does the domain problem properly. Mind you, the OP seems to be oblivious to this.
hughdbrown
+1  A: 

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.

Re: gravity is proportional to `1/d**2` not `1/d**3`. From the wikipedia article (http://en.wikipedia.org/wiki/N-body_problem): "The forces are assumed here to be gravitational and given by Newton's law of universal gravitation; thus, they are proportional to the masses involved, and vary as the inverse square of the distance between the masses. The power of three in the denominator is correct, since it balances the vector difference in the numerator, which is necessary to specify the direction of the force." I am not a physics guy, but this explanation appears to be on point.
hughdbrown
You're absolutely right, I misread the original code.