This script simulates the size shift of a periodic variable star (48 Aurigae).
It uses a B-Spline to fit the star's radius vs. time data points into a smoother curve.
Unfortunately, star brightness and temperature can't be changed from within LUA, so this aspects weren't modelled.
The visual simulation has been exaggerated 5 times in order to make the effect more visible. Please change FACTOR to 1 in the code below. This will show you the real scale of this phenomenon.
The displayed values for radius are real, though.
Code: Select all
--***************************************************************************
--***************************************************************************
-- An open-uniform B-Spline approximation algorithm
-- Applied to (limited) Variable Star Simulation
--
-- Coded by Toti
--
--***************************************************************************
--***************************************************************************
--***************************************************************************
-- Constants
--***************************************************************************
STAR = "48 Aur"
PERIOD = 3 -- period of this star's light curve
DISTANCE = 12 -- * mean star's radius
FACTOR = 5 -- for a more realistic view, this must be 1
ITERATIONS = 100
KM_MLY = 9466411.842
KM_AU = 149597870.7
TIMESPEED = 100000
ANGSPEED = 0.25
UP = celestia:newvector (0, 1, 0)
OBS = celestia:getobserver()
ORDER = 4 -- order of the curve. Increase for additional smoothness
VARIAB = { { 0.000, -310000 , 0}, -- Control points, in the form (phase, delta Radius, 0)
{ 0.040, -90000 , 0}, -- data extracted from http://mb-soft.com/public2/cepheid.html
{ 0.080, 116000 , 0},
{ 0.120, 300000 , 0},
{ 0.160, 459000 , 0},
{ 0.200, 589000 , 0},
{ 0.240, 689000 , 0},
{ 0.280, 760000 , 0},
{ 0.320, 801000 , 0},
{ 0.360, 814000 , 0},
{ 0.400, 798000 , 0},
{ 0.440, 753000 , 0},
{ 0.480, 682000 , 0},
{ 0.520, 582000 , 0},
{ 0.560, 455000 , 0},
{ 0.600, 300000 , 0},
{ 0.640, 121000 , 0},
{ 0.680, -80000 , 0},
{ 0.720, -293000 , 0},
{ 0.760, -506000 , 0},
{ 0.800, -689000 , 0},
{ 0.840, -801000 , 0},
{ 0.880, -806000 , 0},
{ 0.920, -702000 , 0},
{ 0.960, -524000 , 0},
{ 1.000, -310000 , 0} }
-- *********************************************************
-- B-Spline functions
-- *********************************************************
function buildKnots (ord, size)
-- builds a table of open-uniform knots
local j
local K = {}
local kn = size + ord
local step = size / (size - ord + 1)
for j = 1, kn do
if j <= ord then
K[j] = 0
elseif j > size then
K[j] = size
else
K[j] = (j - ord) * step
end
end
return K, size
end
function N (K, i, ord, curr)
-- returns the blending factor for 'curr' position along the curve, using 'K' as knot vector. Uses recursion
local blend
if ord == 1 then -- definition 1st case
if K[i] <= curr and curr < K[i+1] then
blend = 1
else
blend = 0
end
else
local k1 = i+ord-1
local k2 = i+ord
if K[k1] == K[i] and K[k2] == K[i+1] then -- be careful with 'division by zero' errors
blend = 0
elseif K[k1] == K[i] then
blend = (K[k2] - curr) / (K[k2] - K[i+1]) * N(K, i+1, ord-1, curr)
elseif K[k2] == K[i+1] then
blend = (curr - K[i]) / (K[k1] - K[i]) * N(K, i, ord-1, curr)
else -- definition 2nd case
blend = (curr - K[i]) / (K[k1] - K[i]) * N(K, i, ord-1, curr) + (K[k2] - curr) / (K[k2] - K[i+1]) * N(K, i+1, ord-1, curr)
end
end
return blend
end
function buildSpline (K, CP, rows, cols, ord, curr)
-- returns the B-Spline n-tuple value for 'curr' position in knot vector K
local i, j
local S = {}
for j = 1, cols do
S[j] = 0
end
for i = 1, rows do
local w = N(K, i, ord, curr)
for j = 1, cols do
S[j] = CP[i][j] * w + S[j]
end
end
return S
end
--***************************************************************************
-- Other functions
--***************************************************************************
function sphe2cart (rad, rho, phi)
-- converts from spherical to cartesian coordinates. Returns a POSITION
local r = rad / KM_MLY
local x = - r * math.cos (math.rad (rho)) * math.cos (math.rad (phi))
local y = r * math.sin (math.rad (phi))
local z = r * math.sin (math.rad (rho)) * math.cos (math.rad (phi))
return celestia:newposition (x,y,z)
end
--***************************************************************************
-- Main routine
--***************************************************************************
tableSize = table.getn (VARIAB)
datumSize = table.getn (VARIAB[1])
Knots = buildKnots (ORDER, tableSize)
star = celestia:find (STAR)
celestia:select (star)
starName = star:name()
starRadius = star:radius()
solRadius = celestia:find ("Sol"):radius()
dist = DISTANCE * starRadius
startDate = celestia:gettime()
celestia:settimescale (TIMESPEED)
celestia:settime (startDate)
t1 = 0
for iteration = 0, ITERATIONS do
repeat
starPos = star:getposition() -- is the star moving?
t1Rel = t1 / PERIOD * tableSize
rad = buildSpline (Knots, VARIAB, tableSize, datumSize, ORDER, t1Rel)[2]
numberRad = rad + starRadius -- radius for numeric display
visualRad = rad * FACTOR + starRadius -- radius for visual display
-- slow rotation around the star:
radPos = dist - visualRad
phi = ANGSPEED * celestia:getscripttime() -- a short way to avoid big numbers unstability
theta = -45
obsPos = sphe2cart (radPos, phi, theta) + starPos
OBS:setposition (obsPos)
OBS:lookat (starPos, UP)
celestia:flash (string.format("%s's radius: \n %5.3f AU\n %5.1f RSun\n %5.1f percent of mean radius ", starName, numberRad/KM_AU, numberRad/solRadius, numberRad/starRadius*100), 5)
wait()
t1 = celestia:gettime() - startDate -- put this here so old computers (like mine) can run the code fine
until t1 >= PERIOD
t1 = t1 - PERIOD -- cut leftovers
startDate = startDate + PERIOD
end
celestia:flash ("That was all!", 5)
celestia:settimescale (1)
EDIT: should be numerically stable now
Please report bugs, sugestions and comments here.
Bye
