Here is the implementation of the long range interaction.
1. We need to use the void* data that is in every cpBody structure
2. This pointer will point a new structure called Sing (from singularity). Here is an optional code you can do it as you wish.
Code: Select all
typedef struct paramActSingularity {
/** @brief The type of singularity */
char* type;
/** @brief Scalar value of singularity */
cpFloat value;
/** @brief Position in the actor in reference to the CoM (body coordinates)*/
vect position;
} acSingularity;
typedef struct ActorSingularity{
/** @brief Number of singularities*/
int Nsing;
/** @brief Array of singularities*/
acSingularity* sing;
/** @brief Global position of the singularities*/
cpVect* Gpos;
}Sing;
3. Add this code to your function run() or equivalent
Code: Select all
//Reset forces
cpSpaceEachBody(space, &eachBody, NULL);
// Long range interaction N*N/2 instead of N*N
for(int j=0; j<bodies->num; j++)
{
for(int k=j+1; k<bodies->num; k++)
{
LRangeForceApply((cpBody*)bodies->arr[j],(cpBody*)bodies->arr[k]);
}
}
cpSpaceStep(space, dt);
Code: Select all
void LRangeForceApply(cpBody *a, cpBody *b){
Sing *aux = (Sing*)a->data;
Sing *aux2 = (Sing*)b->data;
cpVect delta;
cpFloat r2=0.0;
cpVect f;
// Calculate the forces between the singularities of different bodies
for (int i=0; i<aux->Nsing; i++)
{
for (int j=0; j<aux2->Nsing; j++)
{
delta = cpvsub(aux2->Gpos[j],aux->Gpos[i]);
r2 = cpvlengthsq(delta);
// Coulomb interaction -k*q1*q2/r^2
f=cpvmult(cpvnormalize(delta),-aux->sing[i].value*aux2->sing[j].value/r2);
//Force applied to body A
cpBodyApplyForce(a,f,aux->Gpos[i]);
//Reaction on body B
cpBodyApplyForce(b,cpvneg(f),aux2->Gpos[j]);
}
}
}
4c. Observe that the singularities of the same body do not produce forces.
5. To add the Singularities to the bodies do the following (for example)
Code: Select all
Sing* aux;
aux->Nsing=1; // Here the number of singularities per body
acSingularity singularity;
singularity.type="electrical";
singularity.value=1.0e0;
singularity.position=cpv(0.0e0,1.0e0); //Slightly above the center of mass
cpArray* bodies=space->bodies;
for (int i=0; i<bodies->num; i++)
{
// Memory allocation
aux=(Sing*)malloc(sizeof(Sing));
aux->sing=(acSingularity*)malloc(aux->Nsing*sizeof(acSingularity));
aux->Gpos=(cpVect*)malloc(aux->Nsing*sizeof(cpVect));
cpBody *body = (cpBody *)bodies->arr[i];
// Loop over the singularities of each body
for(int k=0;k<aux->Nsing;k++)
{
aux->sing[k]=singularity[k];
aux->Gpos[k]=cpvadd(body->p,cpv(aux->sing[k].position.x,aux->sing[k].position.y));
}
body->data=aux;
body->position_func=SingularBodyUpdatePosition;
}
Code: Select all
void ChargedBodyUpdatePosition(cpBody *body, cpFloat dt)
{
/** @todo Increase the order of the solvers */
cpVect dp = cpvmult(cpvadd(body->v, body->v_bias), dt);
body->p = cpvadd(body->p, dp);
cpBodySetAngle(body, body->a + (body->w + body->w_bias)*dt);
// Update position of the charges
aux = (Sing*)body->data;
for (int i=0; i<aux->Nsing; i++)
aux->Gpos[i]=cpvadd(body->p,cpvrotate(cpv(aux->sing[i].position.x,
aux->sing[i].position.y), body->rot));
body->v_bias = cpvzero;
body->w_bias = 0.0f;
}
A description of the basic forces you would like to apply can be found here
http://en.wikipedia.org/wiki/Force
Enjoy!