RK4 Solar System

[Fizyk]

Remember my 2 contraptions - holographic Solar System (which I posted in ahref's thread) and RK4 gravity simulator? Well, I connected the two and hereby present the RK4-simulated holographic Solar System:

@name Runge-Kutta Solar System
@inputs 
@outputs R:array
@persist T Tick
#constants
@persist G N Skip SLOW
@trigger 

if(first()|duped()) {
    runOnTick(1)
    SLOW = 4
    Skip = 3
    Tick = 1
	
    G = 33.167
    R:clear()
    N=5 #5 entries per body
    
    #"Sun"
    R[1,number]=1.989e+6
    R[2,vector]=vec(70,70,70) #size
    R[3,vector]=vec(255,255,0) #color
    R[4,vector]=vec(0,0,0) #initial position
    R[5,vector]=vec(0,-1.394e-3,0) #initial velocity
    
    #"Mercury"
    R[6,number]=0.33
    R[7,vector]=vec(0.24,0.24,0.24)
    R[8,vector]=vec(70,60,0)
    R[9,vector]=vec(787.28,110.952,0) #V,p
    R[10,vector]=vec(7.005,48.331,29.124) #I,Omega,omega
    
    #"Venus"
    R[11,number]=4.869
    R[12,vector]=vec(0.61,0.61,0.61)
    R[13,vector]=vec(160,160,0)
    R[14,vector]=vec(552.1,216.4,0)
    R[15,vector]=vec(3.39,76.67,54.852)
    
    #"Earth"
    R[16,number]=5.975
    R[17,vector]=vec(0.64,0.64,0.64)
    R[18,vector]=vec(0,255,255)
    R[19,vector]=vec(469.63,299.2,0)
    R[20,vector]=vec(0,348.74,114.2)
    
    #"Mars"
    R[21,number]=0.64
    R[22,vector]=vec(0.37,0.37,0.37)
    R[23,vector]=vec(255,70,70)
    R[24,vector]=vec(383.78,451.77,0)
    R[25,vector]=vec(1.85,49.562,286.537)
    
    #"Jupiter"
    R[26,number]=1.9e3
    R[27,vector]=vec(7,7,7)
    R[28,vector]=vec(100,100,0)
    R[29,vector]=vec(206.3,1553.73,0)
    R[30,vector]=vec(1.3,100.492,275.066)
    
    T = curtime()/SLOW
    
    for(I=1,R:count()/N) {
        holoCreate(I)
        holoModel(I,"hqicosphere2")
        holoScaleUnits(I,R[I*N-3,vector])
        holoColor(I,R[I*N-2,vector])
        
        #calculate position and velocity
        if(I>1) {
            V = R[I*N-1,vector]:x()
            P = R[I*N-1,vector]:y()
            I2 = R[I*N,vector]:x()
            O = R[I*N,vector]:y()
            O2 = R[I*N,vector]:z()
        
            R[I*N-1,vector] = P*vec(cos(O2)*cos(O)-sin(O2)*sin(O)*cos(I2),sin(O)*cos(O2)+cos(O)*sin(O2)*cos(I2),sin(O2)*sin(I2))
            R[I*N,vector] = V*vec(-cos(O)*sin(O2)-sin(O)*cos(O2)*cos(I2),-sin(O)*sin(O2)+cos(O)*cos(O2)*cos(I2),cos(O2)*sin(I2))
        }
        
        holoPos(I,entity():pos()+vec(0,0,300)+R[I*N-1,vector])
        #ifdef E:setTrails(N,N,N,S,V,N)
            holoEntity(I):setTrails(30,30,R[I*N-1,vector]:length()/20,"trails/laser",R[I*N-2,vector],500)
        #endif
    }
}
elseif(Tick == 0) {
    Tick++
    if(Tick>=Skip) {Tick = 0}
	
#Simulation
    T = curtime()/SLOW
    P1 = array()
    for(I=1,R:count()/N) {      
        P1:pushVector($T*R[I*N,vector])
        P1V = vec()
        for(J=1,R:count()/N)
        {
            if(J!=I) {
                DX1 = R[J*N-1,vector]-R[I*N-1,vector]
                P1V += $T*G*R[J*N-4,number]/(DX1:length())^3*DX1
            }
        }
        P1:pushVector(P1V)
    }
    
    P2 = array()
    for(I=1,R:count()/N) { 
        P2:pushVector($T*(R[I*N,vector]+P1[I*2,vector]/2))
        P2V = vec()
        for(J=1,R:count()/N)
        {
            if(J!=I) {
                DX2 = (R[J*N-1,vector]+P1[J*2-1,vector]/2)-(R[I*N-1,vector]+P1[I*2-1,vector]/2)
                P2V += $T*G*R[J*N-4,number]/(DX2:length())^3*DX2
            }
        }
        P2:pushVector(P2V)
    }
    
    P3 = array()
    for(I=1,R:count()/N) { 
        P3:pushVector($T*(R[I*N,vector]+P2[I*2,vector]/2))
        P3V = vec()
        for(J=1,R:count()/N)
        {
            if(J!=I) {
                DX3 = (R[J*N-1,vector]+P2[J*2-1,vector]/2)-(R[I*N-1,vector]+P2[I*2-1,vector]/2)
                P3V += $T*G*R[J*N-4,number]/(DX3:length())^3*DX3
            }
        }
	       P3:pushVector(P3V)
    }
    
    P4 = array()
    for(I=1,R:count()/N) { 
        P4:pushVector($T*(R[I*N,vector]+P3[I*2,vector]))
        P4V = vec()
        for(J=1,R:count()/N)
        {
            if(J!=I) {
                DX4 = (R[J*N-1,vector]+P3[J*2-1,vector])-(R[I*N-1,vector]+P3[I*2-1,vector])
                P4V += $T*G*R[J*N-4,number]/(DX4:length())^3*DX4
            }
        }
        P4:pushVector(P4V)
    }
    
    for(I=1,R:count()/N) {    
        PX = (P1[I*2-1,vector]+P2[I*2-1,vector]*2+P3[I*2-1,vector]*2+P4[I*2-1,vector])/6
        PV = (P1[I*2,vector]+P2[I*2,vector]*2+P3[I*2,vector]*2+P4[I*2,vector])/6
        
        R[I*N-1,vector] = R[I*N-1,vector]+PX
        R[I*N,vector] = R[I*N,vector]+PV
        
        holoPos(I,entity():pos()+vec(0,0,300)+R[I*N-1,vector])
    }
}
else {
    Tick++
    if(Tick>=Skip) { Tick = 0 }
}

Screenshots in the attachment. The last two pictures show inner planets, from Mercury to Mars (or rather, their trails). The first screenshot is taken from a bit further away to include Jupiter.

I entered data for only 5 planets + Sun, because 6 bodies are near the limit of E2 performance. You can change the speed of the simulation by adjusting the SLOW constant - SLOW = 1 means one year passes in 4 seconds, I set SLOW to 4 here, which means that one year passes in 16 seconds.

You can also add your own bodies in a way similar to the one I described in the RK4 thread, except the vectors for initial position and velocity mean something different now:
1. Initial "position" = (initial velocity, distance from Sun at perihelion, 0)
2. Initial "velocity" = (inclination, longitude of ascending node, argument of perihelion)
(This makes it easier for me to enter the real parameters in the E2.)

What is the difference between this E2 and the previous one? The previous one was just moving the planets along predefined orbits. This one puts planets initially at some positions with some velocities, and then simulates Newtonian gravity - each planet (or Sun) affects every other body! Difference is probably impossible to notice, but it shows that physics works

You can mess with masses of the planets, their positions and velocities and see how all this affects the structure of the Solar System

Have fun!

[Manslaughter]

Very nice. I can finaly simulate how to destroy the world! MUWAHAHAHAAAA!!!
3 's for you.

Newton and Galileo would be prowd.

[Paper Clip]

Neat idea. I like how you simulated the physics instead of just pre-determined orbits. Good job.

[Fizyk]

Thanks

[SystemsLock]

You stole the idea from my thread! Not that I care or anything.

More interesting than preconfigured systems would be to see what happens what you start all the objects in random places with random velocities. It would be interesting to watch them form into systems and observe celestial phenomenon.

[Paper Clip]

You stole the idea from my thread! Not that I care or anything.

Quite a few solar system contraptions have been posted in this section, I wouldn't be surprised if yours wasn't the first of the bunch.

You should make the owner of the chip the sun, then you could have fun messing around with the orbits by moving around

[sebi99]

You could do that chatcmd pops up derma control wchich allows you to spawn new bodies.

[Fizyk]

You stole the idea from my thread! Not that I care or anything.

Unless you mean your old astronomical/gravitational simulations thread, which gave me the idea to make the RK4 simulation at all, I swear I got this idea by myself

You should make the owner of the chip the sun, then you could have fun messing around with the orbits by moving around

I did that once, but that messes with orbits a bit too much

You could do that chatcmd pops up derma control wchich allows you to spawn new bodies.

Nice idea. But as E2 derma is unofficial, I would rather make a chatcmd directly for adding new bodies. I might do that soon

EDIT:
Done.

@name Chat-controlled RK4 Gravitational Simulation
@inputs
@outputs R:array
@persist T Tick
#constants
@persist G N Skip SLOW Simulate
@persist Players:array Center:vector
@persist Highlighted
@trigger 

if(first())
{
    G = 33.167
    N = 5
    Tick = 1
    T = curtime()
    Skip = 1
    SLOW = 1
    
    runOnTick(1)
    runOnChat(1)
    R:clear()
    
    Players:clear()
    Players:pushString(owner():name())
    Highlighte = -1
}

#---------------------- Commands -----------------------------

if(chatClk())
{
    Priv = 0
    for(I=1,Players:count())
    {
        if(lastSpoke():name()==Players[I,string]) { Priv = 1 }
    }
    if(Priv)
    {
        S = lastSpoke():lastSaid()
        if(S[1]!="!") { exit() }
        hideChat(1)
        
        Com = S:explode(" ")
        Bodies = floor(R:count()/N)
        
        #Add privileges for a player
        if(Com[1,string]=="!addplayer")
        {
            Players:pushString(Com[2,string])
            print(Com[2,string]+" can now control the RK4 simulation.")
        }
        
        #Start/Stop simulation
        if(Com[1,string]=="!simulate")
        {
            Simulate = !Simulate
            
            for(B=1,Bodies)
            {
                if(Simulate)
                {
                    holoScaleUnits(B*2+28,vec())
                    holoScaleUnits(B*2+29,vec())
                }
                else
                {
                    holoPos(B*2+28,Center+R[B*N-1,vector]+R[B*N,vector]*5/12)
                    holoPos(B*2+29,Center+R[B*N-1,vector]+R[B*N,vector]*5/6)
                    A = (quat(R[B*N,vector]:toAngle())*qRotation(vec(0,90,0))):toAngle()
                    holoAng(B*2+28,A)
                    holoAng(B*2+29,A)
                    holoScaleUnits(B*2+28,vec(10,10,R[B*N,vector]:length()*5/6))
                    holoScaleUnits(B*2+29,vec(15,15,R[B*N,vector]:length()/6))
                }
            }
        }
        
        #Add a body
        #No parameters
        if(Com[1,string]=="!addbody")
        {
            if(Simulate) { exit() }
            holoCreate(Bodies)
            holoModel(Bodies,"hqicosphere2")
            holoScaleUnits(Bodies,vec(20,20,20))
            holoColor(Bodies,vec(255,255,255))
            holoPos(Bodies,Center)
            
            #velocity representation
            holoCreate(Bodies*2+30)
            holoCreate(Bodies*2+31)
            holoModel(Bodies*2+30,"hqcylinder2")
            holoModel(Bodies*2+31,"hqcone")
            holoScaleUnits(Bodies*2+30,vec())
            holoScaleUnits(Bodies*2+31,vec())
            holoColor(Bodies*2+30,vec(255,255,255))
            holoColor(Bodies*2+31,vec(255,255,255))
            
            R:pushNumber(1) #mass
            R:pushVector(vec(20,20,20)) #size
            R:pushVector(vec(255,255,255)) #color
            R:pushVector(vec(0,0,0)) #position
            R:pushVector(vec(0,0,0)) #velocity
        }
        
        #Highlight a body
        #Parameters: #body
        if(Com[1,string]=="!highlight")
        {
            B = Com[2,string]:toNumber()
            if(B>Bodies) { exit() }
            
            holoColor(B-1,vec(255,0,0))
            Highlighted = B
            timer("unhighlight",2000)
        }
        
        #Set mass
        #Parameters: #body mass
        if(Com[1,string]=="!setmass")
        {
            if(Simulate) { exit() }
            B = Com[2,string]:toNumber()
            if(B>Bodies) { exit() }
            R[B*N-4,number] = Com[3,string]:toNumber()
            print("Mass set.")
        }
        
        #Set size
        #Parameters: #body X Y Z
        if(Com[1,string]=="!setsize")
        {
            B = Com[2,string]:toNumber()
            if(B>Bodies) { exit() }
            X = Com[3,string]:toNumber()
            Y = Com[4,string]:toNumber()
            Z = Com[5,string]:toNumber()
            R[B*N-3,vector] = vec(X,Y,Z)
            holoScaleUnits(B-1,vec(X,Y,Z))
            print("Size set.")
        }
        
        #Set color
        #Parameters: #body R G B
        if(Com[1,string]=="!setcolor")
        {
            B = Com[2,string]:toNumber()
            if(B>Bodies) { exit() }
            X = Com[3,string]:toNumber()
            Y = Com[4,string]:toNumber()
            Z = Com[5,string]:toNumber()
            R[B*N-2,vector] = vec(X,Y,Z)
            holoColor(B-1,vec(X,Y,Z))
            holoColor(B*2+28,vec(X,Y,Z))
            holoColor(B*2+29,vec(X,Y,Z))
            print("Color set.")
        }
        
        #Set position
        #Parameters: #body X Y Z
        if(Com[1,string]=="!setpos")
        {
            if(Simulate) { exit() }
            B = Com[2,string]:toNumber()
            if(B>Bodies) { exit() }
            X = Com[3,string]:toNumber()
            Y = Com[4,string]:toNumber()
            Z = Com[5,string]:toNumber()
            R[B*N-1,vector] = vec(X,Y,Z)
            holoPos(B-1,Center+vec(X,Y,Z))
            holoPos(B*2+28,Center+R[B*N-1,vector]+R[B*N,vector]*5/12)
            holoPos(B*2+29,Center+R[B*N-1,vector]+R[B*N,vector]*5/6)
            print("Position set.")
        }
        
        #Set velocity
        #Parameters: #body X Y Z
        if(Com[1,string]=="!setvel")
        {
            if(Simulate) { exit() }
            B = Com[2,string]:toNumber()
            if(B>Bodies) { exit() }
            X = Com[3,string]:toNumber()
            Y = Com[4,string]:toNumber()
            Z = Com[5,string]:toNumber()
            R[B*N,vector] = vec(X,Y,Z)
            
            holoPos(B*2+28,Center+R[B*N-1,vector]+R[B*N,vector]*5/12)
            holoPos(B*2+29,Center+R[B*N-1,vector]+R[B*N,vector]*21/24)
            A = (quat(R[B*N,vector]:toAngle())*qRotation(vec(0,90,0))):toAngle()
            holoAng(B*2+28,A)
            holoAng(B*2+29,A)
            holoScaleUnits(B*2+28,vec(7,7,R[B*N,vector]:length()*5/6))
            holoScaleUnits(B*2+29,vec(15,15,R[B*N,vector]:length()/6))
            print("Velocity set.")
        }
        
        #Set SLOW rate
        if(Com[1,string]=="!slow")
        {
            if(Simulate) { exit() }
            SLOW = Com[2,string]:toNumber()
            print("SLOW rate set.")
        }
        
        #Set skip rate
        if(Com[1,string]=="!skip")
        {
            if(Simulate) { exit() }
            Skip = Com[2,string]:toNumber()
            print("Skip rate set.")
        }
        
        #Print body data
        if(Com[1,string]=="!data")
        {
            B = Com[2,string]:toNumber()
            if(B>Bodies) { exit() }
            print("Mass: "+R[B*N-4,number]:toString())
            print("Position: "+R[B*N-1,vector]:toString())
            print("Velocity: "+R[B*N,vector]:toString())
        }
    }
}

#------------------------ Simulation --------------------------

if(tickClk())
{
    Center = entity():pos()+vec(0,0,300)
    if(Tick==0)   #tick skipping
    {
        Tick++
        if(Tick>=Skip) {Tick = 0}
        T = curtime()/SLOW

        #Simulation
        if(Simulate)
        {
            P1 = array()
            for(I=1,R:count()/N) {      
                P1:pushVector($T*R[I*N,vector])
                P1V = vec()
                for(J=1,R:count()/N)
                {
                    if(J!=I) {
                        DX1 = R[J*N-1,vector]-R[I*N-1,vector]
                        P1V += $T*G*R[J*N-4,number]/(DX1:length())^3*DX1
                    }
                }
                P1:pushVector(P1V)
            }
    
            P2 = array()
            for(I=1,R:count()/N) { 
                P2:pushVector($T*(R[I*N,vector]+P1[I*2,vector]/2))
                P2V = vec()
                for(J=1,R:count()/N)
                {
                    if(J!=I) {
                        DX2 = (R[J*N-1,vector]+P1[J*2-1,vector]/2)-(R[I*N-1,vector]+P1[I*2-1,vector]/2)
                        P2V += $T*G*R[J*N-4,number]/(DX2:length())^3*DX2
                    }
                }
                P2:pushVector(P2V)
            }
    
            P3 = array()
            for(I=1,R:count()/N) { 
                P3:pushVector($T*(R[I*N,vector]+P2[I*2,vector]/2))
                P3V = vec()
                for(J=1,R:count()/N)
                {
                    if(J!=I) {
                        DX3 = (R[J*N-1,vector]+P2[J*2-1,vector]/2)-(R[I*N-1,vector]+P2[I*2-1,vector]/2)
                        P3V += $T*G*R[J*N-4,number]/(DX3:length())^3*DX3
                    }
                }
                P3:pushVector(P3V)
            }
        
            P4 = array()
            for(I=1,R:count()/N) { 
                P4:pushVector($T*(R[I*N,vector]+P3[I*2,vector]))
                P4V = vec()
                for(J=1,R:count()/N)
                {
                    if(J!=I) {
                        DX4 = (R[J*N-1,vector]+P3[J*2-1,vector])-(R[I*N-1,vector]+P3[I*2-1,vector])
                        P4V += $T*G*R[J*N-4,number]/(DX4:length())^3*DX4
                    }
                }
                P4:pushVector(P4V)
            }
    
            for(I=1,R:count()/N) {    
                PX = (P1[I*2-1,vector]+P2[I*2-1,vector]*2+P3[I*2-1,vector]*2+P4[I*2-1,vector])/6
                PV = (P1[I*2,vector]+P2[I*2,vector]*2+P3[I*2,vector]*2+P4[I*2,vector])/6
        
                R[I*N-1,vector] = R[I*N-1,vector]+PX
                R[I*N,vector] = R[I*N,vector]+PV
        
                holoPos(I-1,Center+R[I*N-1,vector])
            }
        }
    }
    else
    {
        Tick++
        if(Tick>=Skip) { Tick = 0 }
    }
}

if(clk("unhighlight"))
{
    if(Highlighted == -1) { exit() }
    holoColor(Highlighted-1,R[Highlighted*N-2,vector])
    Highlighted = -1
}

Commands:

• !addplayer player - Adds a player to the list of players which can issue commands. You can let your friends modify the system this way, but unfortunately all the info printed by the E2 will be visible only to you.
• !addbody - Adds a body. Body is added at (0,0,0), with velocity (0,0,0), mass 1, size (20,20,20) and white color.
• !setmass B M - Sets mass of body #B to M.
• !setsize B X Y Z - Sets the size of body #B to (X,Y,Z) (you can have elipsoids this way).
• !setcolor B r g b - Sets the color of body #B to (r,g,b)
• !setpos B X Y Z - Sets the position of body #B to (X,Y,Z)
• !setvel B X Y Z - Sets the velocity of body #B to (X,Y,Z) (represented by an arrow )
• !simulate - Resumes/pauses the simulation
• !highlight B - Highlights body #B for 2 seconds (in case you don't remember which is which ).
• !data B - Prints the data (mass, position, velocity) for body #B in the console.
• !slow N - Sets simulation slow-down rate. 1 is normal speed, 2 is 2 times slower, etc.
• !skip N - Sets tick skip rate. Defaults to 1, if you add many bodies increase it, or you will exceed quota limit.

G is 33.167, as in the Solar System simulation. It's set so that 1 mass unit corresponds to 10^24 kg, 1 length unit corresponds to 500 000 km, and 1 second corresponds to 3 months (1/4 of a year). With this information you can see how real-life gravitating systems would behave

I have an idea to add a special body to this, behaving like a spaceship - this way I could make a spaceflight simulation in GMod. It's on my TODO list now, but I don't know when I'll make it.