Figures fractales

' **********************************************************************
' Ensembles de Mandelbrot/Julia : f(z) = z * sinh(z) - c^2 (Catalisator)
' **********************************************************************

' Parametres de base (a modifier eventuellement)

const PicWidth  = 600       ' Largeur de l'image
const PicHeight = 400       ' Hauteur de l'image
const InsideCol = CL_BLANC  ' Couleur de la partie interne

dim Nom$       ' Nom de l'image
dim x0         ' Coord. X du centre
dim y0         ' Coord. Y du centre
dim MaxIter%   ' Nb maxi d'iterations
dim ZoomFact   ' Facteur de zoom
dim DistFact   ' Plus negatif ==> contour plus sombre
dim ColorFact  ' 0 pour noir et blanc, positif pour bandes de couleur
dim cx, cy     ' Constante c pour l'ensemble de Julia
dim Julia&     ' Ensemble de Julia ?

' ----------------------------------------------------------------------
' Exemple
' ----------------------------------------------------------------------

'              Nom       x0 y0  maxit zoom  dist col cx cy
init_params "MandelCat",  0, 0,   200,   2, -0.5, -1, 0, 0

' **********************************************************************

' Parametres supplementaires

const HalfPicWidth  = PicWidth \ 2
const HalfPicHeight = PicHeight \ 2
 
const Esc   = 1.0E+10  ' Rayon d'echappement
const LnEsc = log(Esc)

const p   = 2
const Lnp = log(2)     ' Pour la coloration

dim ColFact, AbsCol, ScaleFact
 
SetParams()

' **********************************************************************

' Programme principal

dim x%, y%, btn%, key$

mode 3, "Clic gauche : zoom + / Clic droit : zoom - / S : sauvegarde / ESC : quitte", PicWidth, PicHeight
  
colormap PicWidth, PicHeight, fadr(mandelbrot)

repeat
  get_mouse x, y, btn

  if btn = 1 or btn = 2 then
    getcoord x0, y0, x, y

    SetZoom(btn = 1)
    SetParams()

    colormap PicWidth, PicHeight, fadr(mandelbrot)
  end_if

  key = inkey()
  if upper(key) = "S" then save()
until key = "ESCAPE"

' **********************************************************************

' Sous-programmes

sub init_params(_Nom$, _x0, _y0, _MaxIter%, _ZoomFact, _DistFact, _ColorFact, _cx, _cy)
' Initialise les parametres

  Nom       = _Nom
  x0        = _x0
  y0        = _y0
  MaxIter   = _MaxIter
  ZoomFact  = _ZoomFact
  DistFact  = _DistFact
  ColorFact = _ColorFact
  cx        = _cx
  cy        = _cy

  Julia = (cx <> 0) or (cy <> 0)
end_sub

sub init_sub (xt, yt, c@, z@, dz@)
' Initialise les iterations au point (xt, yt)

  dim r

  if Julia then
    c  = cmplx(cx, cy)
    z  = cmplx(xt, yt)
    dz = cmplx(-2, 0)
  else
    c  = cmplx(xt, yt)
    z  = C_zero
    dz = C_zero
  end_if
end_sub

sub iter_sub (c@, z@, dz@, zn@, dzn@)
' Calcul d'une iteration

  dim sinhz@, coshz@

  sinhz = csinh(z)
  coshz = ccosh(z)

  zn = z * sinhz - c * c
  dzn = dz * (z * coshz + sinhz)
  
  if not Julia then dzn = dzn - 2 * c
end_sub

sub SetParams ()
  ColFact = 0.01 * ColorFact
  AbsCol = abs(ColFact)
  ScaleFact = 4 / (PicHeight * ZoomFact)
end_sub

sub SetZoom (zoom%)
  if zoom then
    ZoomFact = 5 * ZoomFact
    MaxIter = 1.25 * MaxIter
    ColorFact = 1.1 * ColorFact
  else
    ZoomFact = ZoomFact / 5
    MaxIter = MaxIter / 1.25
    ColorFact = ColorFact / 1.1
  end_if
end_sub

function rgbcol% (iter%, a, v)
' Determine la teinte (Hue) et la saturation d'apres le "Continuous Dwell"

  dim angle, radius, h, s, rr%, gg%, bb%

  if a < 0.5 then
    a = 1 - 1.5 * a
    angle = 1 - a
  else
    a = 1.5 * a - 0.5
    angle = a
  end_if

  radius = sqr(a)

  if (ColFact > 0) and odd(iter) then
    v = 0.85 * v
    radius = 0.667 * radius
  end if

  h = frac(angle * 10) * 360
  s = frac(radius)

  HSVtoRGB h, s, v, rr, gg, bb
  return RGB(rr, gg, bb)
end_function

function mdbcol% (iter%, mz, mdz)
' Coloration pour les ensembles de Mandelbrot/Julia
' Retourne la couleur RGB d'un point en fonction de :
'   iter = nb d'iterations
'   mz, mdz = modules de z et de (dz/dc) a l'iteration iter
' Methode de coloration d'apres R. Munafo (http://mrob.com/pub/muency/color.html)

  dim lmz, dist, dwell, dscale, a, v

  ' Determiner la luminosite (Value) d'apres l'estimateur de distance

  lmz = log(mz)

  dist = p * mz * lmz / mdz
  dscale = log(dist / ScaleFact) / Lnp + DistFact

  if dscale > 0 then
    v = 1
  elseif dscale > -8 then
    v = 1 + dscale / 8
  else
    v = 0
  end_if
 
  dwell = iter + log(LnEsc / lmz) / Lnp
  a = log(abs(dwell)) * AbsCol
  return rgbcol(iter, a, v)
end_function

function mandelbrot% (x%, y%)
' Iteration de la fonction complexe au point (x, y)

  dim iter%
  dim c@, z@, dz@, zn@, dzn@
  dim xt, yt, mz, mdz

  xt = x0 + ScaleFact * (x - HalfPicWidth)
  yt = y0 + ScaleFact * (y - HalfPicHeight)

  init_sub xt, yt, c, z, dz

  iter = 0
  mz = cabs(z)

  while iter < MaxIter and mz < Esc
    iter_sub c, z, dz, zn, dzn
    z = zn
    dz = dzn
    mz = cabs(z)
    iter = iter + 1
  wend

  if iter = MaxIter then return InsideCol

  mdz = cabs(dz)
  return mdbcol(iter, mz, mdz)
end_function

sub getcoord (x0, y0, x%, y%)
' Calcule les coordonnees (x0,y0) du centre
' d'apres la position (x,y) de la souris

  x0 = x0 + ScaleFact * (x - HalfPicWidth)
  y0 = y0 + ScaleFact * (y - HalfPicHeight)
end_sub

sub save ()
' Sauvegarde l'image et les parametres

  img_save Nom + ".png"
  openout #1, Nom + ".cat"
  write #1, Nom, x0, y0, MaxIter, ZoomFact, DistFact, ColorFact, cx, cy
  closeout #1
end_sub

init_params "Septagon",6,1,0,0,0,0,3000,1.75,-2,-3,-0.14,0

' **********************************************************************
' Ce programme trace des ensembles de Julia situes sur un
' cercle entourant l'ensemble de Mandelbrot
' ----------------------------------------------------------------------
' D'apres un programme de Papydall :
' https://panoramic.1fr1.net/t6918-animation-sur-l-ensemble-de-julia#82546
' **********************************************************************

' Parametres de base (a modifier eventuellement)

const Rho       = 1         ' Rayon du cercle
const Delta     = Pi / 180  ' Angle unitaire
const PicWidth  = 500       ' Largeur de l'image
const PicHeight = 500       ' Hauteur de l'image
const InsideCol = CL_BLANC  ' Couleur de la partie interne
const x0        = 0         ' Coord. X du centre
const y0        = 0         ' Coord. Y du centre
const MaxIter   = 1000      ' Nb maxi d'iterations
const ZoomFact  = 1.2       ' Facteur de zoom
const DistFact  = -2        ' Plus negatif ==> contour plus sombre
const ColorFact = 3         ' 0 pour noir et blanc, positif pour bandes de couleur

' **********************************************************************

' Parametres supplementaires

const HalfPicWidth  = PicWidth \ 2
const HalfPicHeight = PicHeight \ 2
 
const Esc   = 1.0E+10  ' Rayon d'echappement
const LnEsc = log(Esc)
const Ln2   = log(2)

dim ColFact, AbsCol, ScaleFact

' **********************************************************************

' Programme principal

dim Nom$, i%, theta, cx, cy

mode 3, "Creation des ensembles de Julia", PicWidth, PicHeight

for i = 0 to 359
  Nom = "julia_" & dec(i, "000")
  theta = i * Delta
  cx = Rho * cos(theta)
  cy = Rho * sin(theta)
  SetParams()
  ? "Creation de l'ensemble "; Nom
  colormap PicWidth, PicHeight, fadr(mandelbrot)
  img_save Nom + ".png"
next i

' **********************************************************************

sub SetParams ()
  ColFact = 0.01 * ColorFact
  AbsCol = abs(ColFact)
  ScaleFact = 4 / (PicHeight * ZoomFact)
end_sub

sub init_sub (xt, yt, c@, z@, dz@)
' Initialise les iterations au point (xt, yt)

  c  = cmplx(cx, cy)
  z  = cmplx(xt, yt)
  dz = C_one
end_sub

sub iter_sub (c@, z@, dz@, zn@, dzn@)
' Calcul d'une iteration

  zn  = z * z + c   ' z^2 + c
  dzn = 2 * z * dz  ' Derivee dz/dc
end_sub

function rgbcol% (iter%, q, v)
' Determine la teinte (Hue) et la saturation d'apres le "Continuous Dwell"

  dim angle, radius, h, s, r%, g%, b%

  if q < 0.5 then
    q = 1 - 1.5 * q
    angle = 1 - q
  else
    q = 1.5 * q - 0.5
    angle = q
  end_if

  radius = sqr(q)

  if (ColFact > 0) and odd(iter) then
    v = 0.85 * v
    radius = 0.667 * radius
  end if

  h = frac(angle * 10) * 360
  s = frac(radius)

  HSVtoRGB h, s, v, r, g, b
  return RGB(r, g, b)
end_function

function mdbcol% (iter%, mz, mdz)
' Coloration pour les ensembles de Mandelbrot/Julia
' Retourne la couleur RGB d'un point en fonction de :
'   iter = nb d'iterations
'   mz, mdz = modules de z et de (dz/dc) a l'iteration iter
' Methode de coloration d'apres R. Munafo (http://mrob.com/pub/muency/color.html)

  dim lmz, dist, dwell, dscale, q, v

  ' Determiner la luminosite (Value) d'apres l'estimateur de distance

  lmz = log(mz)

  dist = 2 * mz * lmz / mdz
  dscale = log(dist / ScaleFact) / Ln2 + DistFact

  if dscale > 0 then
    v = 1
  elseif dscale > -8 then
    v = 1 + dscale / 8
  else
    v = 0
  end_if
 
  dwell = iter + log(LnEsc / lmz) / Ln2
  q = log(abs(dwell)) * AbsCol
  return rgbcol(iter, q, v)
end_function

function mandelbrot% (x%, y%)
' Iteration de la fonction complexe au point (x, y)

  dim iter%
  dim c@, z@, dz@, zn@, dzn@
  dim xt, yt, mz, mdz

  xt = x0 + ScaleFact * (x% - HalfPicWidth)
  yt = y0 + ScaleFact * (y% - HalfPicHeight)

  init_sub xt, yt, c, z, dz

  iter = 0
  mz = cabs(z)

  while iter < MaxIter and mz < Esc
    iter_sub c, z, dz, zn, dzn
    z = zn
    dz = dzn
    mz = cabs(z)
    iter = iter + 1
  wend

  if iter = MaxIter then return InsideCol

  mdz = cabs(dz)
  return mdbcol(iter, mz, mdz)
end_function

' **********************************************************************
' Ce programme trace des ensembles de Julia situes sur un cercle
' entourant le principal "bourgeon" de l'ensemble de Mandelbrot
' ----------------------------------------------------------------------
' D'apres un programme de Papydall :
' https://panoramic.1fr1.net/t6918-animation-sur-l-ensemble-de-julia#82546
' **********************************************************************

' Parametres de base (a modifier eventuellement)

const Xc        = -1        ' Centre du cercle en (Xc, Yc)
const Yc        = 0
const Rho       = 0.27      ' Rayon du cercle
const Delta     = Pi / 180  ' Angle unitaire
const NPics     = 180       ' Nb d'images a tracer
const PicWidth  = 512       ' Largeur de l'image
const PicHeight = 400       ' Hauteur de l'image
const InsideCol = CL_BLANC  ' Couleur de la partie interne
const MaxIter   = 1000      ' Nb maxi d'iterations
const ZoomFact  = 1.5       ' Facteur de zoom
const DistFact  = -1        ' Plus negatif ==> contour plus sombre
const ColorFact = 3         ' 0 pour noir et blanc, positif pour bandes de couleur

' **********************************************************************

' Parametres supplementaires

const HalfPicWidth  = PicWidth \ 2
const HalfPicHeight = PicHeight \ 2

const x0    = 0        ' Centre de l'ensemble de Julia
const y0    = 0        
 
const Esc   = 1.0E+10  ' Rayon d'echappement
const LnEsc = log(Esc)
const Ln2   = log(2)

const Fmt   = "0.000000" ' Format d'affichage des nombres

dim ColFact, AbsCol, ScaleFact

' **********************************************************************

' Programme principal

dim Nom$, i%, theta, cx, cy

mode 3, "Creation des ensembles de Julia", PicWidth, PicHeight

for i = 0 to NPics
  Nom = dec(i, "000")
  theta = i * Delta
  cx = Xc + Rho * cos(theta)
  cy = Yc + Rho * sin(theta)
  SetParams()
  colormap PicWidth, PicHeight, fadr(mandelbrot)
  print Nom & ": X = " & dec(cx, Fmt) & ", Y = " &  dec(cy, Fmt)  
  img_save Nom + ".png"
next i

' **********************************************************************

sub SetParams ()
  ColFact = 0.01 * ColorFact
  AbsCol = abs(ColFact)
  ScaleFact = 4 / (PicHeight * ZoomFact)
end_sub

sub init_sub (xt, yt, c@, z@, dz@)
' Initialise les iterations au point (xt, yt)

  c  = cmplx(cx, cy)
  z  = cmplx(xt, yt)
  dz = C_one
end_sub

sub iter_sub (c@, z@, dz@, zn@, dzn@)
' Calcul d'une iteration

  zn  = z * z + c   ' z^2 + c
  dzn = 2 * z * dz  ' Derivee dz/dc
end_sub

function rgbcol% (iter%, q, v)
' Determine la teinte (Hue) et la saturation d'apres le "Continuous Dwell"

  dim angle, radius, h, s, r%, g%, b%

  if q < 0.5 then
    q = 1 - 1.5 * q
    angle = 1 - q
  else
    q = 1.5 * q - 0.5
    angle = q
  end_if

  radius = sqr(q)

  if (ColFact > 0) and odd(iter) then
    v = 0.85 * v
    radius = 0.667 * radius
  end if

  h = frac(angle * 10) * 360
  s = frac(radius)

  HSVtoRGB h, s, v, r, g, b
  return RGB(r, g, b)
end_function

function mdbcol% (iter%, mz, mdz)
' Coloration pour les ensembles de Mandelbrot/Julia
' Retourne la couleur RGB d'un point en fonction de :
'   iter = nb d'iterations
'   mz, mdz = modules de z et de (dz/dc) a l'iteration iter
' Methode de coloration d'apres R. Munafo (http://mrob.com/pub/muency/color.html)

  dim lmz, dist, dwell, dscale, q, v

  ' Determiner la luminosite (Value) d'apres l'estimateur de distance

  lmz = log(mz)

  dist = 2 * mz * lmz / mdz
  dscale = log(dist / ScaleFact) / Ln2 + DistFact

  if dscale > 0 then
    v = 1
  elseif dscale > -8 then
    v = 1 + dscale / 8
  else
    v = 0
  end_if
 
  dwell = iter + log(LnEsc / lmz) / Ln2
  q = log(abs(dwell)) * AbsCol
  return rgbcol(iter, q, v)
end_function

function mandelbrot% (x%, y%)
' Iteration de la fonction complexe au point (x, y)

  dim iter%
  dim c@, z@, dz@, zn@, dzn@
  dim xt, yt, mz, mdz

  xt = x0 + ScaleFact * (x% - HalfPicWidth)
  yt = y0 + ScaleFact * (y% - HalfPicHeight)

  init_sub xt, yt, c, z, dz

  iter = 0
  mz = cabs(z)

  while iter < MaxIter and mz < Esc
    iter_sub c, z, dz, zn, dzn
    z = zn
    dz = dzn
    mz = cabs(z)
    iter = iter + 1
  wend

  if iter = MaxIter then return InsideCol

  mdz = cabs(dz)
  return mdbcol(iter, mz, mdz)
end_function

init_params "mandel_2", -1.768978635265551, 0.00460540488510505, 10000, 6.0E+11, -1, 5.7

' **********************************************************************
' Ensemble de Mandelbrot : f(z) = z^2 + c
' Calcul en multiprecision avec la bibliotheque MPFR
' **********************************************************************

' Parametres de base (a modifier eventuellement)

set_mpfr_prec 30           ' Nombre de chiffres

const PicWidth  = 500       ' Largeur de l'image
const PicHeight = 500       ' Hauteur de l'image
const InsideCol = CL_BLANC  ' Couleur de la partie interne

dim Nom$       ' Nom de l'image
dim x0!        ' Coord. X du centre
dim y0!        ' Coord. Y du centre
dim ZoomFact!  ' Facteur de zoom
dim DistFact   ' Plus negatif ==> contour plus sombre
dim ColorFact  ' 0 pour noir et blanc, positif pour bandes de couleur
dim MaxIter%   ' Nb maxi d'iterations

' ----------------------------------------------------------------------
' Initialisation a partir d'un fichier de parametres ou directement.
' Les variables en multiprecision sont saisies sous forme de chaines.
' ----------------------------------------------------------------------

Nom = "mandel"  ' Nom du fichier

if file_exists(Nom + ".par") then
  ReadParams()
else
  x0        = "-1.7689786352655529315388638242"
  y0        = " 0.0046054048851050100472176703297"
  ZoomFact  = " 2.5E+20"
  DistFact  = -1
  ColorFact = 11
  MaxIter   = 12500
end_if

' **********************************************************************

' Parametres supplementaires

const HalfPicWidth  = PicWidth \ 2
const HalfPicHeight = PicHeight \ 2
 
dim Esc!, LnEsc!, Ln2!

Esc   = mpfr(1.0E+10)  ' Rayon d'echappement
LnEsc = log(Esc)
Ln2   = log2_const()

dim ColFact, AbsCol, ScaleFact!
 
SetParams()
PrintParams()

' **********************************************************************

' Programme principal

dim x%, y%, btn%, key$, t1, t2

mode 3, "Clic gauche : zoom + / Clic droit : zoom - / S : sauvegarde / ESC : quitte", PicWidth, PicHeight

' Si l'image existe, la charger directement, sinon la creer

if file_exists(Nom + ".png") then
  img_load Nom + ".png"
else
  PlotFractal()
end_if

repeat
  get_mouse x, y, btn
  if btn = 1 or btn = 2 then
    getcoord x0, y0, x, y
    SetZoom(btn = 1)
    SetParams()
    PrintParams()
    PlotFractal()
  end_if
  key = inkey()
  if upper(key) = "S" then save()
until key = "ESCAPE"

' **********************************************************************

' Sous-programmes

sub init_sub (xt!, yt!, c@, z@, dz@)
' Initialise les iterations au point (xt, yt)

  c  = cmplx(xt, yt)
  z  = C_zero
  dz = C_zero
end_sub

sub iter_sub (c@, z@, dz@, zn@, dzn@)
' Calcul d'une iteration

  zn  = z * z + c       ' z^2 + c
  dzn = 2 * z * dz + 1  ' Derivee dz/dc
end_sub

sub SetParams ()
  ColFact = 0.01 * ColorFact
  AbsCol = abs(ColFact)
  ScaleFact = 4 / (PicHeight * ZoomFact)
end_sub

sub SetZoom (zoom%)
  if zoom then
    ZoomFact = 5 * ZoomFact
    MaxIter = 1.25 * MaxIter
    ColorFact = 1.1 * ColorFact
  else
    ZoomFact = ZoomFact / 5
    MaxIter = MaxIter / 1.25
    ColorFact = ColorFact / 1.1
  end_if
  ColorFact = round(ColorFact, 1)
end_sub

function rgbcol% (iter%, q, v)
' Determine la teinte (Hue) et la saturation d'apres le "Continuous Dwell"

  dim angle, radius, h, s
  dim r%, g%, b%

  if q < 0.5 then
    q = 1 - 1.5 * q
    angle = 1 - q
  else
    q = 1.5 * q - 0.5
    angle = q
  end_if

  radius = sqr(q)

  if (ColFact > 0) and odd(iter) then
    v = 0.85 * v
    radius = 0.667 * radius
  end if

  h = frac(angle * 10) * 360
  s = frac(radius)

  HSVtoRGB h, s, v, r, g, b
  return RGB(r, g, b)
end_function

function mdbcol% (iter%, mz!, mdz!)
' Coloration pour les ensembles de Mandelbrot/Julia
' Retourne la couleur RGB d'un point en fonction de :
'   iter = nb d'iterations
'   mz, mdz = modules de z et de (dz/dc) a l'iteration iter
' Methode de coloration d'apres R. Munafo (http://mrob.com/pub/muency/color.html)

  dim lmz!, dist!, dwell!, dscale!, q, v

  ' Determiner la luminosite (Value) d'apres l'estimateur de distance

  lmz = log(mz)

  dist = 2 * mz * lmz / mdz
  dscale = log(dist / ScaleFact) / Ln2 + DistFact

  if dscale > 0 then
    v = 1
  elseif dscale > -8 then
    v = 1 + dscale / 8
  else
    v = 0
  end_if
 
  dwell = iter + log(LnEsc / lmz) / Ln2
  q = log(abs(dwell)) * AbsCol
  return rgbcol(iter, q, v)
end_function

function mandelbrot% (x%, y%)
' Iteration de la fonction complexe au point (x, y)

  dim iter%
  dim c@, z@, dz@, zn@, dzn@
  dim xt!, yt!, mz!, mdz!

  xt = x0 + ScaleFact * (x - HalfPicWidth)
  yt = y0 + ScaleFact * (y - HalfPicHeight)

  init_sub xt, yt, c, z, dz

  iter = 0
  mz = cabs(z)

  while iter < MaxIter and mz < Esc
    iter_sub c, z, dz, zn, dzn
    z = zn
    dz = dzn
    mz = cabs(z)
    iter = iter + 1
  wend

  if iter = MaxIter then return InsideCol

  mdz = cabs(dz)
  return mdbcol(iter, mz, mdz)
end_function

sub PlotFractal ()
' Dessine l'image et calcule le temps
  
  dim t1, t2, t
  
  t1 = timer  ' secondes
  
  colormap PicWidth, PicHeight, fadr(mandelbrot)

  t2 = timer             ' secondes
  t = (t2 - t1) / 86400  ' jours

  ? "Temps de calcul = "; ftime(frac(t))
  ?
end_sub

sub getcoord (x0!, y0!, x%, y%)
' Calcule les coordonnees (x0,y0) du centre
' d'apres la position (x,y) de la souris

  x0 = x0 + ScaleFact * (x - HalfPicWidth)
  y0 = y0 + ScaleFact * (y - HalfPicHeight)
end_sub

sub ReadParams ()
' Lit les parametres dans un fichier .par
' Les parametres en multiprecision sont
' lus sous forme de chaines de caracteres

  dim a$

  openin #1, Nom + ".par"

  input #1, a : x0 = a
  input #1, a : y0 = a
  input #1, a : ZoomFact = a
  input #1, DistFact
  input #1, ColorFact
  input #1, MaxIter
end_sub

sub PrintParams ()
' Ecrit la liste des parametres
' Le point d'interrogation ecrit sur la console

  ? "Nom       = """ & Nom & """"
  ? "x0        = """ & x0  & """"
  ? "y0        = """ & y0  & """"
  
  ? "ZoomFact  = """ & mpfr_to_str(ZoomFact, "%4.2E") & """"

  ? "DistFact  = "   & DistFact  
  ? "ColorFact = "   & ColorFact
  ? "MaxIter   = "   & MaxIter
  ?
end_sub

sub save ()
' Sauvegarde l'image et les parametres

  img_save Nom + ".png"
  openout #1, Nom + ".par"
  
  write #1, x0
  write #1, y0
  write #1, mpfr_to_str(ZoomFact, "%4.2E")
  write #1, DistFact
  write #1, ColorFact
  write #1, MaxIter

  closeout #1
end_sub

' **********************************************************************
' Ensemble de Mandelbrot : f(z) = z^2 + c
' Calcul en multiprecision avec la bibliotheque MPFR
' **********************************************************************

' Parametres de base (a modifier eventuellement)

set_mpfr_prec 50           ' Nombre de chiffres

const PicWidth  = 500       ' Largeur de l'image
const PicHeight = 500       ' Hauteur de l'image
const InsideCol = CL_BLANC  ' Couleur de la partie interne

dim Nom$       ' Nom de l'image
dim x0!        ' Coord. X du centre
dim y0!        ' Coord. Y du centre
dim ZoomFact!  ' Facteur de zoom
dim DistFact   ' Plus negatif ==> contour plus sombre
dim ColorFact  ' 0 pour noir et blanc, positif pour bandes de couleur
dim MaxIter%   ' Nb maxi d'iterations

' ----------------------------------------------------------------------
' Initialisation a partir d'un fichier de parametres ou directement.
' Les variables en multiprecision sont saisies sous forme de chaines.
' ----------------------------------------------------------------------

Nom = "mandel"  ' Nom du fichier

if file_exists(Nom + ".par") then
  ReadParams()
else
  x0        = "-1.768978635265552931538862840826883865374022698613"
  y0        = " 0.004605404885105010047268484803608948913661202185792"
  ZoomFact  = " 2.38E+35"
  DistFact  = -1
  ColorFact = 11
  MaxIter   = 12500
end_if

' **********************************************************************

' Parametres supplementaires

const HalfPicWidth  = PicWidth \ 2
const HalfPicHeight = PicHeight \ 2
 
dim Esc!, LnEsc!, Ln2!

Esc   = mpfr(1.0E+10)  ' Rayon d'echappement
LnEsc = log(Esc)
Ln2   = log2_const()

dim ColFact, AbsCol, ScaleFact!
 
SetParams()
PrintParams()

' **********************************************************************

' Programme principal

dim x%, y%, btn%, key$

mode 3, "Clic gauche : zoom + / Clic droit : zoom - / S : sauvegarde / ESC : quitte", PicWidth, PicHeight

' Si l'image existe, la charger directement, sinon la creer

if file_exists(Nom + ".png") then
  img_load Nom + ".png"
else
  PlotFractal()
end_if

repeat
  get_mouse x, y, btn
  if btn = 1 or btn = 2 then
    getcoord x0, y0, x, y
    SetZoom(btn = 1)
    SetParams()
    PrintParams()
    PlotFractal()
  end_if
  key = inkey()
  if upper(key) = "S" then save()
until key = "ESCAPE"

' **********************************************************************

' Sous-programmes

sub init_sub (xt!, yt!, c@, z@, dz@)
' Initialise les iterations au point (xt, yt)

  c  = cmplx(xt, yt)
  z  = C_zero
  dz = C_zero
end_sub

sub iter_sub (c@, z@, dz@, zn@, dzn@)
' Calcul d'une iteration

  zn  = z * z + c       ' z^2 + c
  dzn = 2 * z * dz + 1  ' Derivee dz/dc
end_sub

sub SetParams ()
  ColFact = 0.01 * ColorFact
  AbsCol = abs(ColFact)
  ScaleFact = 4 / (PicHeight * ZoomFact)
end_sub

sub SetZoom (zoom%)
  dim s = sgn(ColorFact)
  if zoom then
    ZoomFact = 5 * ZoomFact
    MaxIter = 1.25 * MaxIter
    ColorFact = ColorFact + 0.25 * s
  else
    ZoomFact = ZoomFact / 5
    MaxIter = MaxIter / 1.25
    ColorFact = ColorFact - 0.25 * s
  end_if
end_sub

function rgbcol% (iter%, q, v)
' Determine la teinte (Hue) et la saturation d'apres le "Continuous Dwell"

  dim angle, radius, h, s
  dim r%, g%, b%

  if q < 0.5 then
    q = 1 - 1.5 * q
    angle = 1 - q
  else
    q = 1.5 * q - 0.5
    angle = q
  end_if

  radius = sqr(q)

  if (ColFact > 0) and odd(iter) then
    v = 0.85 * v
    radius = 0.667 * radius
  end if

  h = frac(angle * 10) * 360
  s = frac(radius)

  HSVtoRGB h, s, v, r, g, b
  return RGB(r, g, b)
end_function

function mdbcol% (iter%, mz!, mdz!)
' Coloration pour les ensembles de Mandelbrot/Julia
' Retourne la couleur RGB d'un point en fonction de :
'   iter = nb d'iterations
'   mz, mdz = modules de z et de (dz/dc) a l'iteration iter
' Methode de coloration d'apres R. Munafo (http://mrob.com/pub/muency/color.html)

  dim lmz!, dist!, dwell!, dscale!, q, v

  ' Determiner la luminosite (Value) d'apres l'estimateur de distance

  lmz = log(mz)

  dist = 2 * mz * lmz / mdz
  dscale = log(dist / ScaleFact) / Ln2 + DistFact

  if dscale > 0 then
    v = 1
  elseif dscale > -8 then
    v = 1 + dscale / 8
  else
    v = 0
  end_if
 
  dwell = iter + log(LnEsc / lmz) / Ln2
  q = log(abs(dwell)) * AbsCol
  return rgbcol(iter, q, v)
end_function

function mandelbrot% (x%, y%)
' Iteration de la fonction complexe au point (x, y)

  dim iter%
  dim c@, z@, dz@, zn@, dzn@
  dim xt!, yt!, mz!, mdz!

  xt = x0 + ScaleFact * (x - HalfPicWidth)
  yt = y0 + ScaleFact * (y - HalfPicHeight)

  init_sub xt, yt, c, z, dz

  iter = 0
  mz = cabs(z)

  while iter < MaxIter and mz < Esc
    iter_sub c, z, dz, zn, dzn
    z = zn
    dz = dzn
    mz = cabs(z)
    iter = iter + 1
  wend

  if iter = MaxIter then return InsideCol

  mdz = cabs(dz)
  return mdbcol(iter, mz, mdz)
end_function

sub PlotFractal ()
' Dessine l'image et calcule le temps
  
  dim t1, t2, t
  
  t1 = timer  ' secondes
  
  colormap PicWidth, PicHeight, fadr(mandelbrot)

  t2 = timer             ' secondes
  t = (t2 - t1) / 86400  ' jours

  ? "Temps de calcul = "; ftime(frac(t))
  ?
end_sub

sub getcoord (x0!, y0!, x%, y%)
' Calcule les coordonnees (x0,y0) du centre
' d'apres la position (x,y) de la souris

  x0 = x0 + ScaleFact * (x - HalfPicWidth)
  y0 = y0 + ScaleFact * (y - HalfPicHeight)
end_sub

sub ReadParams ()
' Lit les parametres dans un fichier .par
' Les parametres en multiprecision sont
' lus sous forme de chaines de caracteres

  dim a$

  openin #1, Nom + ".par"

  input #1, a : x0 = a
  input #1, a : y0 = a
  input #1, a : ZoomFact = a
  input #1, DistFact
  input #1, ColorFact
  input #1, MaxIter

  closein #1
end_sub

sub PrintParams ()
' Ecrit la liste des parametres
' Le point d'interrogation ecrit sur la console

  ? "Nom       = """ & Nom & """"
  ? "x0        = """ & x0  & """"
  ? "y0        = """ & y0  & """"
  
  ? "ZoomFact  = """ & mpfr_to_str(ZoomFact, "%4.2E") & """"

  ? "DistFact  = "   & DistFact  
  ? "ColorFact = "   & ColorFact
  ? "MaxIter   = "   & MaxIter
  ?
end_sub

sub save ()
' Sauvegarde l'image et les parametres

  img_save Nom + ".png"
  openout #1, Nom + ".par"
  
  write #1, x0
  write #1, y0
  write #1, mpfr_to_str(ZoomFact, "%4.2E")
  write #1, DistFact
  write #1, ColorFact
  write #1, MaxIter

  closeout #1
end_sub

x0 = "-1.768634119157166907234009209392062"
y0 = " 0.002663304521717732927981067912799"
Zoomfact = " 1.70E+30"
DistFact = -1
ColorFact = 10
MaxIter = 20000