Results 1 to 8 of 8

Thread: RK4 Solar System

  1. #1
    Wire Sofaking Fizyk's Avatar
    Join Date
    Jun 2008
    Location
    Łomianki, Poland
    Posts
    738
    Blog Entries
    1

    Default RK4 Solar System

    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:

    Code:
    @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!
    Attached Thumbnails Attached Thumbnails RK4 Solar System-gm_buildgrass0003.jpg   RK4 Solar System-gm_buildgrass0002.jpg   RK4 Solar System-gm_buildgrass0001.jpg  
    Last edited by Fizyk; 05-18-2010 at 08:54 AM. Reason: Added screenshots

    My programs: BIOS - Alcyone - Calculator - Notepad - Movie Player
    My tutorials: applyTorque - Quaternions - PID controllers
    Some other things I made: FT Chip - RK4 Solar System

  2. #2
    Wirererer Manslaughter's Avatar
    Join Date
    Sep 2009
    Location
    Eureka, NT: owner():pos()
    Posts
    150

    Default Re: RK4 Solar System

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

    Newton and Galileo would be prowd.
    Use Dropbox!
    It's Easy, Free, File Sharing

    I'm tired of these ad ridden junk sites, this is free, no ad's, 8!GB Space.
    Build, Fight, Survive
    My FTS Server: 98.237.13.239:27015
    Quote Originally Posted by Albert Einstein
    Two things are infinite: the universe and human stupidity; and I'm not sure about the the universe.

  3. #3
    Wirererer Paper Clip's Avatar
    Join Date
    Sep 2008
    Location
    Canada
    Posts
    368

    Default Re: RK4 Solar System

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

  4. #4
    Wire Sofaking Fizyk's Avatar
    Join Date
    Jun 2008
    Location
    Łomianki, Poland
    Posts
    738
    Blog Entries
    1

    Default Re: RK4 Solar System

    Thanks

    My programs: BIOS - Alcyone - Calculator - Notepad - Movie Player
    My tutorials: applyTorque - Quaternions - PID controllers
    Some other things I made: FT Chip - RK4 Solar System

  5. #5
    Wire Sofaking SystemsLock's Avatar
    Join Date
    Mar 2009
    Posts
    474

    Default Re: RK4 Solar System

    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.
    Make a Small Loan, Make a Big Difference - Check out Kiva.org to Learn How!

  6. #6
    Wirererer Paper Clip's Avatar
    Join Date
    Sep 2008
    Location
    Canada
    Posts
    368

    Default Re: RK4 Solar System

    Quote Originally Posted by SystemsLock View Post
    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

  7. #7
    Wire Amateur sebi99's Avatar
    Join Date
    Dec 2009
    Location
    ring0
    Posts
    33

    Default Re: RK4 Solar System

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

  8. #8
    Wire Sofaking Fizyk's Avatar
    Join Date
    Jun 2008
    Location
    Łomianki, Poland
    Posts
    738
    Blog Entries
    1

    Default Re: RK4 Solar System

    Quote Originally Posted by SystemsLock
    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

    Quote Originally Posted by Paper Clip
    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

    Quote Originally Posted by sebi99
    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.
    Code:
    @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.
    Last edited by Fizyk; 05-18-2010 at 11:36 AM.

    My programs: BIOS - Alcyone - Calculator - Notepad - Movie Player
    My tutorials: applyTorque - Quaternions - PID controllers
    Some other things I made: FT Chip - RK4 Solar System

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •  
unnecessary
unnecessary
unnecessary
unnecessary
linguistic-parrots
linguistic-parrots
linguistic-parrots
linguistic-parrots