Eccentric Orbits using Quaternions

This code has been declared by its author to be Public Domain code.

Eccentric Orbits using Quaternions by Krischan

2017

This code demonstrates how to create a star system with eccentric orbits using Quaternion rotations in miniB3D. The line coordinates for each orbit are stored in a separate Type (TOrbit). The planet pivot is then turned one step further and this goes on and on until a full 360° cycle has been finished.

It uses adaptive segment detail, so closer planet orbits use less segments while planets far away use more segments to have a constant orbit line detail over distance (otherwise, far away planets wouldn't match the orbit line and appear beside or above/below it). Use the mouse to view the scene from different angles (Quarternion rotated, too!) and the Mousewheel to zoom in/out the center.

eccentricity = 0.0

eccentricity = 45.0

SuperStrict

Framework sidesign.minib3d

Graphics3D DesktopWidth(), DesktopHeight(), 32, 2

SeedRnd 1

Const eccentricity:Float = 45		' angle between 0 and 90
Const planets:Int = 10				' number of planets
Const speed:Float = 0.1				' simulation speed
Const maxsegments:Int = 1440		' maximum segments of outer planet orbit

Global MouseWheel:Int = 1

' ------------------------------------------------------------------------------------------------
' Initialize Scene and Camera
' ------------------------------------------------------------------------------------------------
Global pivot:TPivot = CreatePivot()

Global cam:TCamera = CreateCamera(pivot)
PositionEntity cam, 0, 5, -15

Global star:TMesh = CreateSphere(32)
ScaleEntity star, 0.5, 0.5, 0.5
EntityColor star, 255, 192, 0
EntityFX star, 1

PointEntity cam, star

' ------------------------------------------------------------------------------------------------
' Stores Planet properties
' ------------------------------------------------------------------------------------------------
Type TPlanet

	Field id:Int
	Field rotation:Float[]
	Field ent:TEntity
	Field piv:TPivot
	Field radius:Float
	Field rgb:Int[] = [255, 255, 255]
	
	Field orbit:TOrbit = New TOrbit
	
End Type

' ------------------------------------------------------------------------------------------------
' Stores Orbit Coordinates
' ------------------------------------------------------------------------------------------------
Type TOrbit

	Field steps:Float
	Field x:Float[maxsegments + 1]
	Field y:Float[maxsegments + 1]
	Field z:Float[maxsegments + 1]

End Type

' ------------------------------------------------------------------------------------------------
' Initialize Planets
' ------------------------------------------------------------------------------------------------
Global planetlist:TList = CreateList()
Global p:TPlanet

For Local i:Int = 1 To planets
	
	p = New TPlanet
	p.id = i
	p.radius = Normalize(i, 1, planets, 1, 10)
	p.rotation = [Float(Rnd(-eccentricity, eccentricity)), Float(Rnd(360)), 0.0]
	p.rgb = [Rand(255), Rand(255), Rand(255)]
	
	p.orbit.steps = maxsegments / (planets + 1 - p.id)
	
	' planet pivot
	p.piv = CreatePivot()
	
	' create planet
	p.ent = CreateSphere(8, p.piv)
	EntityFX p.ent, 1
	'EntityColor p.ent, p.rgb[0], p.rgb[1], p.rgb[2]
	ScaleEntity p.ent, 0.1, 0.1, 0.1
	
	' pivot initial rotation
	RotateEntity p.piv, 0, 0, 0
	Turn p.piv, p.rotation[0], p.rotation[1], p.rotation[2]
	
	' position planet
	PositionEntity p.ent, 0, 0, 0
	MoveEntity p.ent, 0, 0, p.radius
	
	ListAddLast planetlist, p

	' store current position in Type for later use
	For Local i:Int = 0 To p.orbit.steps
	
		p.orbit.x[i] = EntityX(p.ent, 1)
		p.orbit.y[i] = EntityY(p.ent, 1)
		p.orbit.z[i] = EntityZ(p.ent, 1)
		
		' rotate pivot
		Turn p.piv, 0, 360.0 / p.orbit.steps, 0
				
	Next
		
Next

MoveMouse GraphicsWidth() / 2, GraphicsHeight() / 2

' Mousewheel logic
AddHook EmitEventHook, EventHook

' ------------------------------------------------------------------------------------------------
' Main Loop
' ------------------------------------------------------------------------------------------------
While Not AppTerminate()

	If KeyHit(KEY_ESCAPE) Then End
	
	Turn pivot, Normalize(MouseY(), 0, GraphicsHeight() - 1, -1, 1), Normalize(MouseX(), 0, GraphicsWidth() - 1, -1, 1), 0
				
	CameraZoom cam, MouseWheel
	
	RenderWorld
	
	' enable OpenGL Lines
	' it is IMPORTANT to use this between Renderworld and BeginMax2D or it won't work
	glEnable(GL_COLOR_MATERIAL)
	glBegin(GL_LINES)
	
	Local x1:Float, y1:Float, z1:Float
	Local x2:Float, y2:Float, z2:Float
	
	' cycle all planets
	For Local p:TPlanet = EachIn planetlist
			
		' a full rotation
		For Local i:Int = 0 To p.orbit.steps
		
			' look ahead for next XYZ position
			Local j:Int = (i + 1) Mod p.orbit.steps
			
			' get current coordinates
			x1 = p.orbit.x[i]
			y1 = p.orbit.y[i]
			z1 = p.orbit.z[i]
						
			' get next coordinates
			x2 = p.orbit.x[j]
			y2 = p.orbit.y[j]
			z2 = p.orbit.z[j]
					
			' draw 3D line 
			Line3D(x1, y1, z1, x2, y2, z2, p.rgb[0], p.rgb[1], p.rgb[2], 1)

		Next

		' turn pivot = move planet around star
		Turn p.piv, 0, (planets + 1 - p.radius) * speed, 0
		
		' patch to fix the "running out of sync" bug due Quaternion imprecision
		RotateEntity p.piv, EntityPitch(p.piv), EntityYaw(p.piv), EntityRoll(p.piv)
				
	Next
	
	glEnd
	glDisable(GL_COLOR_MATERIAL)
	
	BeginMax2D()
	
		DrawText "Maximum Eccentricity: " + eccentricity, 0, 0
	
	EndMax2D()
	
	Flip True

Wend

End

' ------------------------------------------------------------------------------------------------
' Normalizes a value to given range
' ------------------------------------------------------------------------------------------------
Function Normalize:Float(Value:Float, vmin:Float, vmax:Float, nmin:Float, nmax:Float, limit:Int = False)

	' normalize	
	Local result:Float = ((Value - vmin) / (vmax - vmin)) * (nmax - nmin) + nmin

	' limit
	If limit Then
	
		If Value >= nmax Then result = nmax
		If Value <= nmin Then result = nmin

	EndIf

	Return result
	
End Function

' ------------------------------------------------------------------------------------------------
' Draw a 3D line using OpenGL functions:
' glEnable(GL_COLOR_MATERIAL) / glBegin(GL_LINES) [...] glEnd() / glDisable(GL_COLOR_MATERIAL)
' ------------------------------------------------------------------------------------------------
Function Line3D:TMesh(x0:Float, y0:Float, z0:Float, x1:Float, y1:Float, z1:Float, r:Int = 255, g:Int = 255, b:Int = 255, a:Float = 1.0)
		
	SetColor r, g, b
	SetAlpha(a)
		
	glVertex3f(x0, y0, -z0)
	glVertex3f(x1, y1, -z1)
	
	SetAlpha(1.0)
			
End Function

' ------------------------------------------------------------------------------------------------
' Turn Entity using Quaternions
' ------------------------------------------------------------------------------------------------
Function Turn(Ent:TEntity, X:Float, Y:Float, Z:Float, Glob:Int = False)

	Local Pitch:Float = 0.0
	Local Yaw:Float = 0.0
	Local Roll:Float = 0.0
	
	Local Quat:TQuaternion = EulerToQuat(0.0, 0.0, 0.0)
	Local Turn_Quat:TQuaternion = EulerToQuat(0.0, 0.0, 0.0)
	
	If Glob = False
	
		Quat = EulerToQuat(EntityPitch(Ent, True), EntityYaw(Ent, True), EntityRoll(Ent, True))
		Turn_Quat = EulerToQuat(X, Y, Z)
		Quat = MultiplyQuats(Quat, Turn_Quat)
		Quat = NormalizeQuat(Quat)
		QuatToEuler2(Quat.x, Quat.y, Quat.z, Quat.w, Pitch, Yaw, Roll)
		RotateEntity Ent, Pitch, Yaw, Roll
	Else
	
		RotateEntity Ent, EntityPitch(Ent) + X, EntityYaw(Ent) + Y, EntityRoll(Ent) + Z
		
	EndIf
	
End Function

' ------------------------------------------------------------------------------------------------
' Euler to Quaternion
' ------------------------------------------------------------------------------------------------
Function EulerToQuat:TQuaternion(pitch:Float, yaw:Float, roll:Float)

	Local cr:Float = Cos(-roll / 2.0)
	Local cp:Float = Cos(pitch / 2.0)
	Local cy:Float = Cos(yaw / 2.0)
	Local sr:Float = Sin(-roll / 2.0)
	Local sp:Float = Sin(pitch / 2.0)
	Local sy:Float = Sin(yaw / 2.0)
	Local cpcy:Float = cp * cy
	Local spsy:Float = sp * sy
	Local spcy:Float = sp * cy
	Local cpsy:Float = cp * sy
	
	Local q:TQuaternion = New TQuaternion
	
	q.w:Float = cr * cpcy + sr * spsy
	q.x:Float = sr * cpcy - cr * spsy
	q.y:Float = cr * spcy + sr * cpsy
	q.z:Float = cr * cpsy - sr * spcy
	
	Return q
	
End Function

' ------------------------------------------------------------------------------------------------
' Quaternion to Euler
' ------------------------------------------------------------------------------------------------
Function QuatToEuler2(x:Float, y:Float, z:Float, w:Float, pitch:Float Var, yaw:Float Var, roll:Float Var)

	Local QuatToEulerAccuracy:Double = 1.0 / 2 ^ 31

	Local sint:Float = (2.0 * w * y) - (2.0 * x * z)
	Local cost_temp:Float = 1.0 - (sint * sint)
	Local cost:Float

	If Abs(cost_temp) > QuatToEulerAccuracy

		cost = Sqr(cost_temp)

	Else

		cost = 0.0

	EndIf

	Local sinv:Float, cosv:Float, sinf:Float, cosf:Float
	
	If Abs(cost) > QuatToEulerAccuracy
	
		sinv = ((2.0 * y * z) + (2.0 * w * x)) / cost
		cosv = (1.0 - (2.0 * x * x) - (2.0 * y * y)) / cost
		sinf = ((2.0 * x * y) + (2.0 * w * z)) / cost
		cosf = (1.0 - (2.0 * y * y) - (2.0 * z * z)) / cost
		
	Else
		
		sinv = (2.0 * w * x) - (2.0 * y * z)
		cosv = 1.0 - (2.0 * x * x) - (2.0 * z * z)
		sinf = 0.0
		cosf = 1.0
		
	EndIf
	
	pitch = ATan2(sint, cost)
	yaw = ATan2(sinf, cosf)
	roll = -ATan2(sinv, cosv)
	
End Function

' ------------------------------------------------------------------------------------------------
' Multiply Quaternion
' ------------------------------------------------------------------------------------------------
Function MultiplyQuats:TQuaternion(q1:TQuaternion, q2:TQuaternion)

	Local q:TQuaternion = New TQuaternion
	
	q.w = q1.w * q2.w - q1.x * q2.x - q1.y * q2.y - q1.z * q2.z
	q.x = q1.w * q2.x + q1.x * q2.w + q1.y * q2.z - q1.z * q2.y
	q.y = q1.w * q2.y + q1.y * q2.w + q1.z * q2.x - q1.x * q2.z
	q.z = q1.w * q2.z + q1.z * q2.w + q1.x * q2.y - q1.y * q2.x

	Return q

End Function

' ------------------------------------------------------------------------------------------------
' Normalize Quaternion
' ------------------------------------------------------------------------------------------------
Function NormalizeQuat:TQuaternion(q:TQuaternion)

	Local uv:Float = Sqr(q.w * q.w + q.x * q.x + q.y * q.y + q.z * q.z)

	q.w = q.w / uv
	q.x = q.x / uv
	q.y = q.y / uv
	q.z = q.z / uv

	Return q

End Function

' ------------------------------------------------------------------------------------------------
' Events
' ------------------------------------------------------------------------------------------------
Function EventHook:Object(id:Int, Data:Object, Context:Object)

	Local Event:TEvent = TEvent(data)
	If Event = Null Then Return Event

	Select Event.id
							
		Case EVENT_MOUSEWHEEL
						
			If Event.Data > 0 Then

				MouseWheel:+1
				If MouseWheel > 5 Then MouseWheel = 5

				Else

					MouseWheel:-1
				If MouseWheel < 1 Then MouseWheel = 1

			EndIf

	End Select
	
End Function

Comments

Krischan

2017

Patched the code above as I noticed some barely noticeable imprecisions after long runs - the planets slightly ran out of the orbit sometimes. The patch was simple: after the Quaternion Pivot rotation using the "Turn" command just read the pivot angles again, rerotate the pivot using these angles and it was gone. At least I didn't notice any visible imprecisions anymore though there could be some.

' turn pivot = move planet around star
Turn p.piv, 0, (planets + 1 - p.radius) * speed, 0

' patch to fix the "running out of sync" bug due Quaternion imprecision
RotateEntity p.piv, EntityPitch(p.piv), EntityYaw(p.piv), EntityRoll(p.piv)

Blitzplotter

2017

This looks really quite something else, trying it out laters.

Code Archives Forum