Ensemble de Mandelbrot z^p+c avec exposant p complexe

' **********************************************************************
' Ensemble de Mandelbrot : z(n) = z(n-1)^p + c  (p complexe)
' **********************************************************************

' Variables definies par l'utilisateur

dim PicWidth%  : PicWidth%  = 640   : ' Taille de l'image
dim PicHeight% : PicHeight% = 480   : '   en pixels
dim p_x        : p_x = 2            : ' Exposant p (partie reelle, > 1)
dim p_y        : p_y = 0.01         : ' Exposant p (partie imaginaire)
dim x0         : x0 = -0.75         : ' Coordonnees du centre
dim y0         : y0 = 0             : '   de l'image
dim MaxIter%   : MaxIter% = 200     : ' Nb maxi d'iterations
dim ZoomFact   : ZoomFact = 1.5     : ' 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_x)

dim flag%         : ' Nature de p (0 : entier, 1 : reel, 2 : complexe)
dim R%, G%, B%    : ' 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

' --------------------------------------------------------------
' 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 : z(n) = z(n-1)^p + c ... Veuillez patienter."

' Image

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

2d_target_is 1

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

Init()

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

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 R%, G%, B%
   2d_point Nx%, Ny%
 next Nx%
next Ny%

caption 0, "Ensemble de Mandelbrot : z(n) = z(n-1)^p + c ... Terminé."

file_save 1, "mandel.bmp"

end

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

sub Init()
' Initialise les variables selon le type de fractale

 if p_y = 0
   if p_x = int(p_x)
     flag% = 0
   else
     flag% = 1
   end_if
 else
   flag% = 2
 end_if
end_sub

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 u_x, u_y

 if flag% = 0
   CIntPower(x, y, int(p_x) - 1)
 else
   if flag% = 1
     CRealPower(x, y, p_x - 1)
   else
     CPower(x, y, p_x - 1, p_y)
   end_if
 end_if
 
 u_x = r_x : u_y = r_y              : ' z^(p-1)
 
 CMul(u_x, u_y, x, y)               : ' z^p
 zn_x = r_x + cx : zn_y = r_y + cy  : ' z^p + c

 CMul(u_x, u_y, dx, dy)             : ' z^(p-1) z'
 dzn_x = r_x : dzn_y = r_y
 
 if flag% = 2
   CMul(p_x, p_y, dzn_x, dzn_y)     : ' p z^(p-1) z'
 else
   dzn_x = p_x * dzn_x
   dzn_y = p_x * dzn_y
 end_if

 dzn_x = dzn_x + 1                  : ' p z^(p-1) z' + 1
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 = 0   : y = 0     : ' z = 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%
   R% = 255
   G% = 255
   B% = 255
 else
   mz = sqr(m2)
   mdz = sqr(dzn_x * dzn_x + dzn_y * dzn_y)
   MdbCol(iter%, mz, mdz)
 end_if
end_sub

sub MdbCol(iter%, mz, mdz)
' Determine la couleur d'un point --> Resultats dans R%, G%, B%

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

 lmz = log(mz)

 Dist = p_x * 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 R%, G%, B%

 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

   R% = int(RR * 255)
   G% = int(GG * 255)
   B% = 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

sub CAbs(a_x, a_y)
' Module : r_mod = |a_x + i a_y|
' Algorithme d'apres "Numerical Recipes"

 ErrCode% = 0

 dim_local AbsX, AbsY, R, C

 AbsX = abs(a_x)
 AbsY = abs(a_y)

 if a_x = 0
   r_mod = abs(a_y)
 else
   if a_y = 0
     r_mod = abs(a_x)
   else
     if AbsX > AbsY
       R = AbsY / AbsX
       C = AbsX
     else
       R = AbsX / AbsY
       C = AbsY
     end_if
     r_mod = C * sqr(1 + R * R)
   end_if
 end_if
end_sub

sub CArg(a_x, a_y)
' Argument : r_arg = arg(a_x + i a_y)
' Resultat dans [-Pi, Pi)
' Equivaut a atan2(a_y, a_x)

 ErrCode% = 0

 if a_x = 0
   r_arg = sgn(a_y) * PiDiv2
 else
   ' 4e / 1er quadrant : -Pi/2..Pi/2
   r_arg = atn(a_y / a_x)
   if a_x < 0
     if a_y > 0
       ' 2e quadrant : Pi/2..Pi
       r_arg = r_arg + Pi
     else
       ' 3e quadrant : -Pi..-Pi/2
       r_arg = r_arg - Pi
     end_if
   end_if
 end_if
end_sub

sub CRealPower(a_x, a_y, p)
' Puissance (exposant reel) : (a_x + i a_y)^p
' Resultat dans r_x, r_y
' Resultat aussi dans r_mod, r_arg si a <> 0

 ErrCode% = 0

 if a_x = 0 and a_y = 0
   if p = 0
     ' 0^0 = lim x^x quand x --> 0 = 1
     r_x = 1
     r_y = 0
   else
     if p > 0
       ' 0^p = 0 si p > 0
       r_x = 0
       r_y = 0
     else
       ' 0^p indefini si p < 0
       ErrCode% = -2
       r_x = MaxNum
       r_y = MaxNum
     end_if
   end_if
 else
   CAbs(a_x, a_y)
   CArg(a_x, a_y)
   r_mod = power(r_mod, p)
   r_arg = r_arg * p
   r_x = r_mod * cos(r_arg)
   r_y = r_mod * sin(r_arg)
 end_if
end_sub

sub CExp(a_x, a_y)
' Exponentielle complexe : r_x + i r_y = exp(a_x + i a_y)

 dim_local ExpX

 if a_x < MinLog
   ErrCode% = -4
   r_x = 0
   r_y = 0
 else
   if a_x > MaxLog
     ErrCode = -3
     ExpX = MaxNum
   else
     ErrCode% = 0
     ExpX = exp(a_x)
   end_if
   r_x = ExpX * cos(a_y)
   r_y = ExpX * sin(a_y)
 end_if
end_sub

sub CPower(a_x, a_y, b_x, b_y)
' Puissance (exposant complexe) : (a_x + i a_y)^(b_x + i b_y)
' Resultat dans r_x, r_y

 ErrCode% = 0

 if a_x = 0 and a_y = 0
   if b_x = 0 and b_y = 0
     ' 0^0 = lim x^x quand x --> 0 = 1
     r_x = 1
     r_y = 0
   else
     ' 0^p = 0 si p > 0
     r_x = 0
     r_y = 0
   end_if
 else
   ' exp(b ln(a))
   CAbs(a_x, a_y)
   CArg(a_x, a_y)
   CMul(b_x, b_y, log(r_mod), r_arg)
   CExp(r_x, r_y)
 end_if
end_sub