I'm attempting to make an animation that models the 3-body problem, with not
much success. I no matter what I set the starting parameters to, nothing
ever goes into orbit, which leads me to believe that my physics are screwed
up.

The source is:
#include "colors.inc"
#include "math.inc"

#declare massA = .5;
#declare massB = .5;
#declare massC = .5;

global_settings {max_trace_level 5}

camera {location <10,-20,10> look_at <-10,20,-10> sky <0, 0, 1>}


light_source {<-140,-200,300> rgb <1.0, 1.0, 0.95>*1.5}
light_source {<140,-200,-300> rgb <0.9, 0.9, 1.00>*0.8 shadowless}

//axes
cylinder {<0, 0, 0>, <2, 0, 0>, .1 pigment {color Red}}
cylinder {<0, 0, 0>, <0, 2, 0>, .1 pigment {color Green}}
cylinder {<0, 0, 0>, <0, 0, 2>, .1 pigment {color Blue}}

#include "orbitdata.inc"

sphere {<px1, py1, pz1>, .5 pigment{color Cyan}}

sphere {<px2, py2, pz2>, .5 pigment{color Magenta}}

sphere {<px3, py3, pz3>, .5 pigment{color Yellow}}

#fopen nextpos "orbitdata.inc" write

/////////////////////////////////////////////////////////////////////ball A
#write (nextpos, "#declare px1= ", str((px1+vx1), 0, -1), ";n")
#write (nextpos, "#declare py1= ", str((py1+vy1), 0, -1), ";n")
#write (nextpos, "#declare pz1= ", str((pz1+vz1), 0, -1), ";nn")

#declare distSB = (px1-px2)*(px1-px2) + (py1-py2)*(py1-py2) +
(pz1-pz2)*(pz1-pz2);
#declare distSC = (px1-px3)*(px1-px3) + (py1-py3)*(py1-py3) +
(pz1-pz3)*(pz1-pz3);

#declare ax = -((massA*massB)/distSB)*sgn(px1-px2) -
((massA*massC)/distSC)*sgn(px1-px3);
#declare ay = -((massA*massB)/distSB)*sgn(py1-py2) -
((massA*massC)/distSC)*sgn(py1-py3);
#declare az = -((massA*massB)/distSB)*sgn(pz1-pz2) -
((massA*massC)/distSC)*sgn(pz1-pz3);

#write (nextpos, "#declare vx1= ", str((vx1+ax), 0, -1), ";n")
#write (nextpos, "#declare vy1= ", str((vy1+ay), 0, -1), ";n")
#write (nextpos, "#declare vz1= ", str((vz1+az), 0, -1), ";nn")

/////////////////////////////////////////////////////////////////////ball B
#write (nextpos, "#declare px2= ", str((px2+vx2), 0, -1), ";n")
#write (nextpos, "#declare py2= ", str((py2+vy2), 0, -1), ";n")
#write (nextpos, "#declare pz2= ", str((pz2+vz2), 0, -1), ";nn")

#declare distSA = (px2-px1)*(px2-px1) + (py2-py1)*(py2-py1) +
(pz2-pz1)*(pz2-pz1);
#declare distSC = (px2-px3)*(px2-px3) + (py2-py3)*(py2-py3) +
(pz2-pz3)*(pz2-pz3);

#declare ax = -((massB*massA)/distSA)*sgn(px2-px1) -
((massB*massC)/distSC)*sgn(px2-px3);
#declare ay = -((massB*massA)/distSA)*sgn(py2-py1) -
((massB*massC)/distSC)*sgn(py2-py3);
#declare az = -((massB*massA)/distSA)*sgn(pz2-pz1) -
((massB*massC)/distSC)*sgn(pz2-pz3);

#write (nextpos, "#declare vx2= ", str((vx2+ax), 0, -1), ";n")
#write (nextpos, "#declare vy2= ", str((vy2+ay), 0, -1), ";n")
#write (nextpos, "#declare vz2= ", str((vz2+az), 0, -1), ";nn")

/////////////////////////////////////////////////////////////////////ball C
#write (nextpos, "#declare px3= ", str((px3+vx3), 0, -1), ";n")
#write (nextpos, "#declare py3= ", str((py3+vy3), 0, -1), ";n")
#write (nextpos, "#declare pz3= ", str((pz3+vz3), 0, -1), ";nn")

#declare distSB = (px3-px2)*(px3-px2) + (py3-py2)*(py3-py2) +
(pz3-pz2)*(pz3-pz2);
#declare distSA = (px3-px1)*(px3-px1) + (py3-py1)*(py3-py1) +
(pz3-pz1)*(pz3-pz1);

#declare ax = -((massC*massB)/distSB)*sgn(px3-px2) -
((massC*massA)/distSA)*sgn(px3-px1);
#declare ay = -((massC*massB)/distSB)*sgn(py3-py2) -
((massC*massA)/distSA)*sgn(py3-py1);
#declare az = -((massC*massB)/distSB)*sgn(pz3-pz2) -
((massC*massA)/distSA)*sgn(pz3-pz1);

#write (nextpos, "#declare vx3= ", str((vx3+ax), 0, -1), ";n")
#write (nextpos, "#declare vy3= ", str((vy3+ay), 0, -1), ";n")
#write (nextpos, "#declare vz3= ", str((vz3+az), 0, -1), ";nn")

#fclose nextpos
(end of file)

distS(X), if you haven't guessed, means distance squared.

Some of the starting parameters I used were:
#declare px1= 0;
#declare py1= 0;
#declare pz1= 0;

#declare vx1= 0;
#declare vy1= 0;
#declare vz1= 0;

#declare px2= 100;
#declare py2= 100;
#declare pz2= 100;

#declare vx2= -.3;
#declare vy2= 0;
#declare vz2= 0;

#declare px3= 0;
#declare py3= 0;
#declare pz3= 8;

#declare vx3= .2;
#declare vy3= 0;
#declare vz3= 0;
(end of file)

and

#declare px1= -3;
#declare py1= 0;
#declare pz1= 3;

#declare vx1= 0;
#declare vy1= 0;
#declare vz1= -.3;

#declare px2= 0;
#declare py2= 0;
#declare pz2= 0;

#declare vx2= 0;
#declare vy2= 0;
#declare vz2= .3;

#declare px3= 3;
#declare py3= 0;
#declare pz3= 3;

#declare vx3= .3;
#declare vy3= 0;
#declare vz3= 0;
(end of file)

PS: Does anybody know how to do formatting in this message board thing?

Post a reply to this message

"Crunchy Frog" <dkl### [at] hotmailcom> wrote in message 
news:web.4431d5d2aab249115b50195f0@news.povray.org...
> I'm attempting to make an animation that models the 3-body problem, with 
> not
> much success. I no matter what I set the starting parameters to, nothing
> ever goes into orbit, which leads me to believe that my physics are 
> screwed
> up.
>
> The source is:
>... snip ...
> #declare distSB = (px1-px2)*(px1-px2) + (py1-py2)*(py1-py2) +
> (pz1-pz2)*(pz1-pz2);
> #declare distSC = (px1-px3)*(px1-px3) + (py1-py3)*(py1-py3) +
> (pz1-pz3)*(pz1-pz3);
>
> #declare ax = -((massA*massB)/distSB)*sgn(px1-px2) -
> ((massA*massC)/distSC)*sgn(px1-px3);
> #declare ay = -((massA*massB)/distSB)*sgn(py1-py2) -
> ((massA*massC)/distSC)*sgn(py1-py3);
> #declare az = -((massA*massB)/distSB)*sgn(pz1-pz2) -
> ((massA*massC)/distSC)*sgn(pz1-pz3);
>
> ... snip ...
>
> PS: Does anybody know how to do formatting in this message board thing?
>

Hi Crunchy,

I see that, for each object, you are calculating the square of the distance 
to the other two objects.
I think it's shortly after this you start to go horribly wrong :-)

I'm not sure what the intention of the next bit of your SDL is, but I think 
you may be thinking that you're calculating the delta in the path that the 
object would follow - but I don't think you are. I see that you have the 
product of two masses over the distance squared, which gives you the 
strength of an accelerative force in a line towards the second object, but 
you seem to be adding (or sometimes subtracting) the whole of that into 
'ax'. You do the same with the force for the third object. I don't think 
this is something you want to do because only a component of the force acts 
in the 'x' direction. At the moment, the instant the direction of object 2 
relative to object 1 passes across the plane x=0 you switch from adding the 
whole of the force in one direction to adding the whole of the force in the 
opposite direction.

I assume from the variable names that the variables you are storing are 
meant to be the position and velocity of each object. I also assume you're 
assuming units that enable you to have an arbitrary gravitational constant 
and that distance and time are in arbitrary unitary slices (because you're 
ignoring them in your calculations).

I think that to calculate the accelerative force acting on an object you 
could take a unit length vector in the direction of the second object and 
scale it to represent the force in that direction, then add it to a scaled 
vector in the direction of the third object, to give you a total vector in 
the direction of the total force whose length represents the force in that 
direction.

You should then be able to use that vector to make some of the unpleasant 
assumptions about gravity, time and distance units that would enable you to 
add components of the vector to the position and velocity vectors. Whether 
this gives you objects that orbit around each other will depend on how all 
those arbitrary assumptions about units pan out. I'd recommend doing the 
calculations on paper to work out positions and velocities that will give 
you viable orbits. I think the likelihood of hitting good combinations just 
by trial and error is probably quite low. These things don't just happen by 
chance you know (or did it?).

As for formatting of postings, I think it's just plain text.

Regards,
Chris B.

Post a reply to this message