Planet simulation with pass by reference - force vector won't update!

[programmingphysicist]

Hi there,

I'm fairly new to this, but I'm using C++ to simulate a basic planet orbiting around a sun. I'm using a 'leapfrog simulation' to update the position and velocity of the planet as time goes on, the details of which aren't super important. Point is, the code doesn't work, with the problem being that the force vector won't update, even as I call the function which updates it using pass by reference in the for loop in main(). Can anybody spot an issue?

#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

#define D 2 // number of dimensions 

struct planet {
	double x[D]; // position
	double v[D]; // velocity
	double F[D]; // force 
	double A; // GMm for the system
	double m; // mass of planet 
	double K; // kinetic energy of planet
	double V; //potential energy of planet
};

double distance (planet a);
double VelocityMag (planet a);
void Force (planet &a);
void energy (planet &a);

int main()
{
	planet a;
	int N=1000; // number of increments in the simulation 
	double t = 0.0; // set initial time t
	double dt = 0.1; // set time increment dt
	
	a.A = -8e44; // set GMm
	a.m = 6e24; // set m
	
	a.x[0] = 0.0; // initial conditions
	a.x[1] = 100000000000;
	a.v[0] = 30000;
	a.v[1] = 0.0; 
	
	ofstream fout;
	fout.open("output.dat");
	
	for (int l=0; l<N; l++)
	{
		
		for (int n=0; n<D; n++)
		{
			a.x[n] = a.x[n] + (a.v[n]*0.5*dt);
		}
		
		t = t + (dt/2);
		
		Force(a);
		
		for (int p=0; p<D; p++)
		{
			a.v[p] = a.v[p] + ((dt/a.m)*a.F[p]);
		}
		
		for (int q=0; q<D; q++)
		{
			a.x[q] = a.x[q] + (a.v[q]*0.5*dt);
		}
		
		t = t + (dt/2);
		
		for (int r=0; r<D; r++)
		{
			fout << a.x[r] << "\t" ;
			fout << a.F[r] << "\t" ;
		}
		
		fout << endl;
		
	}
			
	
	return 0;
}


double distance (planet a)
{
	double dsquare=0;
	for (int i=0; i<D; i++)
	{
		dsquare = dsquare + ((a.x[i])*(a.x[i]));
	}
	double d = sqrt(dsquare);
	return d;
}

double VelocityMag (planet a)
{
	double VMagSquare=0;
	for (int j=0; j<D; j++)
	{
		VMagSquare = VMagSquare + (a.v[j]*a.v[j]);
	}
	double VMag = sqrt(VMagSquare);
	return VMag;
}

void Force (planet &a)
{
	for (int k=0; k<D; k++)
	{
		a.F[k] = (a.A * a.x[k])/(pow(3,distance(a)));
	}
}

void energy (planet &a)
{
	a.K = (1/2*a.m*VelocityMag(a)*VelocityMag(a));
	a.V = a.A / distance(a);
}

[cire]

pow(3, distance(a)) results in a value of inf, which doesn't play well with the rest of your calculation.

[lastchance]

cire's point about the arguments being the wrong way round in the call to pow() is probably what is stopping your code from running.

You don't use it at the moment, but I would fix your expression for kinetic energy:
a.K = (1/2*a.m*VelocityMag(a)*VelocityMag(a));
The integer division implicit in the 1/2 bit will mean that the KE is always zero: not very acceptable if you are a physicist. (The expression is also very inefficient, since you square-rooted to find VelocityMag and now you are squaring back again.)

That also looks like an incredibly small timestep dt, given the distances and velocities involved.