From Micasim
// Star Computer - 3-Body v1.2 by Rob Knop (Prospero Frobozz)
// Copyright (C) 2009 Prospero Frobozz
//
// This script is free software, available under the GNU GPL version 2.
// You are free to use it, modify it, or distribute it. If you distribute
// this script or any modified versions, they must also be under the GNU GPL.
// This means that the script itself must be mod/copy/transfer. Objects that
// use the script, or other components of such objects, do not themselves need
// to be mod/copy/transfer. However, if you distribute those objects, the
// script, or any other scripts derived from it, must, and all derived scripts
// must be made available under the GNU GPL.
// 1m in SL vectors = 1AU in simulation
// time measured in years
// Velocities in AU/year
// Mass in Solar Masses
// G = 39.4776 AU^3 Msun^-1 year^-2
float G=39.4776;
integer nstars = 3; // Set this in state_entry() for state default
// dispscale converts AU to actual SL m displayed : position = ctrpos + displscale * pos
float dispscale = 0.2;
// The simulation runs within the link_message call back.
// Each link_message, it will do subiterations time steps,
// and then send out updates to the star visualizer.
// This is based on the idea that it *should*, ideally, be able
// to do a bunch of iterations between the script delays that will
// be associated with updating the star positions.... one can dream,
// anyway! LSL is very slow, and really the wrong language in which
// to write this sort of thing, but, hey, I'm doing it anyway.
// This version is HARD-CODED for 3 bodies. Reason : LSL has no arrays,
// and I wanted to avoid the overheads of list processing.
integer subiterations = 75; // This many iterations between updates of the star positions
integer energyreport = 20; // This many link messages between energy report
integer doreportenergy = 1; // Flag : report energy out loud or not?
float defaulttimestep = 0.001;
float eps = 0.05; // Potential softening length
// ctrpos is the center of the place where stars should be positioned,
// and stars will never go outside ctrpos +- <1., 1., 1.>*posrange*dispscale
// NOTE -- ctrpos will be set whenever the computer object is rezzed or moved.
// It will be centered on the center position of the computer object, offset by
// tankoffset (to get past the corner of the computer object), offset again
// by <posrange, posrange, posrange>*dispscale.
vector ctrpos = <24., 21., 132.>;
float posrange = 15.; // half-size of display area / posrange
vector tankoffset = <0.75, 0.75, -0.25>;
// The following are the channels on which these listen
// Stars 1-3 listen to "Star Computer", "Star Updater", and the owner
// on channels 1-3.
integer compchan = 16384; // Commands to the star computer here
integer updates = 4096; // Updates about initial conditions shouted here
// State variables and some globals for returning derivative values
// (sad)
vector pos1;
vector pos2;
vector pos3;
vector vel1;
vector vel2;
vector vel3;
vector dpos1dt;
vector dpos2dt;
vector dpos3dt;
vector dvel1dt;
vector dvel2dt;
vector dvel3dt;
vector initpos1;
vector initpos2;
vector initpos3;
vector initvel1;
vector initvel2;
vector initvel3;
float mass1;
float mass2;
float mass3;
float initmass1;
float initmass2;
float initmass3;
float minmass = 0.1;
float maxmass = 10.0;
float initdt;
float dt;
float t;
integer running;
integer energycounter;
// ****************************************
// Reset the state of the simulation
initialize()
{
llSay(0, "Initalizing...");
integer i;
ctrpos = llGetPos() + tankoffset + <1., 1., 1.>*posrange*dispscale;
mass1 = initmass1;
mass2 = initmass2;
mass3 = initmass3;
pos1 = initpos1;
vel1 = initvel1;
pos2 = initpos2;
vel2 = initvel2;
pos3 = initpos3;
vel3 = initvel3;
dt = initdt;
t = 0.;
energycounter = 0;
// Tell stars about their masses, and where
// the vizualization takes place.
llRegionSay(1, "reset");
llRegionSay(1, "CTRPOS " + (string)ctrpos);
llRegionSay(1, "POSRANGE " + (string)posrange);
llRegionSay(1, "DISPSCALE " + (string)dispscale);
llRegionSay(1, "MASS " + (string)mass1);
llRegionSay(2, "reset");
llRegionSay(2, "CTRPOS " + (string)ctrpos);
llRegionSay(2, "POSRANGE " + (string)posrange);
llRegionSay(2, "DISPSCALE " + (string)dispscale);
llRegionSay(2, "MASS " + (string)mass2);
llRegionSay(3, "reset");
llRegionSay(3, "CTRPOS " + (string)ctrpos);
llRegionSay(3, "POSRANGE " + (string)posrange);
llRegionSay(3, "DISPSCALE " + (string)dispscale);
llRegionSay(3, "MASS " + (string)mass3);
// call stopsim to make sure we aren't running,
// and also (incidentally) to position the stars
stopsim();
}
// ****************************************
// Stop the simulation running
stopsim()
{
running = 0;
updatestarpos();
updatefloatingtext();
}
// ****************************************
// Start the simulation running.
startsim()
{
running = 1;
updatestarpos();
updatefloatingtext();
llMessageLinked(LINK_THIS, 1, "", "");
if (energyreport > 0) {
llSay(0, "Initial Energy : " + (string)total_energy());
}
}
// ****************************************
updatefloatingtext()
{
string text = "";
if (running) {
text = "Running...\n";
}
llMessageLinked(LINK_THIS, 10, text + "t = " + (string)t, NULL_KEY);
}
// ****************************************
// Update the stars with where they are
updatestarpos()
{
integer i;
vector pos;
vector vel;
string text = "";
llRegionSay(1 , "POS " + (string)pos1);
llRegionSay(1 , "VEL " + (string)vel1);
llRegionSay(2 , "POS " + (string)pos2);
llRegionSay(2 , "VEL " + (string)vel2);
llRegionSay(3 , "POS " + (string)pos3);
llRegionSay(3 , "VEL " + (string)vel3);
if (energyreport > 0)
{
if (++energycounter >= energyreport)
{
energycounter = 0;
if (doreportenergy)
{
llSay(0, "At t=" + (string)t + ", " +
"Total Energy: " + (string)total_energy());
}
updatefloatingtext();
}
}
}
// ****************************************
// Return the time derivative of each param
// in globals dposNdt and dvelNdt
derivs(vector p1, vector v1, vector p2, vector v2,
vector p3, vector v3)
{
vector fperm;
dpos1dt = v1;
dpos2dt = v2;
dpos3dt = v3;
dvel1dt = <0., 0., 0.>;
dvel2dt = <0., 0., 0.>;
dvel3dt = <0., 0., 0.>;
float d12 = llVecDist(p1, p2);
float d12soft2 = d12*d12 + eps*eps;
float d13 = llVecDist(p1, p3);
float d13soft2 = d13*d13 + eps*eps;
float d23 = llVecDist(p2, p3);
float d23soft2 = d23*d23 + eps*eps;
fperm = G / llPow(d12soft2, 1.5) * (p1-p2);
dvel1dt -= mass2 * fperm;
dvel2dt += mass1 * fperm;
fperm = G / llPow(d13soft2, 1.5) * (p1-p3);
dvel1dt -= mass3 * fperm;
dvel3dt += mass1 * fperm;
fperm = G / llPow(d23soft2, 1.5) * (p2-p3);
dvel2dt -= mass3 * fperm;
dvel3dt += mass2 * fperm;
}
// ****************************************
// Runge-Kutta method
//
// I wish LSL had real arrays
//
// Note : Using globals here since it seems
// that in LSL, things pass by value
advance()
{
vector p1 = pos1;
vector p2 = pos2;
vector p3 = pos3;
vector v1 = vel1;
vector v2 = vel2;
vector v3 = vel3;
vector k1p1;
vector k1p2;
vector k1p3;
vector k2p1;
vector k2p2;
vector k2p3;
vector k3p1;
vector k3p2;
vector k3p3;
vector k4p1;
vector k4p2;
vector k4p3;
vector k1v1;
vector k1v2;
vector k1v3;
vector k2v1;
vector k2v2;
vector k2v3;
vector k3v1;
vector k3v2;
vector k3v3;
vector k4v1;
vector k4v2;
vector k4v3;
derivs(p1, v1, p2, v2, p3, v3);
k1p1 = dpos1dt;
k1p2 = dpos2dt;
k1p3 = dpos3dt;
k1v1 = dvel1dt;
k1v2 = dvel2dt;
k1v3 = dvel3dt;
p1 += 0.5*k1p1*dt;
p2 += 0.5*k1p2*dt;
p3 += 0.5*k1p3*dt;
v1 += 0.5*k1v1*dt;
v2 += 0.5*k1v2*dt;
v3 += 0.5*k1v3*dt;
derivs(p1, v1, p2, v2, p3, v3);
k2p1 = dpos1dt;
k2p2 = dpos2dt;
k2p3 = dpos3dt;
k2v1 = dvel1dt;
k2v2 = dvel2dt;
k2v3 = dvel3dt;
p1 = pos1 + 0.5*k2p1*dt;
p2 = pos2 + 0.5*k2p2*dt;
p3 = pos3 + 0.5*k2p3*dt;
v1 = vel1 + 0.5*k2v1*dt;
v2 = vel2 + 0.5*k2v2*dt;
v3 = vel3 + 0.5*k2v3*dt;
derivs(p1, v1, p2, v2, p3, v3);
k3p1 = dpos1dt;
k3p2 = dpos2dt;
k3p3 = dpos3dt;
k3v1 = dvel1dt;
k3v2 = dvel2dt;
k3v3 = dvel3dt;
p1 = pos1 + k3p1*dt;
p2 = pos2 + k3p2*dt;
p3 = pos3 + k3p3*dt;
v1 = vel1 + k3v1*dt;
v2 = vel2 + k3v2*dt;
v3 = vel3 + k3v3*dt;
derivs(p1, v1, p2, v2, p3, v3);
k4p1 = dpos1dt;
k4p2 = dpos2dt;
k4p3 = dpos3dt;
k4v1 = dvel1dt;
k4v2 = dvel2dt;
k4v3 = dvel3dt;
pos1 += (k1p1 + 2.*k2p1 + 2.*k3p1 + k4p1)/6. * dt;
pos2 += (k1p2 + 2.*k2p2 + 2.*k3p2 + k4p2)/6. * dt;
pos3 += (k1p3 + 2.*k2p3 + 2.*k3p3 + k4p3)/6. * dt;
vel1 += (k1v1 + 2.*k2v1 + 2.*k3v1 + k4v1)/6. * dt;
vel2 += (k1v2 + 2.*k2v2 + 2.*k3v2 + k4v2)/6. * dt;
vel3 += (k1v3 + 2.*k2v3 + 2.*k3v3 + k4v3)/6. * dt;
t += dt;
}
// *****************************************
// Calculate the total energy in the system
float total_energy()
{
float e = 0.;
float mag;
mag = llVecMag(vel1);
e += 0.5 * mass1 * mag * mag;
mag = llVecMag(vel2);
e += 0.5 * mass2 * mag * mag;
mag = llVecMag(vel3);
e += 0.5 * mass3 * mag * mag;
mag = llVecMag(pos1 - pos2);
e -= G * mass1 * mass2 / llSqrt(mag*mag + eps*eps);
mag = llVecMag(pos1 - pos3);
e -= G * mass1 * mass3 / llSqrt(mag*mag + eps*eps);
mag = llVecMag(pos2 - pos3);
e -= G * mass2 * mass3 / llSqrt(mag*mag + eps*eps);
return e;
}
// **********************************************************************
centermomentum()
{
vector cm = ( mass1*initvel1 + mass2*initvel2 + mass3*initvel3 ) / (mass1 + mass2 + mass3);
stopsim();
initvel1 -= cm;
initvel2 -= cm;
initvel3 -= cm;
initialize();
}
centerpositions()
{
vector offset = (mass1*initpos1 + mass2*initpos2 + mass3*initpos3) /
(mass1 + mass2 + mass3);
stopsim();
initpos1 -= offset;
initpos2 -= offset;
initpos3 -= offset;
initialize();
}
shareinits()
{
llRegionSay(updates, "initpos 1 " + (string)initpos1);
llRegionSay(updates, "initpos 2 " + (string)initpos2);
llRegionSay(updates, "initpos 3 " + (string)initpos3);
llRegionSay(updates, "initvel 1 " + (string)initvel1);
llRegionSay(updates, "initvel 2 " + (string)initvel2);
llRegionSay(updates, "initvel 3 " + (string)initvel3);
llRegionSay(updates, "initmass 1 " + (string)initmass1);
llRegionSay(updates, "initmass 2 " + (string)initmass2);
llRegionSay(updates, "initmass 3 " + (string)initmass3);
}
// **********************************************************************
default
{
state_entry()
{
running = 0;
// Gotta have some default: tight binary plus a more distant
// companion. Hope I did this right.
initmass1 = 1.;
initmass2 = 1.;
initmass3 = 0.6;
initpos1 = < 1. , 0. , 0. >;
initvel1 = < 0 , 3.14156 , 0.>;
initpos2 = < -1. , 0. , 0. >;
initvel2 = < 0., -3.14156, 0.>;
initpos3 = < -14. , 0. , 0. >;
initvel3 = < 0. , 0.25 , -0.3>;
initdt = defaulttimestep;
initialize();
llListen(16384, "", "", "");
}
// moving_end() also gets called on a rez (unless it's rezzed in the same
// position, in which case we don't need to call it)
moving_end()
{
initialize();
}
// ****************************************
// Advance on receiving a link message with int 1. Do this
// instead of just calling advance() repeatedly within a loop
// so that we will head out to the main loop every step
// to check for other events that might stop us, or that
// might call the timer for updating the star display
// positions.
link_message(integer sender, integer num, string str, key id)
{
if (num != 1) return;
integer j;
if (running && num == 1)
{
for (j=0 ; j<subiterations ; ++j) {
advance();
}
llMessageLinked(LINK_THIS, 1, "", "");
updatestarpos();
}
}
// ****************************************
listen(integer chan, string name, key id, string message)
{
list params;
string com;
integer num;
integer i;
float val;
vector vec;
params = llParseString2List(message, [" "], []);
num = llGetListLength(params);
if (num == 0)
{
return;
}
com = llList2String(params, 0);
if (com == "stop")
{
stopsim();
}
else if (com == "start")
{
startsim();
}
else if (com == "initialize" || com == "reset")
{
initialize();
}
else if (com == "dt")
{
if (num < 2)
{
llSay(0, "Current dt = "+(string)dt);
if (initdt != dt)
{
llSay(0, "dt to be used after reinitialize = "+
(string)initdt);
}
}
else
{
initdt = llList2Float(params, 1);
llSay(0, "dt to be used after reinitialize = "+
(string)initdt);
}
}
else if (com == "starmass")
{
i = llList2Integer(params, 1);
val = llList2Float(params, 2);
if (val < minmass) val = minmass;
if (val > maxmass) val = maxmass;
if (i == 1) initmass1 = val;
else if (i == 2) initmass2 = val;
else if (i == 3) initmass3 = val;
}
else if (com == "starpos")
{
i = llList2Integer(params, 1);
vec = (vector)llDumpList2String(llList2List(params, 2, -1), " ");
if ( vec.x < -posrange ||
vec.x > posrange ||
vec.y < -posrange ||
vec.y > posrange ||
vec.z < -posrange ||
vec.z > posrange )
{
llSay(0, "Error, position " + (string)vec +
" is outside the simulation. Not setting.");
}
else {
if (i == 1) initpos1 = vec;
else if (i == 2) initpos2 = vec;
else if (i == 3) initpos3 = vec;
}
}
else if (com == "starvel")
{
i = llList2Integer(params, 1);
vec = (vector)llDumpList2String(llList2List(params, 2, -1), " ");
if (i == 1) initvel1 = vec;
else if (i == 2) initvel2 = vec;
else if (i == 3) initvel3 = vec;
}
else if (com == "reportenergy")
{
doreportenergy = 1;
}
else if (com == "noreportenergy")
{
doreportenergy = 0;
}
else if (com == "centermomentum")
{
centermomentum();
shareinits();
}
else if (com == "centerpositions")
{
centerpositions();
shareinits();
}
else if (com == "pushinit")
{
shareinits();
}
else if (com == "dump")
{
llSay(0, "CtrPos : " + (string)ctrpos);
llSay(0, "dispscale : " + (string)dispscale);
llSay(0, "Star masses : ");
llSay(0, " Star 1 : " + (string)mass1);
llSay(0, " Star 2 : " + (string)mass2);
llSay(0, " Star 3 : " + (string)mass3);
llSay(0, "Star positions and velocities : ");
llSay(0, " Star 1 : " + (string)pos1 + " , vel = " + (string)vel1);
llSay(0, " Star 2 : " + (string)pos2 + " , vel = " + (string)vel2);
llSay(0, " Star 3 : " + (string)pos3 + " , vel = " + (string)vel3);
llSay(0, "Time Step = " + (string)dt);
llSay(0, "Time = " + (string)t);
llSay(0, "Energy = " + (string)total_energy());
}
}
}
