UF真是复分形的宝藏

榕坚

UF中的复分形还有很多没被发现，其实包括之前的DEM法也都在UF的范例中。今天看了一下还有这个东西呢：

柳烟

这种分形不知姓甚名谁。

榕坚

在mt系列：
mt-latoocarfian {
; Mark Townsend, 18 August 2001
; Strange attractors from Clifford Pickover's
; book "Chaos in Wonderland" (the "De Jong"
; mutation uses similar equations from A. K.
; Dewdney's "The Armchair Universe".)
global:
  float max = 0
  int i = 0
  float xx = 0, float yy = 0
  float a = 0, float b = 0
  float c = 0, float d = 0
  int seed = @seed
  if !@rnd
    a = @a
    b = @b
    c = @c
    d = @d
  else
  ; We randomize the parameters in a loop
  ; until we find a set that has a reasonably
  ; high Lyapunov exponent.
    float L = -1
    while L < @lthresh
    ; Parameters in the range (-3 < a, b < 3)
    ; and (0.5 < c, d < 1.5)
      seed = random(seed)
      a = (seed / #randomrange) * 3
      seed = random(seed)
      b = (seed / #randomrange) * 3
      seed = random(seed)
      c = (seed / #randomrange) + 0.5
      seed = random(seed)
      d = (seed / #randomrange) + 0.5
      float Lsum = 0
      float lx = 0.1
      float ly = 0.1
      float xe = lx + 0.000001
      float ye = ly
      i = 0
      while i < 1000
      ; Pickover suggests 10 million iterations,
      ; but this is enough to produce decent
      ; attractors
        if @mutation == "Standard"
          xx = sin(ly * b) + c * sin(lx * b)
          yy = (sin(lx * a) + d * sin(ly * a))
        elseif  @mutation == "Alpha"
          xx = sin(ly * b) + sin(lx * b)^2 + sin(lx * b)^3
          yy = sin(lx * a) + sin(ly * a)^2 + sin(ly * c)^3
        elseif  @mutation == "Beta"
          xx = sin(ly * b) + sin(lx * b)^2
          yy = sin(lx * a) + sin(ly * a)^2
        elseif  @mutation == "Gamma"
          xx = abs(sin(ly * b)) + sin(lx * b)^2
          yy = abs(sin(lx * a)) + sin(ly * a)^2
        elseif  @mutation == "De Jong"
          xx = sin(ly * a) - cos(lx * b)
          yy = sin(lx * c) - cos(ly * d)
        endif
        float xsave = xx, float ysave = yy, lx = xe, ly = ye
        i = i + 1
        if @mutation == "Standard"
          xx = sin(ly * b) + c * sin(lx * b)
          yy = (sin(lx * a) + d * sin(ly * a))
        elseif  @mutation == "Alpha"
          xx = sin(ly * b) + sin(lx * b)^2 + sin(lx * b)^3
          yy = sin(lx * a) + sin(ly * a)^2 + sin(ly * c)^3
        elseif  @mutation == "Beta"
          xx = sin(ly * b) + sin(lx * b)^2
          yy = sin(lx * a) + sin(ly * a)^2
        elseif  @mutation == "Gamma"
          xx = abs(sin(ly * b)) + sin(lx * b)^2
          yy = abs(sin(lx * a)) + sin(ly * a)^2
        elseif  @mutation == "De Jong"
          xx = sin(ly * a) - cos(lx * b)
          yy = sin(lx * c) - cos(ly * d)
        endif
        float dLx = xx - xsave, float dLy = yy - ysave
        float dL2 = dLx * dLx + dLy * dLy
        float df = 1000000000000. * dL2
        float rs = 1/sqrt(df)
        xe = xsave + rs * (xx - xsave)
        ye = ysave + rs * (yy - ysave)
        xx = xsave, yy = ysave
        Lsum = Lsum + log(df)
        L = 0.721347 * lsum / i
        lx = xx
        ly = yy  
      endwhile  
    endwhile
  endif
  int pix[#width, #height]
  if 3 * #width < 4 * #height
    float scale = (#width * #magn) / 4
  else
    scale = (#height * #magn) / 3
  endif
  int x = 0, int y = 0
; initialize array  
  while x < #width
    y = 0
    while y < #height
      pix[x, y] = 0
      y = y + 1
    endwhile
    x = x + 1
  endwhile   
  float xnew = 0
  float ynew = 0
  xx = 0.1
  yy = 0.1
  i = 0
  int iters = #width * #height * @max
  while i < iters
    if @mutation == "Standard"
      xnew = sin(yy * b) + c * sin(xx * b)
      ynew = (sin(xx * a) + d * sin(yy * a))
    elseif  @mutation == "Alpha"
      xnew = sin(yy * b) + sin(xx * b)^2 + sin(xx * b)^3
      ynew = sin(xx * a) + sin(yy * a)^2 + sin(yy * c)^3
    elseif  @mutation == "Beta"
      xnew = sin(yy * b) + sin(xx * b)^2
      ynew = sin(xx * a) + sin(yy * a)^2
    elseif  @mutation == "Gamma"
      xnew = abs(sin(yy * b)) + sin(xx * b)^2
      ynew = abs(sin(xx * a)) + sin(yy * a)^2
    elseif  @mutation == "De Jong"
      xnew = sin(yy * a) - cos(xx * b)
      ynew = sin(xx * c) - cos(yy * d)
    endif
    xx = xnew
    yy = ynew
    float xc = real(#center)
    float yc = imag(#center)
    x = round(((xx - xc) * cos(#angle) - (yy + yc) * \
        sin(#angle)) * scale + #width / 2)
    y = round(((xx - xc) * sin(#angle) + (yy + yc) * \
        cos(#angle)) * scale + #height / 2)
  ; plot the point only if inside array
    if x >= 0 && x < #width && y >= 0 && y < #height
      pix[x, y] = pix[x, y] + 1
      if pix[x, y] > max
        max = pix[x, y]
      endif  
    endif
    i = i + 1
  endwhile
final:
  #index = (log(pix[#x, #y])) / log(max)
default:
  title = "Lat鲻carfian Attractors"
  render = false
  param mutation
    caption = "Mutation"
    enum = "Standard" "Alpha" "Beta" "Gamma" "De Jong"
  endparam
  int param max
    caption = "Sample Density"
    default = 30
    min = 1
  endparam  
  float param a
    caption = "a"
    default = -0.966918
    visible = !@rnd
  endparam  
  param b
    caption = "b"
    default = 2.879879
    visible = !@rnd
  endparam  
  float param c
    caption = "c"
    default = 0.765145
    visible = !@rnd && (@mutation == "Standard" || \
              @mutation == "Alpha" || @mutation == "De Jong")
  endparam  
  float param d
    caption = "d"
    default = 0.744728
    visible = !@rnd && (@mutation == "Standard" || \
              @mutation == "De Jong")
  endparam  
  bool param rnd
    caption = "Random Parameters"
    default = false
  endparam
  float param lthresh
    caption = "Lyapunov"
    default = 0.3
    visible = @rnd
  endparam  
  int param seed  
    caption = "Seed"
    default = 12345678
    visible = @rnd
  endparam  
}

xuefeiyang

3# 榕坚


你们用的这些在哪里有？能给我传一份吗？邮箱422161240@qq.com

榕坚

4# xuefeiyang


你有安装UF程序吧。所有的范例都在里面的。

xuefeiyang

我安装了UF程序，但找不到你所说的那个系列：

下载 (18.25 KB)

2014-7-9 11:59

榕坚

6# xuefeiyang


不是在formula面版中找，在outside的面版中的mt系列中。它是一个分形外部着色程序，象那希尔伯特曲线一样。

榕坚

代码太吓人了：
KleinianLimitSets {
; modified june 2009, original June 2006 Ronald Barnett
;
; Based upon the depth-first tree search algorithm in the book \
; Indra's Pearls by Mumford, Series and Wright. It is considerably faster \
; than most other algorithms in generating the limit set.
; modified Feb 2010 to allow mapping transforms to work.
;
global:
  import "common.ulb"
  import "reb.ulb"

  Circle C = new Circle(0,0)
  Circle cir[4]
  Mobius gen[4]
  MobiusArray word = new MobiusArray(@level)
  complex cs[4]
  float rad[4]
  int tag[]
  setlength(tag,@level)

  float magn = #magn
  float scrsize = #width*#height
  float epsilon = @epsilon

  ; this uses a 640x640 screen with magn = 1 as the reference
  if @screen
    epsilon = epsilon*409600/scrsize/magn
  endif

  int w = 0
  int h = 0
  float Inc = @percentInc*0.01
  float wd = Inc*#width
  float ht = Inc*#height
  int wd2 = round(wd*0.5)
  int ht2 = round(ht*0.5)
  int ilev = 0
  complex ta = @ta
  complex tb = @tb
  complex taj = @taj
  complex tbj = @tbj
  if !@expert && @typeb == "Double Cusps"
    tb = (2,0)
    if @tb2 == "1 100 Cusp"
      ta = (1.99901430461837,-3.94468287279129e-005)
    elseif @tb2 == "1 99 Cusp"
      ta = (1.99899431720688,-4.06536410107979e-005)
    elseif @tb2 == "1 98 Cusp"
      ta = (1.9989737159102,-4.19101798220381e-005)
    elseif @tb2 == "1 97 Cusp"
      ta = (1.99895247534081,-4.32190326141746e-005)
    elseif @tb2 == "1 96 Cusp"
      ta = (1.99893056878569,-4.45829500728911e-005)
    elseif @tb2 == "1 95 Cusp"
      ta = (1.99890796812253,-4.60048582566762e-005)
    elseif @tb2 == "1 94 Cusp"
      ta = (1.99888464372953,-4.74878717791897e-005)
    elseif @tb2 == "1 93 Cusp"
      ta = (1.99886056438849,-4.90353081361803e-005)
    elseif @tb2 == "1 92 Cusp"
      ta = (1.99883569718054, -5.06507032899056e-005 )
    elseif @tb2 == "1 91 Cusp"
      ta = (1.99881000737372, -5.23378286366251e-005 )
    elseif @tb2 == "1 90 Cusp"
      ta = (1.99878345830198, -5.41007094967586e-005 )
    elseif @tb2 == "1 89 Cusp"
      ta = (1.99875601123446, -5.59436452831306e-005 )
    elseif @tb2 == "1 88 Cusp"
      ta = (1.99872762523451, -5.78712315206112e-005 )
    elseif @tb2 == "1 87 Cusp"
      ta = (1.99869825700728, -5.98883839104726e-005 )
    elseif @tb2 == "1 86 Cusp"
      ta = (1.99866786073491, -6.20003646555514e-005 )
    elseif @tb2 == "1 85 Cusp"
      ta = (1.99863638789812, -6.42128112880125e-005 )
    elseif @tb2 == "1 84 Cusp"
      ta = (1.99860378708289, -6.6531768270374e-005 )


……

榕坚

这个应该可以用几何画板来做一下：
ESAFractalTree {
;
; version 1.0
; By Etienne Saint-Amant, 2002/07/02
;
; Go to http://ESAfractal.com for examples and tutorial
;
; Based on FractalTree
; By Samuel Monnier, 11.1.00
; Modified with his kind permission
;
init:
  z = 0
  float arg = 0
  int i = 0
  float j = 0
  float d = 0
  float dist = 1e5
  float x = 0
  float y = 0
  float imagCof = @imaginaryCoefficient
  if @mode == 0
    z = #pixel
  endif

  float dr = pi/@rotincrCoefficient*@rotincr/@altitudeTree
  
loop:
  if @mode == 1
    i = i + 1
    if i == @itertr
      z = #z
    endif
  endif
final:
  ; Performs translation, rotation
  ; and magnification
  z = z - @center
  z = z*exp(-1i*@rot*pi/@rotationchange)
  z = z/@size
  i = 0

  while i < @niter
    i = i + 1
   
    imagCof = imagCof + @alterimaginary

    x = abs(real(z))
    y = imag(z)
    z = x + imagCof * 1i * y
    arg = -atan2(z) + pi/@weirdness            
    ; Estimate distance
    if y >= -1 && y <= 0
      d = x^@power
    elseif y < -1
      d = cabs(z + imagCof * 1i)^@power
    elseif y > 0
      d = cabs(z)^@power
    endif
   
    if d < dist
      dist = d
    endif
   
    ; Computes the rotation to apply
    j = @order-1
    rotincr = j*dr
    while j > 0
      j = j - 2
      if arg < (j+1)*dr
        rotincr = j*dr
      endif
    endwhile
   
    ; Performs transformation
    z = z*exp(imagCof*1i*rotincr)
    z = z*@magnincr
    z = z - imagCof*1i
   
  endwhile
  
  #index = dist
   
default:
  title = "ESA Fractal Tree"
  param mode
    caption = "Mode"
    default = 0
    enum = "Pixel" "Distance Estimator"
  endparam
  param power
    caption = "Dist. Est. Power"
    default = 1.0
  endparam
  param order
    caption = "Tree Order"
    default = 4.0
  endparam
  param rotincr
    caption = "Angle Between Branches"
    default = 40.0
  endparam
  param magnincr
    caption = "Tree Magnification Increment"
    default = 1.9
  endparam
  param center
    caption = "Center"
    default = (0,0)
  endparam
  param rot
    caption = "Rotation"
    default = 0.0
  endparam
  param size
    caption = "Size"
    default = 1.0
  endparam
  param niter
    caption = "Number of Iteration"
    default = 10
  endparam
  param itertr
    caption = "Iteration to observe"
    default = 1
  endparam

  ;parameter added by ESA
  param imaginaryCoefficient
    caption = "Imaginary Coefficient"
    default = 0.75
    hint = "Imaginary Coefficient to change the shape of the branches (default .75)"
  endparam

  ;parameter added by ESA
  param alterimaginary
    caption = "Alter Imaginary Coefficient"
    default = 0.1
    hint = "Alter the Imaginary Coefficient for each iteration (default 0.1)"
  endparam

  ;parameter added by ESA
  param altitudeTree
    caption = "Altitude Tree"
    default = 2.0
    hint = "Create a tree developing more vertically (default 2.0)"
  endparam
  
  ;parameter added by ESA
  param weirdness
    caption = "Weirdness"
    default = 2.2
    hint = "Weirdness (default 2.2)"
   endparam
   
  ;parameter added by ESA
  param rotincrCoefficient
    caption = "Rotation Incrementation"
    default = 120
    hint = "Rotation Incrementation Coefficient (default 120)"
   endparam
   
  ;parameter added by ESA
  param rotationchange
    caption = "Rotation Change"
    default = 50
    hint = "Rotation Change (default 50)"
   endparam
}

xuefeiyang

2# 柳烟


flame!