LSL Script : Star Computer - 3 Body

// 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 *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); }      } }