from https://love2d.org/forums/viewtopic.php?t=77487

planet_noise.lua

local lg = love.graphics
local lm = love.math

local pmethods={}

local function earth(u) 
	local cr,cb,cg
	if u < 0.5 then
		cb = u + 0.5 
		cr=0
		cg=0
	else
		cb = 0
		cg = u
		cr = 2*(u -0.5)
	end
	return cr,cg,cb
end

local function ocean(u)
	return earth(u*u)
end

local function bizare1(u)
	return u*u,u * (1-u)*4,(1-u)*(1-u)
end

local function bizare2(u)
	return bizare1(u*u)
end
				

local colFuncs={earth,bizare1,ocean,bizare2}

local function getPlanetImages(radius,nbFrames)

	local colFunc = colFuncs[math.random(1,#colFuncs)]
	
	local ox= 50*math.random() - 25
	local oy= 50*math.random() - 25
	local oz= 50*math.random() - 25

	local angRot = math.pi*2/nbFrames
	
	local sz =8
	while 2*radius > sz do sz = sz *2 end

	-- determination of the rotation axe
	local ux = 2*math.random() - 1
	local uy = 2*math.random() - 1
	local uz = 2*math.random() - 1
	local ul = math.sqrt(ux*ux+uy*uy+uz*uz)
	ux=ux/ul
	uy=uy/ul
	uz=uz/ul
	
	local m11,m12,m13 = 1,0,0
	local m21,m22,m23 = 0,1,0
	local m31,m32,m33 = 0,0,1
	
	local r11,r12,r13,r21,r22,r23,r31,r32,r33
	local n11,n12,n13,n21,n22,n23,n31,n32,n33

	-- determination of the rotation matrix
	do 
		local c = math.cos(angRot)
		local s = math.sin(angRot)
		local cc = 1-c
		local xs,ys,zs = ux*s,uy*s,uz*s
		local xc,yc,zc = ux*cc,uy*cc,uz*cc
		local xyc,yzc,zxc = ux*yc, uy*zc, uz*xc
		r11,r12,r13 = ux*xc+c, xyc -zs, zxc+ys
		r21,r22,r23 = xyc + zs, uy * yc + c, yzc -xs
		r31,r32,r33 = zxc - ys, yzc+xs, uz*zc + c
	end
	local id = love.image.newImageData(sz*nbFrames,sz)
	local xc = (2*radius-1)/2
	local yc = xc
	local r2=radius*radius
	for iframe=0,nbFrames-1 do
	for x=0,2*radius-1 do
		local u2 = (x-xc)*(x-xc)
		local x0=(x-xc)/radius
		for y=0,2*radius-1 do
			local d2=u2+(y-yc)*(y-yc)
			if d2 <=r2 then
				z=math.sqrt(r2-d2)
				local l=((xc-x)+(yc-y)+z*math.sqrt(2))*0.5/radius
				if l<0 then
					l=0
				end
				l=(0.2 + 0.8 * l) * 255
				
				local y0 = (y-yc)/radius 
				local z0 = z/radius 
				
				local xx = ox + m11 * x0 + m12 * y0 + m13 * z0
				local yy = oy + m21 * x0 + m22 * y0 + m23 * z0
				local zz = oz + m31 * x0 + m32 * y0 + m33 * z0
				
				local u1 = lm.noise(xx,yy,zz)
				local u2 = lm.noise(2*xx,2*yy,2*zz)
				local u3 = lm.noise(4*xx,4*yy,4*zz)
				local u4 = lm.noise(8*xx,8*yy,8*zz)
				local u = u1 + u2/2 + u3/4 + u4/8
				u=u/1.875
				--u=u*u
				local cr,cg,cb = colFunc(u)
				--[[
				cr = u*u
				cb = (1-u)*(1-u)
				cg = u * (1-u)*4
				
				if u < 0.5 then
					cb = u + 0.5 
					cr=0
					cg=0
				else
					cb = 0
					cg = u
					cr = 2*(u -0.5)
				end
				--]]
				id:setPixel(x+sz*iframe,y,cr *l,cg*l,cb*l,255)
			end
		end
	end
	
	    -- matrix multiplication porn !
	    n11=m11*r11+m12*r21+m13*r31
		n12=m11*r12+m12*r22+m13*r32
		n13=m11*r13+m12*r23+m13*r33

		n21=m21*r11+m22*r21+m23*r31
		n22=m21*r12+m22*r22+m23*r32
		n23=m21*r13+m22*r23+m23*r33

		n31=m31*r11+m32*r21+m33*r31
		n32=m31*r12+m32*r22+m33*r32
		n33=m31*r13+m32*r23+m33*r33
		
		m11,m12,m13 = n11,n12,n13
		m21,m22,m23 = n21,n22,n23
		m31,m32,m33 = n31,n32,n33
		
	end
	return lg.newImage(id)
end


function pmethods:draw()
	if self.img then
		lg.draw(self.img,self.quad,self.x-self.radius,self.y-self.radius)
	end
end

function pmethods:update(dt)
	if self.img then
		
		local x = self.x + self.vx * dt
		local y = self.y + self.vy * dt
		local r = self.radius
		
		if x-r < 0 then
			self.vx = -self.vx
			x = r
		end
		if y-r < 0 then
			self.vy = -self.vy
			y=r
		end
		if x+r > SW then
			self.vx=-self.vx
			x=SW-r
		end
		if y+r > SH then
			self.vy = -self.vy
			y=SH-r
		end
		self.x=x
		self.y=y
		self.timer=self.timer+dt
		if self.timer > self.trig then
			self.timer = self.timer - self.trig
			self.iframe = self.iframe + 1
			if self.iframe >= self.nbFrames then
				self.iframe = 0
			end
			self:setQuad()
		end
	end
end


function pmethods:setQuad()
	local sz = self.sz
	self.quad = lg.newQuad(self.iframe*sz,0,sz,sz,sz*self.nbFrames,sz)
end

function pmethods:generate()
	self.nbFrames = 64
	self.radius = math.random(30,45)
	self.img =  getPlanetImages(self.radius,self.nbFrames)
	self.sz = self.img:getHeight()
	self.iframe=0
	self:setQuad()
	self.timer=0
	self.trig = self.delay/self.nbFrames
	
	local alpha = 2*math.pi*math.random()
	local speed = math.random(50,150)
	self.vx = speed * math.cos(alpha)
	self.vy = speed * math.sin(alpha)

end

local function newPlanet(x,y)
	local p=setmetatable({x=x,y=y},{__index=pmethods})
	p.delay=3
	return p
end

return newPlanet