• 1

• 11.9k
• 1
• 20
• 46

I'll second everything that @Ilmari_Karonen said. In addition, I have the following suggestions.

#Avoid Parallel Arrays Any time you find yourself using parallel arrays (as you have with x, y, theta, e_x, and e_y), that's a sign that you really need to have a data type. OOP is all about creating a data type and methods to act on that data type to separate it out from all the other code you write to make it easier to maintain.

So you should have a Point class (or use an existing one from a library of some sort). I recommend something like this:

class Point {
public:
 Point(double newX, double newY) : x(newX), y(newY) {}
 ~Point() {}
 double getX() const { return x; }
 double getY() const { return y; }
 void setX(double newX) { x = newX; }
 void setY(double newY) { y = newY; }
private:
 double x;
 double y;
};

You might then want the class to contain such useful methods as operator+(), operator-(), operator=(), operator==(), dot(const Point& x).

You might further create class like for a particle. Something like this:

class Particle {
 public:
 Particle(const Point& location, const Point& velocity);
 Particle(const Point& location, double speed, double angle);
 ~Particle();
 private:
 Point loc;
 Point vel;
 double theta;
};

Then you could create an array (or std::vector) of Particle objects. This would allow you to further simplify your code by using range-based loops to iterate over the collection.

Doing all of the above will make reading your code enormously easier.

#Performance As always, the best way to determine why your code is slow is to profile it. Anything else is just a guess. Having said that, I will take an educated guess at what might be the problem, as I've written similar code in the past and hit this problem.

You're using the pow() function to compute squares and cubes. The pow() function is very slow. You should pull out the terms that you want to square and cube and manually do it. Something like this:

double r_ij_squared = r_ij * r_ij;
double r_ij_cubed = r_ij_squared * r_ij;
double er_min = e_x[i] * r_y_ij - e_y[i] * r_x_ij;
double er_min_squared = er_min * er_min;
double er_plus = e_x[j] * r_x_ij + e_y[j] * r_y_ij;
double er_plus_squared = er_plus * er_plus;

You also have some common terms that you are continuously recalculating. It looks like you reuse r_x_ij * e_x[j] + r_y_ij * e_y[j] in the addition to S_theta, then again in the next line when setting S.

These lines can be simplified:

 while(theta[i]>pi)
 theta[i]-=2*pi;
 while(theta[i]<-pi)
 theta[i]+=2*pi;

This is basically a modulus operator. You can do something like:

theta[i] = fmod(theta[i], 2 * pi);
if (theta[i] < 0.0)
{
 theta[i] += 2 * pi;
}

I believe that accomplishes the same thing without a loop. If theta is 27000π, it will execute in the same amount of time as if theta was between 0 and 2π.