Runge-Kutta 4 E2 (Gravity Simulator)

[Fizyk]

I got the idea from this thread

This E2 uses Runge-Kutta integration method to simulate the gravity in a multiple body system. Basically that means it is pretty stable

To test it, I made a test E2 that simulates a Sun-Earth-Moon-like system. Here is the code:

@name Runge-Kutta Test #2
@inputs 
@outputs R:array
@persist T Tick
#constants
@persist G N Skip
@trigger 

if(first()|duped()) {
    runOnTick(1)
    Skip = 3
    Tick = 1
	
    G = 6.67259e-10
    R:clear()
    N=5 #5 entries per body
    
    #"Sun"
    R[1,number]=1.989e+16
    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,-3.6,0) #initial velocity
    
    #"Earth"
    R[6,number]=5.975e+14
    R[7,vector]=vec(20,20,20)
    R[8,vector]=vec(0,255,255)
    R[9,vector]=vec(1000,0,0)
    R[10,vector]=vec(0,120,0)
    
    #"Moon"
    R[11,number]=7e12
    R[12,vector]=vec(10,10,10)
    R[13,vector]=vec(100,100,100)
    R[14,vector]=vec(1050,0,0)
    R[15,vector]=vec(0,195,30)
    
    #"Mercury"
    R[16,number]=7e10
    R[17,vector]=vec(8,8,8)
    R[18,vector]=vec(100,100,0)
    R[19,vector]=vec(0,200,0)
    R[20,vector]=vec(-260,0,40)
    
    #"Black Hole"
    #R[21,number]=1e19
    #R[22,vector]=vec(120,120,120)
    #R[23,vector]=vec(0,0,0)
    #R[24,vector]=vec(0,-10000,10000)
    #R[25,vector]=vec(0,5,5)
    
    T = curtime()
    
    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])
        holoPos(I,entity():pos()+vec(0,0,300)+R[I*N-1,vector])
        holoEntity(I):setTrails(30,30,2,"trails/laser",R[I*N-2,vector],500)
    }
}
elseif(Tick == 0) {
    Tick++
    if(Tick>=Skip) {Tick = 0}
	
#Simulation
    T = curtime()
    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 }
}

The RK4 simulation part starts at the "#Simulation" comment. It's possible that it could be optimized (multiple for-loops nested in a for-loop), but I'm too sleepy to try it today

To define the bodies, use the array initialization part. Each body takes 5 entries - mass, size, color, initial position, initial velocity. Just change the entries in the array and have fun But beware of high ops usage - current 3-body setup eats 3000 ops.

You can also try messing with the gravitational constant G to see what it changes

I may post some pics tomorrow, but no picture will show you what fun it is to change masses, positions and velocities and accidentally put Earth on a hyperbolic trajectory :P

EDIT: I posted updated code. I fixed the RK4 method a bit, it had some errors before. I also changed runOnTick to interval, which reduced the ops usage almost 4 times, while being still quite precise. It can handle up to 5 bodies now (maybe 6, but with 5 it eats 2600 ops, adding one more body might be too much).

Also, I took some pics (haven't had time for the video, sorry):





You can see the "Sun", "Earth" and "Moon" there.

Also, I left a "Black Hole" part in the code - if you uncomment it, you can see what a very massive body can do to the system, even if it's far away

EDIT2: Syranide suggested that I still use runOnTick but just skip some executions. I modified the code to do that, but I didn't test it, so it may fail - if it does, please tell me. If it works, you can adjust how many executions it should skip by changing the Skip variable in the initialization part - Skip=1 means no skipping, Skip=2 makes it run each 2 executions, etc.

EDIT3: I did some more testing. Skipping ticks works and 6 bodies are doable - with skip 3 it goes up to 5055 ops sometimes, but with skip=4, 6 bodies are no problem

[hzzzln]

somehow, your avatar fits the things you make. i know what they are, what they do, but i fucking dont know how.

this is awesome. (the e2, not the fact that i dont understand it.)

[d3cr1pt0r]

I actually started gmod after a long time, just to see this in action. Awsome work!

[SystemsLock]

I am very impressed! It seems like you took my idea and made it far better. I'd love to see how it would react to larger systems. Perhaps non predetermined ones.

I'm working on my own version with holos as well but it looks like I'll be playing catch up.

[Fizyk]

Thanks I just realized I made one thing a bit wrong and it may cause small loss of precision, so I'm going to rewrite this E2 today and make a video of it

[Fizyk]

Update bump: corrected code, also it is now capable of handling up to 5, maybe 6 bodies (I reduced ops by making it run less often ). I also added some pics, see the first post.

EDIT: Modified it a bit more, details are in EDIT2 and EDIT3 in the first post

[SystemsLock]

I'm still working on my own version, but your is coming along nicely. Have your tried a multi solar orbit yet?

I'm pretty sure you noticed this buy your code could be shrunk, a lot, if you used a nested loop.

You should really label you variables, it's a bad habit.

[Fizyk]

Have your tried a multi solar orbit yet?

Nope, but it should work. It could be a problem to find a stable orbit though.

I'm pretty sure you noticed this buy your code could be shrunk, a lot, if you used a nested loop.

Unless you can make an array of arrays, no, it can't. Or maybe I don't see how, can you elaborate?

You should really label you variables, it's a bad habit.

If you mean lack of comments, then yeah, I probably could comment my code more.
And if you mean names of the variables, the ones that could be descriptive without making my code really huge already are descriptive.

[SystemsLock]

Alright well the whole lower half of your code looks like gibberish to me. But from what I can gather you repeat a lot of it. All I know is that my code doesn't do that, has a nested loop, and is much smaller. It's not finished but I'll post it on completion.

[Fizyk]

You're right, it is repeated in a way, but unfortunately that's a piece of code which can't be put in a loop currently. This is exactly the part where I do the Runge-Kutta method. It requires me to perform 4 steps on the positions and velocities, and the result is 4 arrays, P1 - P4, which contain some information about changes in positions and velocities of all the bodies. I need all 4 of them, to calculate PX and PV (overall change in positions and velocities), so unless I can make an array of 4 arrays, there is no way I can loop through them.

As I think now, I could change names of PX and PV, they aren't really the best names possible, but it is more consistent with P1-P4 (which are just what these values were called in a lecture about RK4)

EDIT: I just realized that with some workarounds I could put that in a loop, but it would require quite a bit of work and would probably eat all the remaining ops, so I prefer to leave it as it is