Improvements...

[Antoni Gual]

The version you once posted at progboards works ok.
Maybe it's just a matter of cut & paste.. If you lost it in your pc crash i can mail it to you this afternoon.

[Jark]

I have kept a back-up of all the versions of the TC-Lib, and this routine was in version 1.3...

The best I can do is to test the routine as I did before: random generation of the coefficients, and check the result for each root...

These wrong pixels (there are still a few of them) really puzzle me...

I have also initiated TC-Land, trying to understand Hugo Elias logics for Cloud Covers...

and I have found THE trick to raytrace a perlin moutain without that horrible slow algorithm I used for QB Ground... A pixel above another pixel is always at a longer range!

But that's another story, and I still have the moebius strip to prog (I have a solution that works on paper, but I must test that...).

[Antoni Gual]

The thest you suggest is what I have done, in fact I was trying to optimize the solver...

For the Moebius Band, are you planning to switch to a numerical equation solver? Newton method or similar?

[Jark]

The initial equation coefficient are modified by Equa4!

That's when I switched from a, b, c, d to a4, a3, a2, a1, etc...

There were already local variables called a1 & a2, so the test routine was messed up, while the result was OK...

Still, there are some wrong pixels...

[Jark]

The parametric definition of the moebius strip is:
x = (r + l * Cos(Theta/2) * cos(Theta)
y = (r + l * Cos(Theta/2) * sin(Theta)
z= l * Sin (Theta)

Let's define T:
T = r + l * Cos(Theta/2),
i.e.: T = r+ z /Tan(Theta/2)

Since Tan (Theta) = 2 * Tan(Theta/2)/(1 - Tan(Theta/2)^2), then:
y/x = 2 * t / (1-t^2) with Tan(Theta/2) = t
and
T^2 = r^2 + 2 * r * z/t + z^2/t^2 = x^2+ y^2

That can be re-written:
t^2 * y + 2* x* t - y = 0
t^2*(x^2+y^2-r^2) - 2* r* z* t - z^2 = 0

If Alpha = t^2 and Beta = t then:
y * Alpha + 2* x * Beta = y
(x^2+y^2-r^2) * Alpha - 2 * r * z* Beta = z^2

That's a two equations linear system. So:
Alpha = (x* z^2+r * z* y)/(x * (x^2 + y^2 - r^2) + r * y * z)
Beta = (y - y* Alpha)/(2* x)


We're in a raytracing program:
x = f(Dist): y = g(Dist): z = h(Dist)

So we have:
1) Alpha = Phi(Dist)
2) Beta = Psi (Dist)

and we are looking for Dist that matches Alpha = Beta^2. This can be done with a dichotomy solving method.

So we have Dist, so we have x, y, and z...

The normal vector is fairly simple to calculate:
S = +/- SQR (Beta^2/(1+Beta^2))
xn(1) = S (1-S^2)
xn(2) = 1 * S^2 * SQR (1-S^2)
xn(3) = +/- SQR(1-S^2)

The last step is to set a limit to the strip width:

Sin (Theta/2) = +/- SQR(Beta^2/(1-Beta^2))

The limit is ABS(z) <= r * SIN (Theta/2) with Theta = Arg(x,y)


' =================================
' Well, I did this after a couple of coffees waiting for the bus to come... I checked it and checked it and checked it... Finding Alpha and Beta such as Alpha = Beta ^2 is the trick. Hope it will work !

Write that in a natural maths handwriting, you'll see it's not that complicated!

[Jark]

I forgot to post the corrected code !

Code:

SUB Equa4 (a4#, a3#, a2#, a1#, a0#, x1#, x2#, x3#, x4#, Err4%, nSol4%)
' Fourth degree equation resolution - Ferrari method

DIM u#(3)

IF a4# = 0 THEN
Equa3 a3#, a2#, a1#, a0#, x1#, x2#, x3#, Err4%, nSol4%
EXIT SUB
END IF

bdim# = a3# / a4#: cdim# = a2# / a4#: ddim# = a1# / a4#: edim# = a0# / a4#

Fa4# = cdim# - 3 / 8 * bdim# ^ 2
Fa3# = bdim# ^ 3 / 8 - bdim# * cdim# / 2 + ddim#
Fa2# = bdim# ^ 2 * cdim# / 16 - 3 / 256 * bdim# ^ 4 - ddim# * bdim# / 4 + edim#


' BiQuartic case
IF Fa3# = 0 THEN
Equa2 1, Fa4#, Fa2#, x21#, x22#, Err4%, nSol2%

IF Err4% <> 0 THEN
nSol4% = 0
EXIT SUB
END IF

IF nSol2% = 0 THEN
nSol4% = 0
EXIT SUB
END IF

IF x21# < 0 AND x22# < 0 THEN
nSol4% = 0
EXIT SUB
END IF

nSol4% = 4

IF x21# >= 0 THEN
x1# = SQR(x21#) - a3# / a4# / 4
x2# = -x1# - a3# / a4# / 4
END IF

IF x22# >= 0 THEN
x3# = SQR(x22#) - a3# / a4# / 4
x4# = -x3# - a3# / a4# / 4
END IF

IF x21# < 0 AND x22# >= 0 THEN
nSol4% = 2
x1# = x3#: x2# = x4#
x3# = 0: x4# = 0
END IF

IF x22# < 0 AND x21# >= 0 THEN
nSol4% = 2
x3# = 0: x4# = 0
END IF

END IF

IF Fa3# <> 0 THEN

Equa3 1, -Fa4#, -4 * Fa2#, 4 * Fa4# * Fa2# - Fa3# * Fa3#, u#(1), u#(2), u#(3), Err3%, nSol3%
Choice% = 0

FOR k% = 1 TO nSol3%
IF u#(k%) > Fa4# THEN
SF# = SQR(u#(k%) - Fa4#)
ELSE
GOTO Skip:
END IF

d1# = SF# ^ 2 - 4 * (u#(k%) / 2 + Fa3# / 2 / SF#)
d2# = SF# ^ 2 - 4 * (u#(k%) / 2 - Fa3# / 2 / SF#)

IF d1# >= 0 OR d2# >= 0 THEN Choice% = k%: EXIT FOR
Skip:
NEXT k%

IF Choice% = 0 THEN nSol4% = 0: EXIT SUB

U0# = u#(Choice%)
SF# = SQR(U0# - Fa4#)
a11# = SF#
b11# = Fa4# - U0#
C11# = SF# * U0# / 2 + Fa3# / 2
d11# = b1# ^ 2 - 4 * a1# * C1#
a12# = SF#
b12# = U0# - Fa4#
C12# = SF# * U0# / 2 - Fa3# / 2
d12# = b2# ^ 2 - 4 * a2# * C2#
Equa2 a11#, b11#, C11#, fx1#, fx2#, Err21%, nSol21%
Equa2 a12#, b12#, C12#, fx3#, fx4#, Err22%, nSol22%
nSol4% = nSol21% + nSol22%
IF nSol21% = 2 AND nSol22% = 0 THEN
x1# = fx1# - a3# / a4# / 4
x2# = fx2# - a3# / a4# / 4
END IF
IF nSol21% = 0 AND nSol22% = 2 THEN
x1# = fx3# - a3# / a4# / 4
x2# = fx4# - a3# / a4# / 4
END IF
IF nSol21% = 2 AND nSol22% = 2 THEN
x1# = fx1# - a3# / a4# / 4
x2# = fx2# - a3# / a4# / 4
x3# = fx3# - a3# / a4# / 4
x4# = fx4# - a3# / a4# / 4
END IF

END IF

END SUB

As for Equa3:

Code:

SUB Equa3 (a3#, a2#, a1#, a0#, x1#, x2#, x3#, Err3%, nSol3%)
' Third degree equation resolution - Cardan method
Err3% = 0

IF a3# = 0 AND a2# = 0 AND a1# = 0 AND a0# = 0 THEN
nSol3% = 0
Err3% = 1
EXIT SUB
END IF

IF a3# = 0 AND a2# = 0 AND a1# = 0 AND a0# <> 0 THEN
nSol3% = 0
EXIT SUB
END IF

IF a3# = 0 AND a2# = 0 THEN
Equa1 a1#, a0#, x1#, Err3%, nSol3%
EXIT SUB
END IF

IF a3# = 0 THEN
Equa2 a2#, a1#, a0#, x1#, x2#, Err3%, nSol3%
EXIT SUB
END IF

IF a3# <> 0 THEN
p# = a1# / a3# - 1 / 3 * (a2# / a3#) ^ 2
q# = 2 / 27 * (a2# / a3#) ^ 3 + a0# / a3# - a2# / 3 / a3# ^ 2 * a1#
Delta3# = q# ^ 2 + 4 / 27 * p# ^ 3

IF Delta3# > 0 THEN
x1# = Sqr3(-q# / 2 + SQR(Delta3#) / 2) + Sqr3(-q# / 2 - SQR(Delta3#) / 2) - a2# / 3 / a3#
nSol3% = 1
EXIT SUB
END IF

nSol3% = 3

IF Delta3# = 0 THEN
IF p# <> 0 THEN
x1# = q# / p# * 3 - a2# / 3 / a3#
x2# = -SQR(-p# / 3) - a2# / 3 / a3#
x3# = x2#
IF x1# > x2# THEN SWAP x1#, x2#
IF x1# > x3# THEN SWAP x1#, x3#
IF x2# > x3# THEN SWAP x2#, x3#
ELSE
x1# = -a2# / 3 / a3#: x2# = x1#: x3# = x1#
IF x1# > x2# THEN SWAP x1#, x2#
IF x1# > x3# THEN SWAP x1#, x3#
IF x2# > x3# THEN SWAP x2#, x3#
END IF
EXIT SUB
END IF

IF Delta3# < 0 THEN
r# = SQR(-p# ^ 3 / 27)
CosTeta3# = -q# / 2 / r#
Teta3# = Acos(CosTeta3#)

x1# = 2 * SQR(-p# / 3) * COS(Teta3# / 3) - a2# / 3 / a3#
x2# = 2 * SQR(-p# / 3) * COS(Teta3# / 3 + 2 * Pi / 3) - a2# / 3 / a3#
x3# = 2 * SQR(-p# / 3) * COS(Teta3# / 3 + 4 * Pi / 3) - a2# / 3 / a3#
IF x1# > x2# THEN SWAP x1#, x2#
IF x1# > x3# THEN SWAP x1#, x3#
IF x2# > x3# THEN SWAP x2#, x3#
END IF

END IF


END SUB

Equa2 and Equa1 are as follows:

Code:

SUB Equa1 (a1#, a0#, x1#, Err1%, nSol1%)
' First degree equation resolution
IF a1# = 0 AND a0# = 0 THEN Err1% = 1: nSol1% = 0: EXIT SUB
IF a1# = 0 AND a0# <> 0 THEN Err1% = 0: nSol1% = 0: EXIT SUB
IF a1# <> 0 THEN Err1% = 0: x1# = -a0# / a1#: nSol1% = 1: EXIT SUB
END SUB

SUB Equa2 (a2#, a1#, a0#, x1#, x2#, Err2%, nSol2%)
' Second degree equation resolution - Discriminant method
Err2% = 0
IF a2# = 0 AND a1# = 0 AND a0# = 0 THEN Err2% = 1: nSol2% = 0: EXIT SUB
IF a2# = 0 AND a1# = 0 AND a0# <> 0 THEN nSol2% = 0: EXIT SUB

IF a2# = 0 THEN
  Equa1 a1#, a0#, x1#, Err2%, nSol2%
ELSE
  Delta2# = a1# * a1# - 4 * a2# * a0#
  IF Delta2# >= 0 THEN
    nSol2% = 2
    x1# = -(a1# + SQR(Delta2#)) * .5 / a2#
    x2# = (SQR(Delta2#) - a1#) * .5 / a2#
  ELSE
    nSol2% = 0: EXIT SUB
  END IF
END IF
END SUB

The important feature in these routines is the error flag, which is set to "1" if the equation has no meaning, like 0*x = 0 for example.

In a raytracer, that corresponds to probing a plan with a ray that belongs to that plan, or to probing a mathematical cone with a ray that hits the particular point of the cone (the point where the two half cones meet).

[Jark]

I hope there won't be too many corrections to this version :lol:

The equation solvers are corrected, and picture 8 now corresponds to the "Bike Mirrors" picture, since the "Power of Three" became useless...

So the latest version of Tc-Ray is:

http://mandelbrot.dazibao.free.fr/Tcray/TcRay21.zip

[Agamemnus]

Sigh......

Once you have a 100% full working version, edit what you have... that code is way inefficient...... you're duplicating results...... using ANDs........

[na_th_an]

Nah. If you want more speed, talk to me when you finish and I'll help you to port it to C / Allegro.

Amazing work, btw.

[Jark]

When I wrote these routines, my "philosophy" was to prog the maths formulas the way someone would write them, for lisibility.

I have to admit I'm horrified when I look back to these old pieces of code...

Aga : I have saved your posts, and I will analyse that, I promise :wink:

As for today, I have started writing the documentation for TC-Ray, (in english). The heaviest part is to draw sketches that help visualising the way the prog works...