Ensemble de Mandelbrot : [c(z^p - z^q)-1]^2

' **********************************************************************
' Ensemble de Mandelbrot : [c(z^p - z^q) - 1]^r
' **********************************************************************

' Variables definies par l'utilisateur

dim PicWidth%  : PicWidth%  = 640   : ' Taille de l'image
dim PicHeight% : PicHeight% = 480   : '   en pixels
dim p%         : p% = 6             : ' Exposant p
dim q%         : q% = 5             : ' Exposant q
dim r%         : r% = 2             : ' Exposant r
dim x0         : x0 = -15           : ' Coordonnees du centre
dim y0         : y0 = 0             : '   de l'image
dim MaxIter%   : MaxIter% = 200     : ' Nb maxi d'iterations
dim ZoomFact   : ZoomFact = 0.15    : ' Facteur de zoom
dim DistFact   : DistFact = -2      : ' Distance estimee --> V
dim ColorFact  : ColorFact = -3     : ' Continuous dwell --> H, S

' Variables supplementaires

dim HalfPicWidth%  : HalfPicWidth%  = PicWidth% / 2
dim HalfPicHeight% : HalfPicHeight% = PicHeight% / 2

dim Esc   : Esc   = 10000000000  : ' Rayon d'échappement (10^10)
dim Esc2  : Esc2  = Esc * Esc
dim LnEsc : LnEsc = log(Esc)
dim Lnp   : Lnp   = log(p%)

dim Red%, Green%, Blue% : ' Parametres de coloration
dim AbsColor            : ' abs(ColorFact)
dim ScaleFact           : ' Facteur d'echelle = distance entre 2 pixels
dim Nx%, Ny%            : ' Coordonnées d'un point (pixels)
dim xt, yt              : ' Coordonnées d'un point (algébriques)
dim zn_x, zn_y          : ' Nombre complexe z(n)
dim dzn_x, dzn_y        : ' Dérivée dz(n)/dc
dim z0_x, z0_y          : ' Point critique

' --------------------------------------------------------------
' Variables globales de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

' Constantes mathematiques

dim MaxNum, MinNum, MaxLog, MinLog, Pi, PiDiv2

MaxLog =  709.78     : ' Argument max. pour EXP
MinLog = -708.39     : ' Argument min. pour EXP
MaxNum = exp(MaxLog) : ' Nb reel max. ~ 2^1024
MinNum = exp(MinLog) : ' Nb reel min. ~ 2^(-1022)
Pi     = 4 * atn(1)
PiDiv2 = Pi / 2

' Resultats des calculs
' Partie reelle, partie imaginaire, module, argument, signe

dim r_x, r_y, r_mod, r_arg, r_sgn

' Code d'erreur
'  0 = pas d'erreur
' -1 = argument hors bornes
' -2 = singularite
' -3 = overflow
' -4 = underflow

dim ErrCode%

' ----------------------------------------------------------------------
' Description des objets
' ----------------------------------------------------------------------

' Fenêtre principale

left    0, 50
top     0, 50
width   0, PicWidth% + 70
height  0, PicHeight% + 120
caption 0, "Ensemble de Mandelbrot : [c(z^p - z^q) - 1]^r ... Veuillez patienter."

' Image

picture 1
left    1, 30
top     1, 40
width   1, PicWidth%
height  1, PicHeight%

2d_target_is 1

' ----------------------------------------------------------------------
' Programme principal
' ----------------------------------------------------------------------

ColorFact = 0.01 * ColorFact
AbsColor  = abs(ColorFact)
ScaleFact = 4 / (PicHeight% * ZoomFact)

z0_x = power(q% / p%, 1 / (p% - q%))
z0_y = 0

for Ny% = 0 to PicHeight% - 1
 yt = y0 - ScaleFact * (Ny% - HalfPicHeight%)
 for Nx% = 0 to PicWidth% - 1
   xt = x0 + ScaleFact * (Nx% - HalfPicWidth%)
   Mandelbrot(xt, yt)
   2d_pen_color Red%, Green%, Blue%
   2d_point Nx%, Ny%
 next Nx%
next Ny%

caption 0, "Ensemble de Mandelbrot : [c(z^p - z^q) - 1]^r ... Terminé."

file_save 1, "mandel.bmp"

end

' ----------------------------------------------------------------------
' Sous-programmes
' ----------------------------------------------------------------------

sub Func(x, y, cx, cy, dx, dy)
' Calcule le nombre complexe z(n) et la derivee dz(n)/dc
' Resultats dans (zn_x, zn_y) et (dzn_x, dzn_y)

 dim_local zpm1_x, zpm1_y   : ' z^(p-1)
 dim_local zp_x, zp_y       : ' z^p
 dim_local zqm1_x, zqm1_y   : ' z^(q-1)
 dim_local zq_x, zq_y       : ' z^q
 dim_local g_x, g_y         : ' g = z^p - z^q
 dim_local u_x, u_y         : ' c * g - 1
 dim_local urm1_x, urm1_y   : ' u^(r-1)
 dim_local dg_x, dg_y       : ' dg/dc

 CIntPower(x, y, p% - 1)    : zpm1_x = r_x : zpm1_y = r_y
 CMul(x, y, zpm1_x, zpm1_y) : zp_x   = r_x : zp_y   = r_y
 CIntPower(x, y, q% - 1)    : zqm1_x = r_x : zqm1_y = r_y
 CMul(x, y, zqm1_x, zqm1_y) : zq_x   = r_x : zq_y   = r_y

 g_x = zp_x - zq_x
 g_y = zp_y - zq_y

 CMul(cx, cy, g_x, g_y)         : u_x = r_x - 1 : u_y = r_y
 CIntPower(u_x, u_y, r% - 1)    : urm1_x = r_x  : urm1_y  = r_y
 CMul(u_x, u_y, urm1_x, urm1_y) : zn_x = r_x    : zn_y = r_y
 
 dg_x = p% * zpm1_x - q% * zqm1_x
 dg_y = p% * zpm1_y - q% * zqm1_y

 CMul(dx, dy, dg_x, dg_y)         : dg_x  = r_x      : dg_y = r_y
 CMul(cx, cy, dg_x, dg_y)         : dg_x = r_x + g_x : dg_y = r_y + g_y
 CMul(dg_x, dg_y, urm1_x, urm1_y) : dzn_x = r% * r_x : dzn_y = r% * r_y
end_sub

sub Mandelbrot(xt, yt)
' Calcul des iterations au point (xt, yt)

 dim_local iter%, x, y, cx, cy, dx, dy, m2, mz, mdz

 x = z0_x : y = z0_y  : ' z_0 = pt. critique
 dx = 0   : dy = 0    : ' z'_0 = 0
 cx = xt  : cy = yt   : ' c = z_pixel

 m2 = x * x + y * y   : ' |z|^2

 iter% = 0

 while iter% < MaxIter% and m2 < Esc2

   Func(x, y, cx, cy, dx, dy)

   x = zn_x   : y = zn_y
   dx = dzn_x : dy = dzn_y

   m2 = x * x + y * y

   iter% = iter% + 1

 end_while

 if iter% = MaxIter%
   Red% = 255
   Green% = 255
   Blue% = 255
 else
   mz = sqr(m2)
   mdz = sqr(dx * dx + dy * dy)
   MdbCol(iter%, mz, mdz)
 end_if
end_sub

sub MdbCol(iter%, mz, mdz)
' Determine la couleur d'un point --> Resultats dans Red%, Green%, Blue%

 dim_local lmz, Dist, Dwell, D, Q, Angle, Radius, H, S, V

 lmz = log(mz)

 Dist = p% * mz * lmz / mdz

' Determine la luminosite (Value, V) d'apres la distance estimee (Dist)

 D = log(Dist / ScaleFact) / Lnp + DistFact

 if D > 0
   V = 1
 else
   if D > -8
     V = 1 + D / 8
   else
     V = 0
   end_if
 end_if

' Determine la teinte (Hue, H) et la Saturation (S)
' d'apres le "Continuous dwell"

 Dwell = iter% + log(LnEsc / lmz) / Lnp

 Q = log(Dwell) * AbsColor

 if Q < 0.5
   Q = 1 - 1.5 * Q
   Angle = 1 - Q
 else
   Q = 1.5 * Q - 0.5
   Angle = Q
 end_if

 Radius = sqr(Q)

' Si ColorFact > 0, assombrir une bande sur 2

 if (ColorFact > 0) and (odd(iter%) > 0)
   V = 0.85 * V
   Radius = 0.667 * Radius
 end_if

 H = Angle * 10
 H = H - int(H)
 H = H * 360
 
 S = Radius - int(Radius)
 
' Convertir HSV en RGB  

 HSVtoRGB(H, S, V)
end_sub

sub HSVtoRGB(H, S, V)
' Conversion HSV --> RGB. Resultats dans Red%, Green%, Blue%

 dim_local II%, ZZ, FF, PP, QQ, TT, RR, GG, BB

 if S = 0
   R% = int(V * 255)
   G% = R%
   B% = R%
 else
   ZZ  = H / 60
   II% = int(ZZ)
   FF  = ZZ - int(ZZ)
   PP  = V * (1 - S)
   QQ  = V * (1 - S * FF)
   TT  = V * (1 - S * (1 - FF))

   select II%
     case 0
       RR = V  : GG = TT : BB = PP
     case 1
       RR = QQ : GG = V  : BB = PP
     case 2
       RR = PP : GG = V  : BB = TT
     case 3
       RR = PP : GG = QQ : BB = V
     case 4
       RR = TT : GG = PP : BB = V
     case 5
       RR = V  : GG = PP : BB = QQ
   end_select

   Red% = int(RR * 255)
   Green% = int(GG * 255)
   Blue% = int(BB * 255)
 end_if
end_sub

' --------------------------------------------------------------
' Procedures de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

sub CMul(a_x, a_y, b_x, b_y)
' Multiplication : r_x + i r_y = (a_x + i a_y) * (b_x + i b_y)

 ErrCode% = 0

 r_x = a_x * b_x - a_y * b_y
 r_y = a_x * b_y + a_y * b_x
end_sub

sub CSquare(a_x, a_y)
' Carre : r_x + i r_y = (a_x + i a_y)^2

 ErrCode% = 0

 r_x = a_x * a_x - a_y * a_y
 r_y = 2 * a_x * a_y
end_sub

sub CInv(a_x, a_y)
' Inverse : r_x + i r_y = 1 / (a_x + i a_y)

 dim_local Temp

 if a_x = 0 and a_y = 0
   ErrCode% = -3
   r_x = MaxNum
   r_y = MaxNum
 else
   ErrCode% = 0
   Temp = a_x * a_x + a_y * a_y
   r_x = a_x / Temp
   r_y = 0 - a_y / Temp
 end_if
end_sub

sub CIntPower(a_x, a_y, n%)
' Puissance entiere : r_x + i r_y = (a_x + i a_y)^n

 dim_local m%, b_x, b_y, res_x, res_y

 ErrCode% = 0

 if a_x = 0 and a_y = 0
   if n% = 0
     ' 0^0 = lim x^x quand x --> 0 = 1
     r_x = 1
     r_y = 0
   else
     if n% > 0
       ' 0^n = 0 si n > 0
       r_x = 0
       r_y = 0
     else
       ' 0^n indefini si n < 0
       ErrCode% = -2
       r_x = MaxNum
       r_y = MaxNum
     end_if
   end_if
 else
   if n% < 0
     m% = abs(n%)
     CInv(a_x, a_y)
     b_x = r_x
     b_y = r_y
   else
     m% = n%
     b_x = a_x
     b_y = a_y
   end_if

   res_x = 1 : res_y = 0

   while m% > 0
     if odd(m%) = 1
       CMul(b_x, b_y, res_x, res_y)
       res_x = r_x
       res_y = r_y
     end_if
     CSquare(b_x, b_y)
     b_x = r_x
     b_y = r_y
     m% = int(m% / 2)
   end_while

   r_x = res_x
   r_y = res_y
 end_if
end_sub