Rotating Earth build 2019-04-14

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Post Reply
UEZ
Posts: 988
Joined: May 05, 2017 19:59
Location: Germany

Rotating Earth build 2019-04-14

Post by UEZ »

I've already coded this using GDI+ but now my try to code it without GDI+ using native FB code. That means it should run also on Linux (not tested by me).

Image


Source code:

Code: Select all

'Main code by UEZ
'Image to sphere projection by andalmeida -> https://www.codeproject.com/Articles/19712/Mapping-Images-On-Spherical-Surfaces-Using-C
'Image resize function by Saleth Prakash -> https://www.codeproject.com/Articles/33838/Image-Processing-Using-C
'ASM sinus / cosinus functions code by eukalyptus
'CreateGradientSphere function by Martijn van Iersel -> http://www.helixsoft.nl/articles/sphere/sphere.html
'Thanks to dodicat for tweaking the code and the Regulate function!

#include "vbcompat.bi"

Declare Function _ASM_Cos6th(fX As Double) As Double
Declare Function _ASM_Sin6th(fX As Double) As Double

Type Ang 'FLOATS For angles
    As Double sx,sy,sz
    As Double cx,cy,cz
End Type

Function Construct(x As Double,y As Double,z As Double) As Ang
    Return Type (_ASM_Sin6th(x),_ASM_Sin6th(y),_ASM_Sin6th(z), _
                 _ASM_Cos6th(x),_ASM_Cos6th(y),_ASM_Cos6th(z))
End Function

Declare Sub CreateGlobe(radius As Single, rx As Single, ry As Single, dx As Ushort, dy As Ushort, w As Ushort, h As Ushort, destw As Ushort, desth As Ushort, fEclipticAngle As Single = 0)
Declare Sub RotX(angle As ang, Byref y As Double, Byref z As Double)
Declare Sub RotY(angle As ang, Byref x As Double, Byref z As Double)
Declare Sub RotZ(angle As ang, Byref x As Double, Byref y As Double)
Declare Sub AA2(x As Ushort, y As Ushort, w As Ushort, h As Ushort, iScale As Ubyte = 2)
Declare Sub CreateSphere(imageS As Any Ptr, imageD As Any Ptr, w As Ushort, h As Ushort)
Declare Function ResizeImage(pImage As Any Pointer, iW_new As Ushort, iH_new As Ushort) As Any Pointer
Declare Function RandomRange(fStart As Single, fEnd As Single) As Single
Declare Function CreateGradientSphere(cx As Single, cy As Single, radius As Ushort, lx As Single, ly As Single, zl As Single = 255) As Any Ptr
Declare Function Ecliptic(iDayOfYear As Ushort) As Double
Declare Function Regulate(Byval MyFps As Long, Byref fps As Ushort) As Long

#Define PixelGet(_x, _y)                   	*Cptr(Ulong Ptr, imgDataS + (_y) * pitchS + (_x) Shl 2)
#Define PixelGet2(_x, _y)                  	*Cptr(Ulong Ptr, imgDataD + (_y) * pitch + (_x) Shl 2)
#Define PixelGetM(_x, _y)                  	*Cptr(Ulong Ptr, imgDataM + (_y) * pitchM + (_x) Shl 2)
#Define PixelSet(_x, _y, colour)           	*Cptr(Ulong Ptr, imgDataD + (_y) * pitch + (_x) Shl 2) = (colour)
#Define PixelSetM(_x, _y, colour)          	*Cptr(Ulong Ptr, imgDataM + (_y) * pitchM + (_x) Shl 2) = (colour)
#Define Min(a, b)                      		(Iif(a < b, a, b))
#Define Max(a, b)                      		(Iif(a > b, a, b))
#Define MapCoordinate(i1, i2, w1, w2, p)    (((p - i1) / (i2 - i1)) * (w2 - w1) + w1)
#Define Red(colors)                       	((colors Shr 16) And 255)
#Define Green(colors)                     	((colors Shr 8) And 255)
#Define Blue(colors)                      	(colors And 255)
#Define Alpha(colors)                     	((colors Shr 24) And 255)
#Define Round(x)                      		((x + 0.5) Shr 0)
#Define Floor(x)                      		(((x) * 2.0 - 0.5) Shr 1)


Const phi0 = 0.0, phi1 = Acos(-1), theta0 = 0.0, theta1 = 2.0 * Acos(-1), rad = Acos(-1) / 180, deg = 180 / Acos(-1)
Const iW = 1000, iH = 500, iW2 = iW \ 2, iH2 = iH \ 2, srotx = 1.5 * deg, r = 125, dm = 2 * r, md = (iW - 200)  / 2
Const iWe = 600, iHe = 300, iW2e = iWe \ 2, iH2e = iHe \ 2
Dim As Double fEclipticAngle = Ecliptic(Int(Now() - DateSerial(Year(Now()),1,1)))

'needed for Sub CreateGlobe
Const iMin = Min(iWe, iHe), iMin2 = iMin * 0.33, deltac = 1 + (r \ iMin) Shl 1

Screenres iW, iH, 32, 2
Screenset 1, 0

Windowtitle "FB Rotating Earth coded by UEZ build 2019-04-14"

Dim Shared As Any Ptr pImageS, pImageD, pImageB, pImageG, pImageGM, pImageM, pImageMs, imgDataS, imgDataD, imgDataM
pImageS = Imagecreate(iWe, iHe, 0, 32)
pImageD = Imagecreate(iW, iH, 0, 32)
pImageB = Imagecreate(iW, iH, 0, 32)
pImageM = Imagecreate(200, 200, 0, 32)
pImageGM = Imagecreate(200, 200, 0, 32)


Bload("Earth4_600x300.bmp", pImageS)
Bload("Galaxy_1000x500.bmp", pImageB)
Bload("Moon_200x200.bmp", pImageM)
Dim Shared As Integer pitch, pitchS, pitchG, pitchM

Imageinfo(pImageS, , , , pitchS, imgDataS)  'earth map
Imageinfo(pImageD, , , , pitch, imgDataD)  'final image which will be drawn To Screen
Imageinfo(pImageM, , , , pitchM, imgDataM) 'moon image

'make moon transparent
For xx As Ushort = 0 To 199
	For yy As Ushort = 0 To 199
		If PixelGetM(xx, yy) = &hFFFF0000 Then PixelSetM(xx, yy, 0)
	Next
Next

'blur moon To Get rid of sharp edges
Dim As Ulong resultAlpha, resultRed, resultGreen, resultBlue, col, iScale = 2, gridSize = iScale * iScale
For iY As Ushort = 0 To 199
	For iX As Ushort = 0 To 199
		resultRed = 0: resultGreen = 0: resultBlue = 0: resultAlpha = 0
		For xx As Ubyte = 0 To iScale - 1
			For yy As Ubyte = 0 To iScale - 1
				col = PixelGetM(iX + xx, iY + yy)
				resultRed += Red(col)
				resultGreen += Green(col)
				resultBlue += Blue(col)
				resultAlpha += Alpha(col)
			Next
		Next
		PixelSetM(iX, iY, Rgba(resultRed / gridSize, resultGreen / gridSize, resultBlue / gridSize, resultAlpha / gridSize))
	Next
Next

pImageG = CreateGradientSphere(r, r, r, 80 * rad, 0)


Dim As Single fDeg = -2 * deg, z = 90, z2 = 90, zz, zzz, fz, fs, fx, ca1, fy = iH2 - r / 2,  fSpeedMoon = 0.2

Type tComet
	As Single x, y, vx, vy, len
End Type

Dim As tComet Comet
Dim As Short iLife = 255
Dim As Byte iDir = 1
Randomize , 2

Comet.x = RandomRange(-1, 1)
If Comet.x < 0 Then 
	Comet.x = 0
Else
	Comet.x = iW
	iDir *= -1
End If

Comet.y = 10 + Rnd() * (iH - 20)
Comet.vx = RandomRange(1, 10) * iDir
Comet.vy = RandomRange(1, 10) * iDir
Comet.len = RandomRange(5, 10) * iDir


Dim As Ushort iFPS = 0, iFPS_current
Dim As Single pdx, pdy, dt = 5 + Rnd() * 5

Dim As Double fTimer1, fTimer2
fTimer1 = Timer
fTimer2 = Timer

Do
	Put pImageD, (0, 0), pImageB, Pset 'copy background (galaxy)
	If Timer - fTimer1 > dt Then 'draw a simple meteor flying thru the screen
		pdx = Comet.x + Comet.vx * Comet.len
		pdy = Comet.y + Comet.vy * Comet.len
		Line pImageD, (Comet.x, Comet.y)-(pdx, pdy), Rgb(iLife, iLife, iLife)
		iLife -= Abs(Comet.len \ 2)
		If iLife < 1 Or (Comet.x < 0 And pdx < 0) Or (Comet.x > iW And pdx > iW) Or (Comet.y < 0 And pdy < 0) Or (Comet.y > iH And pdy > iH) Then
			iDir = 1
			Comet.x = RandomRange(-1, 1)
			If Comet.x < 0 Then 
				Comet.x = 0
			Else
				Comet.x = iW
				iDir *= -1
			End If
			Comet.y = 20 + Rnd() * (iH - 40)
			Comet.vx = 1 + Rnd() * 10 * iDir
			Comet.vy = 1 + Rnd() * 10 * iDir
			Comet.len = RandomRange(5, 10) * iDir
			iLife = 255
			dt = 5 + Rnd() * 5
			fTimer1 = Timer
		Else
			Comet.x += Comet.vx
			Comet.y += Comet.vy
		End If
	End If
	
	zz = z * rad
	fx = iW2 - _ASM_Cos6th(zz) * md
	z -= fSpeedMoon
	fz = _ASM_Sin6th(zz)
	fy += _ASM_Sin6th(-zz) * fSpeedMoon * 2.25  
    fS = 100 - fz * 50
	pImageMs = ResizeImage(pImageM, fS, fS) 'resize moon
	ca1 = fs / 2
	PImageGM = CreateGradientSphere(ca1, ca1, (fS + 2) / 2, iW, 0, 255 + fz)
	
	If fz > 0 Then 'moon is in on back side of earth
		Put pImageD, (fx - ca1, fy), pImageMs, Alpha
		Put pImageD, (fx - ca1, fy), pImageGM, Alpha
		CreateGlobe(r - 3, srotx, fDeg, iW2, iH2, iWe, iHe, iW, iH, fEclipticAngle)   
		Put pImageD, (iW2 - r, iH2 - r), pImageG, Alpha
	Else
		CreateGlobe(r - 3, srotx, fDeg, iW2, iH2, iWe, iHe, iW, iH, fEclipticAngle)
		Put pImageD, (iW2 - r, iH2 - r), pImageG, Alpha
		Put pImageD, (fx - ca1, fy), pImageMs, Alpha
		Put pImageD, (fx - ca1, fy), pImageGM, Alpha
	End If

	fDeg -= 0.75

	Put (0, 0), pImageD, Pset 
	Imagedestroy(pImageMs)
	Imagedestroy(pImageGM)

	Draw String(8, 8), iFPS_current & " fps", Rgb(&hF0, &hF0, &hF0)

	If Timer - fTimer2 > 0.99 Then
	  iFPS_current = iFPS
	  iFPS = 0
	  fTimer2 = Timer
	End If

	Flip
	Sleep Regulate(60, iFPS), 1
Loop Until Inkey = Chr(27) 'press ESC To exit

Imagedestroy(pImageS)
Imagedestroy(pImageD)
Imagedestroy(pImageB)
Imagedestroy(pImageG)
Imagedestroy(pImageM)

Function Regulate(Byval MyFps As Long, Byref fps As Ushort) As Long    'code by dodicat
	Static As Double timervalue, _lastsleeptime, t3, frames
	Var t = Timer
	frames += 1
	If (t - t3) >= 1 Then t3 = t : fps = frames : frames = 0
	Var sleeptime = _lastsleeptime + ((1 / myfps) - T + timervalue) * 1000
	If sleeptime < 1 Then sleeptime = 1
	_lastsleeptime = sleeptime
	timervalue = T
	Return sleeptime
End Function

'Obliquity of the ecliptic, only an approximation -> see https://en.wikipedia.org/wiki/Axial_tilt#Earth
Function Ecliptic(iDayOfYear As Ushort) As Double 'iDayOfYear -> 1 to 365
	Return -Cos((iDayOfYear - 1) * 2 * ACos(-1) / 365.25 + 0.205) * 23.44
End Function

Sub CreateGlobe(radius As Single, rx As Single, ry As Single, dx As Ushort, dy As Ushort, w As Ushort, h As Ushort, destw As Ushort, desth As Ushort, fEclipticAngle As Single = 0)
	Dim As Double theta, phi, px, py, pz, c1
	Dim As Short posx, posy
	Dim As Ulong col
	Dim As Single col2
	Dim As Ang a = Construct(rx * rad, ry * rad, fEclipticAngle * rad)
    
	For y As Ushort = 0 To h - 1
		For x As Ushort = 0 To w - 1
			theta = MapCoordinate(0, w - 1, theta1, theta0, x)
			phi = MapCoordinate(0, h - 1, phi0, phi1, y)
			c1 = radius * _ASM_Sin6th(phi)
			px = c1 * _ASM_Cos6th(theta)
			py = c1 * _ASM_Sin6th(theta)
			pz = radius * _ASM_Cos6th(phi)
			RotX(a, py, pz) 'rx
			RotY(a, px, pz) 'ry
			RotZ(a, px, py)
			If pz > 0 Then
				posx = (px + dx)
				posy = (py + dy)
				If posx >= 0 And posx < destw And posy >= 0 And posy < desth Then 
					col = PixelGet(x, y) 
					If radius < iMin2 Then
						PixelSet(posx, posy, col)
					Else
						'Line pImageD, (posx, posy) - (posx + c, posy + c), col, BF 'too slow
						For yy As Ushort = posy To posy + deltac
							For xx As Ushort = posx To posx + deltac
								If xx < destw And yy < desth Then PixelSet(xx, yy, col)
							Next
						Next
						'PixelSet(posx, posy, col)
					End If
				End If
			End If
		Next
	Next   
End Sub

Sub RotX(angle As ang, Byref y As Double, Byref z As Double)
    Dim As Double y1 = y * angle.cx - z * angle.sx, _
				  z1 = y * angle.sx + z * angle.cx
	y = y1
	z = z1
End Sub

Sub RotY(angle As ang, Byref x As Double, Byref z As Double)  
	Dim As Double x1 = x * angle.cy - z * angle.sy, _
				  z1 = x * angle.sy + z * angle.cy             
	x = x1
	z = z1
End Sub

Sub RotZ(angle As ang, Byref x As Double, Byref y As Double) 'not used yet
	Dim As Double x1 = x * angle.cz - y * angle.sz, _
				  y1 = x * angle.sz + y * angle.cz
	x = x1
	y = y1
End Sub

Function RandomRange(fStart As Single, fEnd As Single) As Single
   Return Rnd() * (fEnd - fStart) + fStart
End Function

Function CreateGradientSphere(cx As Single, cy As Single, radius As Ushort, lx As Single, ly As Single, zl As Single = 255) As Any Ptr
	Dim As Any Ptr pImageG, imgDataG
	Dim As Integer pitchG
	pImageG = Imagecreate(radius Shl 1, radius Shl 1, 0, 32)
	Imageinfo(pImageG, , , , pitchG, imgDataG)
	#Define PixelSetG(_x, _y, colour)           	*Cptr(Ulong Ptr, imgDataG + (_y) * pitchG + (_x) Shl 2) = (colour)
	Dim As Single x, y, z, lightx, lighty, lightz, q_cos, light, c, c1 = _ASM_Cos6th(ly)
	lightx = _ASM_Sin6th(lx) * c1
	lighty = _ASM_Sin6th(ly)
	lightz = _ASM_Cos6th(lx) * c1
	For y = -radius To radius
		q_cos = _ASM_Cos6th(-Asin(y / radius)) * radius
		For x = -q_cos To q_cos
			z = Sqr(radius * radius - x * x - y * y)
			c = (x / radius * lightx + y / radius * lighty + z / radius * lightz)
			light = Iif(c < 0, 0, c) * 255
			PixelSetG(Cshort(x + cx), Cshort(y + cy), Rgba(0, 0, 0, Max(Min(255, zl - light), 0)))
		Next
	Next
	Return pImageG
End Function

Function _ASM_Sin6th(fX As Double) As Double 
	'By Eukalyptus 
	Asm
		jmp 0f
		1: .Double 683565275.57643158 
		2: .Double -0.0000000061763971109087229 
		3: .Double 6755399441055744.0 
		  
		0: 
			movq xmm0, [fX] 
			mulsd xmm0, [1b] 
			addsd xmm0, [3b] 
			movd ebx, xmm0 

			lea  eax, [ebx*2+0x80000000] 
			sar  eax, 2 
			imul eax 
			sar  ebx, 31 
			lea  eax, [edx*2-0x70000000] 
			lea  ecx, [edx*8+edx-0x24000000] 
			imul edx 
			Xor  ecx, ebx 
			lea  eax, [edx*8+edx+0x44A00000]
			imul ecx 

			cvtsi2sd xmm0, edx 
			mulsd xmm0, [2b] 
			movq [Function], xmm0 
	End Asm 
End Function

Function _ASM_Cos6th(fX As Double) As Double 
	'By Eukalyptus 
	Asm 
		jmp 0f 
		1: .Double 683565275.57643158 
		2: .Double -0.0000000061763971109087229 
		3: .Double 6755399441055744.0 

		0: 
			movq xmm0, [fX] 
			mulsd xmm0, [1b] 
			addsd xmm0, [3b] 
			movd ebx, xmm0 

			Add ebx, 0x40000000 'SinToCos 

			lea  eax, [ebx*2+0x80000000] 
			sar  eax, 2 
			imul eax 
			sar  ebx, 31 
			lea  eax, [edx*2-0x70000000] 
			lea  ecx, [edx*8+edx-0x24000000] 
			imul edx 
			Xor  ecx, ebx 
			lea  eax, [edx*8+edx+0x44A00000] 
			imul ecx 

			cvtsi2sd xmm0, edx 
			mulsd xmm0, [2b] 
			movq [Function], xmm0 
	End Asm 
End Function

'https://www.codeproject.com/Articles/33838/Image-Processing-Using-C
Private Function ResizeImage(pImage As Any Pointer, iW_new As Ushort, iH_new As Ushort) As Any Pointer
	#Define GetPixelRI(_x, _y)          *Cptr(Ulong Ptr, imgData + _y * pitch + _x Shl 2)
	#Define SetPixelRI(_x, _y, _color)  *Cptr(Ulong Ptr, imgData_Resized + _y * pitch_Resized + _x Shl 2) = (_color)

	Dim pImage_Resized As Any Ptr = Imagecreate(iW_new, iH_new, 0, 32)
	Dim As Integer w, h, wr, hr, pitch, pitch_Resized
	Dim As Any Pointer imgData, imgData_Resized

	Imageinfo(pImage, w, h, , pitch, imgData)
	Imageinfo(pImage_Resized, wr, hr, , pitch_Resized, imgData_Resized)

	Dim As Single fWidthFactor = w / wr, fHeightFactor = h / hr, fx, fy, nx, ny
	Dim As Ulong cx, cy, fr_x, fr_y, color1, color2, color3, color4
	Dim As Ubyte nRed, nGreen, nBlue, nAlpha, bp1, bp2

	For x As Ushort = 0 To wr - 1
		For y As Ushort = 0 To hr - 1

			fr_x = Floor(x * fWidthFactor)
			fr_y = Floor(y * fHeightFactor)

			cx = fr_x + 1
			If cx >= w Then cx = fr_x
			cy = fr_y + 1
			If cy >= h Then cy = fr_y

			fx = x * fWidthFactor - fr_x
			fy = y * fHeightFactor - fr_y

			nx = 1.0 - fx
			ny = 1.0 - fy

			color1 = GetPixelRI(fr_x, fr_y)
			color2 = GetPixelRI(cx, fr_y)
			color3 = GetPixelRI(fr_x, cy)
			color4 = GetPixelRI(cx, cy)

			'red
			bp1 = nx * Red(color1) + fx * Red(color2)
			bp2 = nx * Red(color3) + fx * Red(color4)         
			nRed = ny * bp1 + fy * bp2         

			'green
			bp1 = nx * Green(color1) + fx * Green(color2)
			bp2 = nx * Green(color3) + fx * Green(color4)         
			nGreen = ny * bp1 + fy * bp2

			'blue
			bp1 = nx * Blue(color1) + fx * Blue(color2)
			bp2 = nx * Blue(color3) + fx * Blue(color4)
			nBlue = ny * bp1 + fy * bp2

			'Alpha
			bp1 = nx * Alpha(color1) + fx * Alpha(color2)
			bp2 = nx * Alpha(color3) + fx * Alpha(color4)         
			nAlpha = ny * bp1 + fy * bp2

			SetPixelRI(x, y, Rgba(nRed, nGreen, nBlue, nAlpha))
		Next
	Next
	Return pImage_Resized
End Function
To run the code properly you need the images which can be downloaded here: FB Rotating Earth build 2019-04-11.zip.

On my old notebook it runs best with -gen gcc -Wc -Ofast (x64) at ~68 fps.


Btw, it is not a realistic simulation. ^^
Last edited by UEZ on Apr 15, 2019 7:21, edited 9 times in total.
dodicat
Posts: 7983
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Rotating Earth build 2019-04-08

Post by dodicat »

Thanks UEZ.
Very nice.
If you don't mind, I have tweaked a bit to give shorter days (and nights)
I have left rotz alone, you do not use it.
Also I note there is no Mercator moon bitmap to spread out and see the far side, but I nit pick.

Code: Select all

'Main code by UEZ
'Image to sphere projection by andalmeida -> https://www.codeproject.com/Articles/19712/Mapping-Images-On-Spherical-Surfaces-Using-C
'Image resize function by Saleth Prakash -> https://www.codeproject.com/Articles/33838/Image-Processing-Using-C

Type ang 'FLOATS for angles
    As Single sx,sy,sz
    As Single cx,cy,cz
End Type

Function construct(x As Single,y As Single,z As Single) As Ang
    Return   Type (Sin(x),Sin(y),Sin(z), _
                   Cos(x),Cos(y),Cos(z))
End Function

Declare Sub CreateGlobe(imageS As Any Ptr, imageD As Any Ptr, radius As Single, rx As Single, ry As Single, dx As Short, dy As Short, w As Ushort, h As Ushort)
Declare Sub RotX(angle As ang, Byref y As Double, Byref z As Double)
Declare Sub RotY(angle As ang, Byref x As Double, Byref z As Double)
Declare Sub RotZ(angle As Single, Byref x As Double, Byref y As Double)
Declare Sub AA2(x As Ushort, y As Ushort, w As Ushort, h As Ushort, iScale As Ubyte = 2)
Declare Sub CreateSphere(imageS As Any Ptr, imageD As Any Ptr, w As Ushort, h As Ushort)
Declare Function ResizeImage(pImage As Any Pointer, iW_new As Ushort, iH_new As Ushort) As Any Pointer
Declare Function _ASM_Cos6th(fX As Double) As Double
Declare Function _ASM_Sin6th(fX As Double) As Double


#Define PixelGet(_x, _y)            		*Cptr(Ulong Ptr, imgDataS + (_y) * pitch + (_x) Shl 2)
#Define PixelGet2(_x, _y)            		*Cptr(Ulong Ptr, imgDataD + (_y) * pitch + (_x) Shl 2)
#Define PixelGetM(_x, _y)            		*Cptr(Ulong Ptr, imgDataM + (_y) * pitchM + (_x) Shl 2)
#Define PixelSet(_x, _y, colour)    		*Cptr(Ulong Ptr, imgDataD + (_y) * pitch + (_x) Shl 2) = (colour)
#Define PixelSetG(_x, _y, colour)    		*Cptr(Ulong Ptr, imgDataG + (_y) * pitchG + (_x) Shl 2) = (colour)
#Define PixelSetM(_x, _y, colour)    		*Cptr(Ulong Ptr, imgDataM + (_y) * pitchM + (_x) Shl 2) = (colour)
#Define Min(a, b) 							(Iif(a < b, a, b))
#Define Max(a, b) 							(Iif(a > b, a, b))
#Define MapCoordinate(i1, i2, w1, w2, p) 	(((p - i1) / (i2 - i1)) * (w2 - w1) + w1)
#Define Red(colors)                 		((colors Shr 16) And 255)
#Define Green(colors)               		((colors Shr 8) And 255)
#Define Blue(colors)                		(colors And 255)
#Define Alpha(colors)                		((colors Shr 24) And 255)
#Define Round(x) 							((x + 0.5) Shr 0)
#Define Floor(x) 							(((x) * 2.0 - 0.5) Shr 1)

Const phi0 = 0.0, phi1 = Acos(-1), theta0 = 0.0, theta1 = 2.0 * Acos(-1), rad = Acos(-1) / 180, deg = 180 / Acos(-1)

Const iW = 600, iH = 300, iW2 = iW \ 2, iH2 = iH \ 2, srotx = 1.5 * deg, r = 135, dm = 2 * r
Screenres iW, iH, 32, 2
Screenset 1, 0

Windowtitle "FB Rotating Earth coded by UEZ build 2019-04-08"

Dim Shared As Any Ptr pImageS, pImageD, pImageB, pImageG, pImageM, pImageMs, imgDataS, imgDataD, imgDataG, imgDataM
pImageS = Imagecreate(iW, iH, 0, 32)
pImageD = Imagecreate(iW, iH, 0, 32)
pImageB = Imagecreate(iW, iH, 0, 32)
pImageM = Imagecreate(200, 200, 0, 32)
pImageG = Imagecreate(dm + 2, dm + 2, 0, 32)


Bload("Earth4_600x300.bmp", pImageS)
Bload("Galaxy_MilkyWay.bmp", pImageB)
Bload("Moon_200x200.bmp", pImageM)
Dim Shared As Integer pitch, pitchG, pitchM

Imageinfo(pImageS, , , , pitch, imgDataS) 'earth map
Imageinfo(pImageD, , , , pitch, imgDataD) 'final image which will be drawn to screen
Imageinfo(pImageG, , , , pitchG, imgDataG) 'final image of transformation of earth map
Imageinfo(pImageM, , , , pitchM, imgDataM) 'moon image

'make moon transparent
For xx As Ushort = 0 To 199
	For yy As Ushort = 0 To 199
		If PixelGetM(xx, yy) = &hFFFF0000 Then PixelSetM(xx, yy, 0)
	Next
Next

'blur moon to get rid of sharp edges
Dim As Ulong resultAlpha, resultRed, resultGreen, resultBlue, col, iScale = 2, gridSize = iScale * iScale
For iY As Ushort = 0 To 199
	For iX As Ushort = 0 To 199
		resultRed = 0: resultGreen = 0: resultBlue = 0: resultAlpha = 0
		For xx As Ubyte = 0 To iScale - 1
			For yy As Ubyte = 0 To iScale - 1
				col = PixelGetM(iX + xx, iY + yy)
				resultRed += Red(col)
				resultGreen += Green(col)
				resultBlue += Blue(col)
				resultAlpha += Alpha(col)
			Next
		Next
		PixelSetM(iX, iY, Rgba(resultRed / gridSize, resultGreen / gridSize, resultBlue / gridSize, resultAlpha / gridSize))
	Next
Next

'create gradient ellipse
Dim As Single r2 = (r + 3) * (r + 3), X0 = r, Y0 = r
Dim As Long X, Y, D2
Dim As Ulong iCol
For y = -r To r
	For x = -r To r
		D2 = x * x + y * y
		If D2 <= r2 Then                              
			iCol = Max(0, 320 - Sqr(r2 - D2) * 3.5)
			PixelSetG(Cushort(x + X0), Cushort(y + Y0), iCol Shl 24)
		End If
	Next 
Next 


Dim As Ushort iFPS = 0, iFPS_current
Dim As Single fDeg = -2 * deg, z = 90, zz, fz, fs, fx, fy = 30
Dim As Double fTimer2
fTimer2 = Timer


Do
	Put pImageD, (0, 0), pImageB, Pset 'copy background (galaxy)
	
	zz = z * rad
	fx = iW2 - _ASM_Cos6th(zz) * iW2
	z += 0.5
	fz = _ASM_Sin6th(zz)
	fS = 100 - fz * 50
	pImageMs = ResizeImage(pImageM, fs, fs) 'resize moon
	
	If fz > 0 Then 'moon is in the back of earth
		Put pImageD, (fx - fs \ 2, fy), pImageMs, Alpha
		CreateGlobe(pImageS, pImageD, r - 1, srotx, fDeg, iW2, iH2, iW, iH)	
		Put pImageD, (iW2 - r, iH2 - r), pImageG, Alpha
	Else
		CreateGlobe(pImageS, pImageD, r - 1, srotx, fDeg, iW2, iH2, iW, iH)
		Put pImageD, (iW2 - r, iH2 - r), pImageG, Alpha
		Put pImageD, (fx - fs \ 2, fy), pImageMs, Alpha
	End If

	fDeg -= 0.1
	
	Put (0, 0), pImageD, Pset 
	Imagedestroy(pImageMs)
	
	Draw String(8, 8), iFPS_current & " fps", Rgb(&hF0, &hF0, &hF0)

	If Timer - fTimer2 > 0.99 Then
		iFPS_current = iFPS
		iFPS = 0
		fTimer2 = Timer
	End If
	iFPS += 1
   
	Flip
	Sleep(1, 1)
Loop Until Inkey = Chr(27) 'press ESC to exit

Imagedestroy(pImageS)
Imagedestroy(pImageD)
Imagedestroy(pImageB)
Imagedestroy(pImageG)
Imagedestroy(pImageM)


Sub CreateGlobe(imageS As Any Ptr, imageD As Any Ptr, radius As Single, rx As Single, ry As Single, dx As Short, dy As Short, w As Ushort, h As Ushort)
	Dim As Double theta, phi, px, py, pz, c1
	Dim As Short posx, posy, iMin = Min(w, h), iMin2 = iMin * 0.33, c = 1 + (radius \ iMin) Shl 1
	Dim As Ulong col
	Dim As Single col2
    dim as ang a=construct(rx*rad,ry*rad,0)
    
	For y As Ushort = 0 To h - 1
		For x As Ushort = 0 To w - 1
			theta = MapCoordinate(0, w - 1, theta1, theta0, x)
			phi = MapCoordinate(0, h - 1, phi0, phi1, y)
			c1 = radius * _ASM_Sin6th(phi)
			px = c1 * _ASM_Cos6th(theta)
			py = c1 * _ASM_Sin6th(theta)
			pz = radius * _ASM_Cos6th(phi)
			RotX(a, py, pz)'rx
			RotY(a, px, pz)'ry
			If pz > 0 Then
				posx = (px + dx)
				posy = (py + dy)
				If posx >= 0 And posx < w And posy >= 0 And posy < h Then 
					col = PixelGet(x, y) 
					If radius < iMin2 Then
						PixelSet(posx, posy, col)
					Else
						'Line pImageD, (posx, posy) - (posx + c, posy + c), col, BF 'too slow
						For yy As Ushort = posy To posy + c
							For xx As Ushort = posx To posx + c
								If xx < w And yy < h Then PixelSet(xx, yy, col)
							Next
						Next
						'PixelSet(posx, posy, col)
					End If
				End If
			End If
		Next
	Next	
End Sub

Sub RotX(angle As ang, Byref y As Double, Byref z As Double)
	'Dim As Single y1 = y * Cos(angle) - z * Sin(angle), _
				  'z1 = y * Sin(angle) + z * Cos(angle)
                  Dim As Single y1 = y *angle.cx - z * angle.sx, _
				  z1 = y * angle.sx + z * angle.cx
	y = y1
	z = z1
End Sub

Sub RotY(angle As ang, Byref x As Double, Byref z As Double)
	'Dim As Single x1 = x * Cos(angle) - z * Sin(angle), _
				  'z1 = x * Sin(angle) + z * Cos(angle)
                  
     Dim As Single x1 = x * angle.cy - z * angle.sy, _
				  z1 = x * angle.sy + z * angle.cy             
	x = x1
	z = z1
End Sub

Sub RotZ(angle As Single, Byref x As Double, Byref y As Double)
	Dim As Single x1 = x * Cos(angle) - y * Sin(angle), _
				  y1 = x * Sin(angle) + y * Cos(angle)
	x = x1
	y = y1
End Sub

Function _ASM_Sin6th(fX As Double) As Double 
	'By Eukalyptus 
   Asm
      jmp 0f
      1: .Double 683565275.57643158 
      2: .Double -0.0000000061763971109087229 
      3: .Double 6755399441055744.0 
          
      0: 
         movq xmm0, [fX] 
         mulsd xmm0, [1b] 
         addsd xmm0, [3b] 
         movd ebx, xmm0 

         lea  eax, [ebx*2+0x80000000] 
         sar  eax, 2 
         imul eax 
         sar  ebx, 31 
         lea  eax, [edx*2-0x70000000] 
         lea  ecx, [edx*8+edx-0x24000000] 
         imul edx 
         Xor  ecx, ebx 
         lea  eax, [edx*8+edx+0x44A00000]
         imul ecx 

         cvtsi2sd xmm0, edx 
         mulsd xmm0, [2b] 
         movq [Function], xmm0 
   End Asm 
End Function

Function _ASM_Cos6th(fX As Double) As Double 
   'By Eukalyptus 
   Asm 
      jmp 0f 
         1: .Double 683565275.57643158 
         2: .Double -0.0000000061763971109087229 
         3: .Double 6755399441055744.0 
       
      0: 
         movq xmm0, [fX] 
         mulsd xmm0, [1b] 
         addsd xmm0, [3b] 
         movd ebx, xmm0 
          
         Add ebx, 0x40000000 'SinToCos 
    
         lea  eax, [ebx*2+0x80000000] 
         sar  eax, 2 
         imul eax 
         sar  ebx, 31 
         lea  eax, [edx*2-0x70000000] 
         lea  ecx, [edx*8+edx-0x24000000] 
         imul edx 
         Xor  ecx, ebx 
         lea  eax, [edx*8+edx+0x44A00000] 
         imul ecx 
          
         cvtsi2sd xmm0, edx 
         mulsd xmm0, [2b] 
         movq [Function], xmm0 
   End Asm 
End Function

'https://www.codeproject.com/Articles/33838/Image-Processing-Using-C
Private Function ResizeImage(pImage As Any Pointer, iW_new As Ushort, iH_new As Ushort) As Any Pointer
	#Define GetPixelRI(_x, _y)          *Cptr(Ulong Ptr, imgData + _y * pitch + _x Shl 2)
	#Define SetPixelRI(_x, _y, _color)  *Cptr(Ulong Ptr, imgData_Resized + _y * pitch_Resized + _x Shl 2) = (_color)

	Dim pImage_Resized As Any Ptr = Imagecreate(iW_new, iH_new, 0, 32)
	Dim As Integer w, h, wr, hr, pitch, pitch_Resized
	Dim As Any Pointer imgData, imgData_Resized

	Imageinfo(pImage, w, h, , pitch, imgData)
	Imageinfo(pImage_Resized, wr, hr, , pitch_Resized, imgData_Resized)

	Dim As Single fWidthFactor = w / wr, fHeightFactor = h / hr, fx, fy, nx, ny
	Dim As Ulong cx, cy, fr_x, fr_y, color1, color2, color3, color4
	Dim As Ubyte nRed, nGreen, nBlue, nAlpha, bp1, bp2

	For x As Ushort = 0 To wr - 1
		For y As Ushort = 0 To hr - 1

			fr_x = Floor(x * fWidthFactor)
			fr_y = Floor(y * fHeightFactor)

			cx = fr_x + 1
			If cx >= w Then cx = fr_x
			cy = fr_y + 1
			If cy >= h Then cy = fr_y

			fx = x * fWidthFactor - fr_x
			fy = y * fHeightFactor - fr_y

			nx = 1.0 - fx
			ny = 1.0 - fy

			color1 = GetPixelRI(fr_x, fr_y)
			color2 = GetPixelRI(cx, fr_y)
			color3 = GetPixelRI(fr_x, cy)
			color4 = GetPixelRI(cx, cy)

			'red
			bp1 = nx * Red(color1) + fx * Red(color2)
			bp2 = nx * Red(color3) + fx * Red(color4)         
			nRed = ny * bp1 + fy * bp2         

			'green
			bp1 = nx * Green(color1) + fx * Green(color2)
			bp2 = nx * Green(color3) + fx * Green(color4)         
			nGreen = ny * bp1 + fy * bp2

			'blue
			bp1 = nx * Blue(color1) + fx * Blue(color2)
			bp2 = nx * Blue(color3) + fx * Blue(color4)
			nBlue = ny * bp1 + fy * bp2

			'Alpha
			bp1 = nx * Alpha(color1) + fx * Alpha(color2)
			bp2 = nx * Alpha(color3) + fx * Alpha(color4)         
			nAlpha = ny * bp1 + fy * bp2

			SetPixelRI(x, y, Rgba(nRed, nGreen, nBlue, nAlpha))
		Next
	Next
   Return pImage_Resized
End Function 
UEZ
Posts: 988
Joined: May 05, 2017 19:59
Location: Germany

Re: Rotating Earth build 2019-04-08

Post by UEZ »

Thanks for your feedback dodicat. :-)

No, I don't mind of course and I also wanted to ask how to tweak the code that it runs faster. That's a nice idea to use the construct which give a great performance boost to fps - now at ~ 77 fps (doubled+). I need to analyse and check it to understand the idea of that coding technique.

I will update the code soon.

Now I can add some more elements...
Last edited by UEZ on Apr 08, 2019 20:34, edited 1 time in total.
Tourist Trap
Posts: 2958
Joined: Jun 02, 2015 16:24

Re: Rotating Earth build 2019-04-08

Post by Tourist Trap »

dodicat wrote:Thanks UEZ.
Very nice.
If you don't mind, I have tweaked a bit to give shorter days (and nights)
Hi all,

@dodicat, your snippet doesn't use GDI+ it seems? If so, I'll try to learn from it because I was looking for a planet for my game !
UEZ
Posts: 988
Joined: May 05, 2017 19:59
Location: Germany

Re: Rotating Earth build 2019-04-08

Post by UEZ »

Tourist Trap wrote:
dodicat wrote:Thanks UEZ.
Very nice.
If you don't mind, I have tweaked a bit to give shorter days (and nights)
Hi all,

@dodicat, your snippet doesn't use GDI+ it seems? If so, I'll try to learn from it because I was looking for a planet for my game !
Neither the code I posted nor the tweaked version from dodicat uses GDI+. I only mentioned that my previous version was based on GDI+.
Tourist Trap
Posts: 2958
Joined: Jun 02, 2015 16:24

Re: Rotating Earth build 2019-04-08

Post by Tourist Trap »

UEZ wrote: Neither the code I posted nor the tweaked version from dodicat uses GDI+. I only mentioned that my previous version was based on GDI+.
Hi UEZ,

thanks, my fear was to use a system specific tool for a game I would I would rather keep portable. Otherwise, I like learning gdi+, it's very useful on windows and I think it should be used for professional projects.
UEZ
Posts: 988
Joined: May 05, 2017 19:59
Location: Germany

Re: Rotating Earth build 2019-04-11

Post by UEZ »

Little update on 2019-04-11:
* added simple meteor movement
* sun light direction is now from the right side
* ecliptic calculation (earth axis)

If you are interested have a look to the 1st post. You need to download the zip-archive because I changed the background size.
Post Reply