Gravity Simulator - new beta version

Tony YaBB Administrator Posts: 1054 Gender:	Re: new beta version Reply #30 - 01/25/11 at 18:09:20 Edit the object's name and put an asterisk as the last character. Then you will get red as your negative z color. For example, Ceres* instead of Ceres.
Back to top	IP Logged

phoenixshade

Uploader

Proud to be an Ape

Posts: 26
Gender: male

Re: new beta version
Reply #31 - 01/30/11 at 21:09:57

My wish-list just keeps growing and growing. Almost every time I use this program, I think of one or two minor tweaks I'd like to see. Most of them end up being forgotten, but here's a simple one that wasn't.

Perhaps this would fall under the "Preferences" menu. I'd like to be able to change the value of the gravitational constant G. In the simplest case, it's just a text field that allows you to type in any number you wish, along with a "Default" button. Even better would be to have radio buttons for units that work like the radio buttons for mass units.

Normally, I could achieve the same effect just by changing my masses or distances appropriately (for example, if I want to double G, I could just multiply all my masses by √2) except in one case: Changing the sign of G. This would make gravity act something like electricity, with opposite-signed masses attracting and like-signed ones repelling. If I only wanted two-body versions, again no problem — just make one of the masses negative — but this won't work for n-body problems.

One of these days I might just break down and write my own simulator. Not that there's anything wrong with this one; I'm just a tinkerer by nature.

Back to top

IP Logged

EDG Ultimate Member oh, crumbs!!! Posts: 611 Gender:	Re: new beta version Reply #32 - 01/30/11 at 22:11:15 Huh. That'd be kinda interesting
Back to top	(formerly known as Mal) IP Logged

frankuitaalst

Ultimate Member

Great site

Posts: 1507
Gender: male

Re: new beta version
Reply #33 - 01/31/11 at 10:11:30

Quote from phoenixshade on 01/30/11 at 21:09:57:

Normally, I could achieve the same effect just by changing my masses or distances appropriately (for example, if I want to double G, I could just multiply all my masses by √2) except in one case: Changing the sign of G. This would make gravity act something like electricity, with opposite-signed masses attracting and like-signed ones repelling. If I only wanted two-body versions, again no problem — just make one of the masses negative — but this won't work for n-body problems.
.

Be prepared to see some unexpected results if you're playing with negative G or negative masses . It may be good to define two kinds of mass in his case : inertial and gravitational mass . I guess Gravsim has in the .gsim files some space left to do so .

Back to top

IP Logged

Tony

YaBB Administrator

Posts: 1054
Gender: male

Re: new beta version
Reply #34 - 02/01/11 at 10:57:26

There's still lots of rffu's (reserved for future use) in the gsim file. That's something I didn't do with Gravity Simulator 1.0, and adding more stuff to save became a huge chore.

At the moment, I'm not able to compile C++, so any tinkering with G would have to be executed in the VB code. That means ignoring codeoptomizer.dll. The sims will run slower, but depending on what you're trying to demonstrate, that might not matter.

Back to top

IP Logged

Tony YaBB Administrator Posts: 1054 Gender:	Re: new beta version Reply #35 - 02/01/11 at 15:56:51 http://orbitsimulator.com/gravity/beta/GravitySimulatorBeta2011Feb01.exe This has a Preferences option for changing G. Gravity Simulator ignores the dll if a custom value is chosen. Let me know if it's buggy.
Back to top	IP Logged

Tony

YaBB Administrator

Posts: 1054
Gender: male

Re: new beta version
Reply #36 - 05/30/11 at 12:47:46

http://orbitsimulator.com/gravity/beta/GravitySimulatorBeta2011May30.exe

This beta introduces a 4th order Runge Kutta engine to Gravity Simulator. In the Beta menu you can choose between the two. I'm working on an analysis to compare the new RK4 method with the existing Euler-Cromer method. They both have their strengths and weaknesses. After minimal testing, it seems that RK4 is the best choice for short-duration simulations that use small time steps, for example, a recreation of Apollo missions. And it seems like EC is the best choice for simulations spanning years into the future, and/or using large time steps. I hope everyone will play around with it and let me know what you think.

This beta does not use the external dll, so the integrators are done in Visual Basic rather than C++. Later I'll make a new C++ dll for this one.

Back to top

IP Logged

Tony

YaBB Administrator

Posts: 1054
Gender: male

Re: new beta version
Reply #37 - 06/01/11 at 17:51:58

http://orbitsimulator.com/gravity/beta/GravitySimulatorBeta2011June01.exe
This fixes a bug in the new 4th order Runge Kutta routine. This new routine now passes every test I've thrown at it. When run at the same time step as the Euler-Cromer routine, it is a bit slower, but to make up for that, you can take much larger time steps without compromising accuracy. At the moment, collisions don't work with the RK4 method.

Back to top

IP Logged

Tony

YaBB Administrator

Posts: 1054
Gender: male

Re: new beta version
Reply #38 - 06/01/11 at 17:52:30

Code:

'=======================  RK4 ==========================================
'Use current positions to compute accelerations.  Store them in arrays a1x(), a1y(), a1z().
	
    'zero out arrays a1x(), a1y(), a1z()
    For k = 1 To objects
	  a1x(k) = 0: a1y(k) = 0: a1z(k) = 0
    Next k
	
    For k = 1 To objects
	  Mu = ObjMass(k) * G
	  For j = k + 1 To objects
		Mu2 = ObjMass(j) * G
		dx = (objx(j) - objx(k))
		dy = (objy(j) - objy(k))
		dz = (objz(j) - objz(k))
		D = Sqr(dx * dx + dy * dy + dz * dz)
		If Mu > 0 Then 'ignore massless particles
		    a1 = -(1 / D) * (1 / D) * Mu
		    a1x(j) = a1x(j) + (dx / D) * a1
		    a1y(j) = a1y(j) + (dy / D) * a1
		    a1z(j) = a1z(j) + (dz / D) * a1
		End If
		
		If Mu2 > 0 Then 'ignore massless particles
		    a1 = -(1 / D) * (1 / D) * Mu2
		    a1x(k) = a1x(k) - (dx / D) * a1
		    a1y(k) = a1y(k) - (dy / D) * a1
		    a1z(k) = a1z(k) - (dz / D) * a1
		End If
	  Next j
    Next k

'Copy the current velocities into an arrays v1x(), v1y(), v1z().
    For k = 1 To objects
	  v1x(k) = objvx(k)
	  v1y(k) = objvy(k)
	  v1z(k) = objvz(k)
    Next k

'Use the velocities stored in v1x(), v1y(), v1z() to predict positions of objects one-half a step into the future.
'Store the positions in arrays p2x(), p2y(), p2z()."
    halfstep = timescale / 2
    For k = 1 To objects
	  p2x(k) = objx(k) + v1x(k) * halfstep
	  p2y(k) = objy(k) + v1y(k) * halfstep
	  p2z(k) = objz(k) + v1z(k) * halfstep
    Next k
   
'Use the accelerations stored in a1x(), a1y(), a1z to predict future velocities at half a step into the future.
'Store the velocities in arrays v2x(), v2y(), v2z().
    For k = 1 To objects
	  v2x(k) = objvx(k) + a1x(k) * halfstep
	  v2y(k) = objvy(k) + a1y(k) * halfstep
	  v2z(k) = objvz(k) + a1z(k) * halfstep
    Next k

'Use the positions stored in arrays p2x(), p2y(), p2z() (these are the positions at half a step in future)
'to compute accelerations on the bodies a half timestep into the future.
'Store these accelerations in arrays a2x(), a2y(), a2z()

    'zero out arrays a2x(), a2y(), a2z()
    For k = 1 To objects
	  a2x(k) = 0: a2y(k) = 0: a2z(k) = 0
    Next k
	
    For k = 1 To objects
	  Mu = ObjMass(k) * G
	  For j = k + 1 To objects
		Mu2 = ObjMass(j) * G
		dx = (p2x(j) - p2x(k))
		dy = (p2y(j) - p2y(k))
		dz = (p2z(j) - p2z(k))
		D = Sqr(dx * dx + dy * dy + dz * dz)
		If Mu > 0 Then 'ignore massless particles
		    a2 = -(1 / D) * (1 / D) * Mu
		    a2x(j) = a2x(j) + (dx / D) * a2
		    a2y(j) = a2y(j) + (dy / D) * a2
		    a2z(j) = a2z(j) + (dz / D) * a2
		End If
		
		If Mu2 > 0 Then 'ignore massless particles
		    a2 = -(1 / D) * (1 / D) * Mu2
		    a2x(k) = a2x(k) - (dx / D) * a2
		    a2y(k) = a2y(k) - (dy / D) * a2
		    a2z(k) = a2z(k) - (dz / D) * a2
		End If
	  Next j
    Next k
   
'Use the accelerations stored in arrays a2x(), a2y(), a2z() to compute velocities at half a step into the future.
'Store these velocities in arrays v3x(), v3y(), v3z().
    For k = 1 To objects
	  v3x(k) = objvx(k) + a2x(k) * halfstep
	  v3y(k) = objvy(k) + a2y(k) * halfstep
	  v3z(k) = objvz(k) + a2z(k) * halfstep
    Next k

'Use the velocities stored in arrays v2x(), v2y(), v2z() to predict positions of objects one-half a step into the future.
'Store these positions in arrays p3x(), p3y(), p3z().
    For k = 1 To objects
	  p3x(k) = objx(k) + v2x(k) * halfstep
	  p3y(k) = objy(k) + v2y(k) * halfstep
	  p3z(k) = objz(k) + v2z(k) * halfstep
    Next k

'Use the positions in arrays p3x(), p3y(), p3z() (these positions are at half a step in future)
'to compute the accelerations on the bodies a half time step into the future.
'Store them in arrays a3x(), a3y(), a3z()
    For k = 1 To objects
	  a3x(k) = 0: a3y(k) = 0: a3z(k) = 0
    Next k
	
    For k = 1 To objects
	  Mu = ObjMass(k) * G
	  For j = k + 1 To objects
		Mu2 = ObjMass(j) * G
		dx = (p3x(j) - p3x(k))
		dy = (p3y(j) - p3y(k))
		dz = (p3z(j) - p3z(k))
		D = Sqr(dx * dx + dy * dy + dz * dz)
		If Mu > 0 Then 'ignore massless particles
		    a3 = -(1 / D) * (1 / D) * Mu
		    a3x(j) = a3x(j) + (dx / D) * a3
		    a3y(j) = a3y(j) + (dy / D) * a3
		    a3z(j) = a3z(j) + (dz / D) * a3
		End If
		
		If Mu2 > 0 Then 'ignore massless particles
		    a3 = -(1 / D) * (1 / D) * Mu2
		    a3x(k) = a3x(k) - (dx / D) * a3
		    a3y(k) = a3y(k) - (dy / D) * a3
		    a3z(k) = a3z(k) - (dz / D) * a3
		End If
	  Next j
    Next k

'Use the accelerations in arrays a3x(), a3y(), a3z() to compute velocities at one FULL step into the future.
'Store them in arrays v4x(), v4y(), v4z().
    For k = 1 To objects
	  v4x(k) = objvx(k) + a3x(k) * timescale
	  v4y(k) = objvy(k) + a3y(k) * timescale
	  v4z(k) = objvz(k) + a3z(k) * timescale
    Next k

'Use the veolcities stored in arrays v3x(), v3y(), v3z() to compute the positions of objects one FULL step into the future.
'Store them in arrays p4x(), p4y(), p4z().
    For k = 1 To objects
	  p4x(k) = objx(k) + v3x(k) * timescale
	  p4y(k) = objy(k) + v3y(k) * timescale
	  p4z(k) = objz(k) + v3z(k) * timescale
    Next k
	
'Use the positions stored in arrays p4x(), p4y(), p4z()(these are the positions at one FULL step in future)
'to compute accelerations on the bodies in the future.
'Store these accelerations in arrays a4x(), a4y(), a4z()
    For k = 1 To objects
	  a4x(k) = 0: a4y(k) = 0: a4z(k) = 0
    Next k
	
    For k = 1 To objects
	  Mu = ObjMass(k) * G
	  For j = k + 1 To objects
		Mu2 = ObjMass(j) * G
		dx = (p4x(j) - p4x(k))
		dy = (p4y(j) - p4y(k))
		dz = (p4z(j) - p4z(k))
		D = Sqr(dx * dx + dy * dy + dz * dz)
		If Mu > 0 Then 'ignore massless particles
		    a4 = -(1 / D) * (1 / D) * Mu
		    a4x(j) = a4x(j) + (dx / D) * a4
		    a4y(j) = a4y(j) + (dy / D) * a4
		    a4z(j) = a4z(j) + (dz / D) * a4
		End If
		
		If Mu2 > 0 Then 'ignore massless particles
		    a4 = -(1 / D) * (1 / D) * Mu2
		    a4x(k) = a4x(k) - (dx / D) * a4
		    a4y(k) = a4y(k) - (dy / D) * a4
		    a4z(k) = a4z(k) - (dz / D) * a4
		End If
	  Next j
    Next k


'Apply the Runge-Kutta 4th Order method to these arrays to determine the slope of the position and velocities
'Then add these slopes to the current positions and velocities to get the new positions and velocities.
'Velocity is the derivative of position.
'Acceleration is the derivative of velocity.

'The RK4 method states that the slope for velocity is
'(1/6) * timestep * (the acceleration at the beginning of the interval, stored in a1x(), a1y(), a1z()
'			  + 2 * the acceleration at the first estimate of the half time step, stored in a2x(), a2y(), a2z()
'			  + 2 * the acceleration at the second estimate of the half time step, stored in a3x(), a3y(), a3z()
'			  + the acceleration at the estimate of the end of a full time step, stored in a4x(), a4y(), a4z() ).
	
'The RK4 method states that the slope for position is
'(1/6) * timestep * (the velocity at the beginning of the interval, stored in v1x(), v1y(), v1z()
'			  + 2 * the velocity at the first estimate of the half time step, stored in v2x(), v2y(), v2z()
'			  + 2 * the velocity at the second estimate of the half time step, stored in v3x(), v3y(), v3z()
'			  + the velocity at the estimate of the end of a full time step, stored in v4x(), v4y(), v4z() ).

	
    For k = 1 To objects
	  deltavx = (timescale / 6) * (a1x(k) + 2 * a2x(k) + 2 * a3x(k) + a4x(k))
	  deltavy = (timescale / 6) * (a1y(k) + 2 * a2y(k) + 2 * a3y(k) + a4y(k))
	  deltavz = (timescale / 6) * (a1z(k) + 2 * a2z(k) + 2 * a3z(k) + a4z(k))
	  
	  objvx(k) = objvx(k) + deltavx
	  objvy(k) = objvy(k) + deltavy
	  objvz(k) = objvz(k) + deltavz
	  
	  deltapx = (timescale / 6) * (v1x(k) + 2 * v2x(k) + 2 * v3x(k) + v4x(k))
	  deltapy = (timescale / 6) * (v1y(k) + 2 * v2y(k) + 2 * v3y(k) + v4y(k))
	  deltapz = (timescale / 6) * (v1z(k) + 2 * v2z(k) + 2 * v3z(k) + v4z(k))
	  
	  objx(k) = objx(k) + deltapx
	  objy(k) = objy(k) + deltapy
	  objz(k) = objz(k) + deltapz
    Next k

Back to top

IP Logged

Tony

YaBB Administrator

Posts: 1054
Gender: male

Re: new beta version
Reply #39 - 06/08/11 at 20:34:55

http://orbitsimulator.com/gravity/beta/GravitySimulatorBeta2011June08.exe

It seems my last beta version wiped out the rotating frames. That's fixed now. The "distance and velocity box" now shows more digits, and the radial velocity reads negative if the objects are approaching each other.

This version also introduces two new integration methods: Leapfrog and Verlet. These are 2nd order routines whose accuracy rivals RK4, but whose execution speed rivals Euler-Cromer. You can choose between them in the new Integrator menu.

Back to top

IP Logged

abyssoft YaBB Administrator I love YaBB 1G - SP1! Posts: 302	Re: new beta version Reply #40 - 06/09/11 at 15:18:03 Create object appears to have broken with the latest beta. I do create object and it doesn't adhere to the parameters I end up with stuff scattered all over the place, have a system i can occasionally use to do runs 8)
Back to top	IP Logged

Tony

YaBB Administrator

Posts: 1054
Gender: male

Re: new beta version
Reply #41 - 06/09/11 at 17:39:48

Quote from abyssoft on 06/09/11 at 15:18:03:

Create object appears to have broken with the latest beta.

I do create object and it doesn't adhere to the parameters I end up with stuff scattered all over the place, have a system i can occasionally use to do runs 8)

I just tried it and it works fine for me. Can you give me some details of what you're creating and the steps you're using?

Back to top

IP Logged

abyssoft YaBB Administrator I love YaBB 1G - SP1! Posts: 302	Re: new beta version Reply #42 - 06/11/11 at 00:00:51 starting with planets only sim paused tried all 4 engines creating 600 obs, 1/2400 mass earth, earth as central obj, +/- % at 25 for all variances, inc 5, SMA 80000, radius 250 km, orbital params un modified 180 mod 100%,
Back to top	IP Logged

Tony

YaBB Administrator

Posts: 1054
Gender: male

Re: new beta version
Reply #43 - 06/11/11 at 10:53:19

Thanks for the report.
You'll probably find that there is a problem with onlyplanets.gsim even if you don't add more items. There are some problems opening simulations created with previous versions that are saved with time steps other than 1 second in the new beta. This might sound like a pain to try, but open onlyplanets in an older beta (prior to May 2011), reduce the time step to 1, resave it. Then it should open fine in the new beta and you can create your additional objects and it should work fine. I'll correct this soon so you can open older simulations with the Beta, but for the time being, you can only open existing sims saved at timestep = 1.

Back to top

IP Logged

wil

Uploader

I Love YaBB 2!

Posts: 39

Re: new beta version
Reply #44 - 06/12/11 at 13:39:08

I propose new option:
show orbital inclinations and orientations (line of nodes) in invariant plane of system.

Invariant plane is perpendicular to total rotational momentum vector.
Simple sum: L = m r x v = [Lx, Ly,Lz].
Normalize it and we have normal vector of invariant plane.
----

And other one yet:

8 - 10-th order integrator.
For example RK10 will be good enough, mayby some symplectic version...

Back to top

« Last Edit: 06/14/11 at 13:37:33 by wil »

IP Logged