comment{
;These ere my first try with
;coloring formulas.
;I hope you enjoy these!!!
}

1(BOTH) {
init:
  zold = @seed
  i = 0
loop:
  zold = #z*@Func1(zold) + @c
  i = i + 1
  
final:
  #index = (|zold|+ #maxiter)^@power
  
default:
  title = "1"
  
  Param seed
    Caption = "seed"
    Default = (0.0,4.0)
  Endparam
  
  Param c
    Caption = "c"
    Default = (2.0,0.0)
  Endparam
  
  Param power
    Caption = "power"
    Default = .03
  Endparam
  
  
  Func Func1
    Caption = "Function 1"
    Default = cos()
  Endfunc
  
}

2(BOTH) {
init:
  a = @para
  b = @parb
  float p = @power
  float s = 0
loop:
  If @scelta == 0
    zold = #z + #angle
  Elseif @scelta == 1
    zold = |#z|
  Elseif @scelta == 2
    zold = #angle / #y
  Elseif @scelta == 3
    zold = #x / #z
  Endif
  
  q = a * @Func1 (zold)/@Func2(b * #z)
    
  If @iterazione == 0
    s = 0
  Else
    s = 1
  Endif
  
final:
  #index = s*#maxiter^p + |q| 
default:
  title = "2"
  
  Param para
    caption = "a"
    Default = (1.0,0.0)
  Endparam
  
  Param parb
    caption = "b"
    Default = (1.0,0.0)
  Endparam
  
  Param scelta
    Caption = "Type"
    enum = "z + angle" "|z|" "angle/y" "x/z"
    Default = 1
  Endparam
  
  Param iterazione
    Caption = "Iteration?"
    enum = "no" "yes"
    Default = 0
  Endparam
  
  Param power
    Caption = "power"
    Default = .03
  Endparam
  
  
  Func Func1
    Caption = "Function 1"
    Default = log()
  Endfunc
  
  Func Func2
    Caption = "Function 2"
    Default = asin()
  Endfunc
}

3(BOTH) {
init:
  a = @para
  b = @parb
  c = @parc
  d = @pard
  
  If @scelta == 0
    iterazione = #z
  Elseif @scelta == 1
    iterazione = #pixel
  Elseif @scelta == 2
    iterazione = #angle
  Elseif @scelta == 3
    iterazione = #z + #random*@quanto
  Elseif @scelta == 4
    iterazione = #x
  Elseif @scelta == 5
    iterazione = #y
  Elseif @scelta == 6
    iterazione = |#z|
  Elseif @scelta == 7
    iterazione = #z/#angle
  Endif
  
  
loop:
  z_old = iterazione + #z
  medio = a * @Func1(b * iterazione) + c * @Func2(d * z_old)
  float fine = |@Func3(medio - z_old)|
  
  
final:
  #index = |@Func4(fine)|^@pow
default:
  title = "3"
  
  Param scelta
    Caption = "Type?"
    enum = "z" "pixel" "angle" "z + random" "real" "imag" "|z|" "z / angle"
    Default = 0
  Endparam
  
  Param para
    Caption = "a"
    Default = (1.0,0.0)
  Endparam
  
  Param parb
    Caption = "b"
    Default = (1.0,0.0)
  Endparam
  
  Param parc
    Caption = "c"
    Default = (1.0,0.0)
  Endparam
  
  Param pard
    Caption = "d"
    Default = (1.0,0.0)
  Endparam
  
  Param pow
    Caption = "power"
    Default = 0.5
  Endparam
  
  Param quanto
    Caption = "How many rand?"
    Hint = "This works only with the choice of random"
    Default = .52
  Endparam
  
  Func Func1
    Caption = "Function 1"
    Default = asin()
  Endfunc
  
  Func Func2
    Caption = "Function 2"
    Default = round()
  Endfunc
  
  Func Func3
    Caption = "Function 3"
    Default = atanh()
  Endfunc
  
  Func Func4
    Caption = "Function 4"
    Default = sinh()
  Endfunc
  
    
}


4(BOTH) {
init:
  a = @para
  b = @parb
  c = @parc
  d = @pard
  z_old = 0
  
  If @scelta == 0
    iterazione = #z
  Elseif @scelta == 1
    iterazione = #pixel
  Elseif @scelta == 2
    iterazione = #angle
  Elseif @scelta == 3
    iterazione = #z + #random*@quanto
  Elseif @scelta == 4
    iterazione = #x
  Elseif @scelta == 5
    iterazione = #y
  Elseif @scelta == 6
    iterazione = |#z|
  Elseif @scelta == 7
    iterazione = #z/#angle
  Endif
  
  
loop:
  z_old = iterazione + #z
  medio = a * @Func1(b * iterazione) + c * @Func2(d * z_old)
  fine = @Func3(medio - z_old)
  
  
final:
  
  magma = z_old * medio
  
  fine = fine/magma
  
  
  
  
  #index = |@Func4(fine)|^@pow
default:
  title = "4"
  
  Param scelta
    Caption = "Type?"

    enum = "z" "pixel" "angle" "z + random" "real" "imag" "|z|" "z / angle"
    Default = 0
  Endparam
  
  Param para
    Caption = "a"
    Default = (1.0,0.0)
  Endparam
  
  Param parb
    Caption = "b"
    Default = (1.0,0.0)
  Endparam
  
  Param parc
    Caption = "c"
    Default = (1.0,0.0)
  Endparam
  
  Param pard
    Caption = "d"
    Default = (1.0,0.0)
  Endparam
  
  Param pow
    Caption = "power"
    Default = 0.5
  Endparam
  
  Param quanto
    Caption = "How many rand?"
    Hint = "This works only with the choice of random"
    Default = .52
  Endparam
  
  Func Func1
    Caption = "Function 1"
    Default = cos()
  Endfunc
  
  Func Func2
    Caption = "Function 2"
    Default = exp()
  Endfunc
  
  Func Func3
    Caption = "Function 3"
    Default = atanh()
  Endfunc
  
  Func Func4
    Caption = "Function 4"
    Default = sinh()
  Endfunc
  
    
}

FunctionsTrap(BOTH) {

;Daniele Lupo
;alcamese@libero.it
;
;A very special thank to Mark Townsend
;for his basic explanation about traps.
;
;This formula create traps based on shapes
;of many function, as parabola, lines etc.
;you can define the function, how many function,
;his position, rotation and others.
;
;Useful with sqrt or cuberoot transfer functions, 
;but others work well, too.

;Sorry, but it's a little slow with many functions. :-(



init:

;Basic variables of the formula. With these variables
;you can describe the formula 
  
  float closest = 1e15   ; starting of calculation to min dist
  float farthest = 0     ; starting of calculation of max dist
  float a = @para        ; a parameter of functions
  float b = @parb        ; b parameter of functions
  float c = @parc        ; c parameter of functions 
  float d = @pard        ; Try to discover his function!!! :-)
  float f = @parf        ; This too!!!
  float g = @parg        ; The last parameter!!!
  float rotation = @rot  ; rotation angle of functions
  float distance = 1e15     ; distance from z to the shape
  float temp_distance = 0
  float tolerance = @tol ; tolerance to numeric algorithm to find
                         ; the minimun distance
  float q = 0            ; the variable to plot the function.
                         ; I prefer use it, instead of x,
                         ; to solve problems with the builtin #x
  
  complex origin = @center_of_axis    ; translation of origin of
                                      ; axis

;These variables describe the ripetition of formula pattern:
;number of function and the position in the complex plane.
    
  int functions = @num_func      ; the number of functions that create the trap
  float x_translation = @x_tran  ; x translation of every function in the formula
  float y_translation = @y_tran  ; y translation of every function in the formula  
  float group_rotation = @g_rot  ; rotations of other functions added to the original
                                 ; rotation of original function
;others variables
  
  float qold = 0         ; newton algorithm use it
  float diff = 1e10      ; |q - qold| in newton algorithm
  float z_x = 0.0        ; x value of z
  float z_y = 0.0        ; y value of z
  float fq = 0.0         ; value of f(q)
  float der1fq = 0.0     ; value of the first derivate of f(q)
  float der2fq = 0.0     ; value of the second derivate of f(q)
  int max_ite =  @maxite ; max iteration of newton algorithms
  float ang = 0          ; I use it to rotate the function
  float length = 0
  int cont = 0
  int iter_closest = 0       
  int iter_farthest = 0
  float sum = 0
  float difference = 0
  float multiplication = 1
  int num_ite = 0
  
  
loop:
  
;initialization of q

  If (@type_q == 0)
    q = real(#z)
  Elseif (@type_q == 1)
    q = imag(#z)
  Elseif (@type_q == 2)
    q = real(#z)+imag(#z)
  Elseif (@type_q == 3)
    q = real(#z)-imag(#z)
  Elseif (@type_q == 4)
    q = imag(#z)-real(#z)
  Elseif (@type_q == 5)
    q = real(#z)*imag(#z)
  Elseif (@type_q == 6)
    q = real(#z)/imag(#z)
  Elseif (@type_q == 7)
    q = imag(#z)/real(#z)
  Elseif (@type_q == 8)
    q = real(#z)^imag(#z)
  Elseif (@type_q == 9)
    q = imag(#z)^real(#z)
  Elseif (@type_q == 10)
    q = real(@init_func(real(#z)))
  Elseif (@type_q == 11)
    q = real(@init_func(imag(#z)))
  Elseif (@type_q == 12)
    q = real(@init_func(real(#z)+imag(#z)))
  Elseif (@type_q == 13)
    q = real(@init_func(real(#z)-imag(#z)))
  Elseif (@type_q == 14)
    q = real(@init_func(imag(#z)-real(#z)))
  Elseif (@type_q == 15)
    q = real(@init_func(real(#z)*imag(#z)))
  Elseif (@type_q == 16)
    q = real(@init_func(real(#z)/imag(#z)))
  Elseif (@type_q == 17)
    q = real(@init_func(imag(#z)/real(#z)))
  Elseif (@type_q == 18)
    q = real(@init_func(real(#z)^imag(#z)))
  Elseif (@type_q == 19)
    q = real(@init_func(imag(#z)^real(#z)))
  Elseif (@type_q == 20)
    q = real(#pixel)
  Elseif (@type_q == 21)
    q = imag(#pixel)
  Elseif (@type_q == 22)
    q = real(#pixel)+imag(#pixel)
  Elseif (@type_q == 23)
    q = real(#pixel)-imag(#pixel)
  Elseif (@type_q == 24)
    q = imag(#pixel)-real(#pixel)
  Elseif (@type_q == 25)
    q = real(#pixel)*imag(#pixel)
  Elseif (@type_q == 26)
    q = real(#pixel)/imag(#pixel)
  Elseif (@type_q == 27)
    q = imag(#pixel)/real(#pixel)
  Elseif (@type_q == 28)
    q = real(#pixel)^imag(#pixel)
  Elseif (@type_q == 29)
    q = imag(#pixel)^real(#pixel)
  Elseif (@type_q == 30)
    q = @seed
  Endif

; I use tranformation values, like rotation etc.
; in z calculation, instead of in the function.
; This simplifies calculus  
  
 
  z_x = real(#z)-real(origin)
  z_y = imag(#z)-imag(origin)
  ang = atan2(z_x+z_y*1i)
  length = cabs(z_x+z_y*1i)
  z_x = length * cos(ang + rotation*#pi/180)
  z_y = length * sin(ang + rotation*#pi/180)
  
  
;Now, I must calculate the distance from every function

  cont = 0

  While (cont < functions)
    cont = cont+1
    diff = 1e10

; calculation of q point that gives the minimun distance
; formula obtained derivating the distance from the function formula.
;  dist = (x-q)^2 + (y-f(q)^2, and appliyng newton
; formula to this. When derivate is equal to 0, then I've found
; the q value that minimize the distance.
; Then I calculate the distance from z to the point (q, f(q))  
  
    num_ite = max_ite
    
    While (diff > tolerance && num_ite > 0)
      
      num_ite = num_ite-1
      qold = q

;The type of function    
    
      If (@type_f == 0)
        fq = a*q+b
        der1fq = a
        der2fq = 0
      Elseif (@type_f == 1)
        fq = a*q^2+b*q+c
        der1fq = 2*a*q+b
        der2fq = 2*a
      Elseif (@type_f == 2)
        fq = a*q^3+b*q^2+c*q+d
        der1fq = 3*a*q^2+2*b*q+c
        der2fq = 6*a*q+2*b
      Elseif (@type_f == 3)
        fq = a*q^b+c*q^d
        der1fq = a*b*q^(b-1)+c*d*q^(d-1)
        der2fq = a*b*(b-1)*q^(b-2)+c*d*(d-1)*q^(d-2)
      Elseif (@type_f == 4)
        fq = a*sin(b*q+c)+d*cos(f*q+g)
        der1fq = a*b*cos(c+b*q)-d*f*sin(g+f*q)
        der2fq = -d*f*f*cos(g+f*q)-a*b*b*sin(c+b*q)
      Elseif (@type_f == 5)
        fq = a*q*sin(b*q+c)
        der1fq = a*b*q*cos(b*q+c)+a*sin(b*q+c)
        der2fq = 2*a*b*cos(b*q+c)-a*b*b*q*sin(b*q+c)
      Elseif (@type_f == 6)
        fq = a*sin(b*q+c)/q
        der1fq = a*(b*q*cos(c+b*q)-sin(c+b*q))/q^2
        der2fq = -a*(2*b*q*cos(c+b*q)+(-2+b^2*q^2)*sin(c+b*q))/q^3
      Elseif (@type_f == 7)
        fq = a*q^2+b*q+d*sin(f*q+g)
        der1fq = 2*a*q+b+d*f*cos(f*q+g)
        der2fq = 2*a-d*f*f*sin(f*q+g)
      Elseif (@type_f == 8)
        fq = a*tan(b*q+c)*exp(d*q+f)
        der1fq = a*exp(f+d*q)*(2*b+d*sin(2*c+2*b*q))/(1+cos(2*c+2*b*q))
        der2fq = 0.25*a*exp(f+d*q)*(8*b*d*cos(c+b*q)+2*(4*b^2+d^2+d^2*cos(2*(c+b*q)))*sin(c+b*q))/(cos(c+b*q)^3)
      Elseif (@type_f == 9)
        fq = a*cosh(b*q+c)+d*sinh(f*q+g)
        der1fq = d*f*cosh(g+f*q)+a*b*sinh(c+b*q)
        der2fq = a*b*b*cosh(c+b*q)+d*f*f*sinh(g+f*q)
      Endif

;Newton formula
    
      q = q - (((z_y-fq+(q-z_x)*der1fq)^2-(z_y-fq)*((q-z_x)^2+z_y^2+fq^2-2*z_y*fq)*der2fq)/(((q-z_x)^2+z_y^2+fq^2-2*z_y*fq)^(3/2)))
      diff = |q-qold|
  
    Endwhile   ;newton loop
    
    
  temp_distance = real(sqr(z_x-q)+sqr(z_y-fq))
    
;calculation of variables for the new function, with his trantaltion and rotazion

    z_x = z_x-x_translation
    z_y = z_y-y_translation
    ang = atan2(z_x+z_y*1i)
    length = cabs(z_x+z_y*1i)
    z_x = length * cos(ang + group_rotation*#pi/180)
    z_y = length * sin(ang + group_rotation*#pi/180)
    
    
    
    If (distance > temp_distance)
      distance = temp_distance
    Endif
    
  Endwhile    ;number of functions loop
  

;Now I've the distance. Now I verify the trap


  If (distance < closest)
    closest = distance
    iter_closest = iter_closest+1
  Endif
  
  If (distance > farthest)
    farthest = distance
    iter_farthest = iter_farthest+1
  Endif
  
    
  sum = sum + distance
  
  difference = difference - distance
  
  multiplication = multiplication * distance
            
final:

  If (@type_trap == 0)
    #index = sqrt(closest)
  Elseif (@type_trap == 1)
    #index = sqrt(farthest)
  Elseif (@type_trap == 2)
    #index = iter_closest*.1
  Elseif (@type_trap == 3)
    #index = iter_farthest*1
  Elseif (@type_trap == 4)
    #index = sum
  Elseif (@type_trap == 5)
    #index = sum/#numiter
  Elseif (@type_trap == 6)
    #index = difference
  Elseif (@type_trap == 7)
    #index = multiplication
  Endif
  
default:
  title = "Functions Trap"
  
  Param type_q
    Caption = "seed"
    Enum = "x" "y" "x+y" "x-y" "y-x" "x*y" "x/y" "y/x" "x^y" "y^x" "Q(x)" "Q(y)" \
           "Q(x+y)" "Q(x-y)" "Q(y-x)" "Q(x*y)" "Q(x/y)" "Q(y/x)" "Q(x^y)" "Q(y^x)" \
           "pix_x" "pix_y" "pix_x+pix_y" "pix_x-pix_y" "pix_y-pix_x" "pix_x*pix_y" \
           "pix_x/pix_y" "pix_y/pix_x" "pix_x^pix_y" "pix_y^pix_x" "manual"
    Default = 0
    Hint = "This is the initial point to start the algorithm to find the minimum \
            distance. Cause of this, many points can give the same solution, and work \
            in the same way, depend of the formula."
  Endparam
  
  Param type_f
    Caption = "Function"
    Enum = "ax+b" "ax^2+bx+c" "ax^3+bx^2+cx+d" "ax^b+cx^d" "a*sin(bx+c)+d*cos(fx+g)" \
           "ax*sin(bx+c)" "asin(bx+c)/x)" "ax^2+bx+c+d*sin(fx+g)" "a*tan(b*cx)*exp(dx+f)" \
           "a*cosh(bx+c)+d*sinh(fx+g)"
    Hint = "There are all the functions used as trap. Try, try and try!"
    Default = 0
  Endparam
  
  Param type_trap
    Caption = "Trap"
    Enum = "Closest" "Farthest" "iter cl." "iter far." "sum" "average" "difference" "multiplication"
    Default = 0
    Hint = "This specifies the type of trap of function."
  Endparam
  
  Param center_of_axis
    Caption = "Center"
    Default = (0.0,0.0)
    Hint = "This is the centre of xy axes of function"
  Endparam
  
  Param rot
    Caption = "Rotation"
    Default = 0.0
    Hint = "This is the rotation of xy axes respect to imaginary plane."
  Endparam
  
  Param num_func
    Caption = "N° of funcs"
    Default = 1
    Min = 1
    Hint = "This number is number of functions that must be used as trap. \
            You must use traslations values, or you can't use others, because \
            They're one on top of the other."
  Endparam
  
  Param x_tran
    Caption = "x translation"
    Default = 0.5
    Hint = "This is the x translation of every function from the previosus. It doesn't \
            work with only one function"
  Endparam
  
  Param y_tran
    Caption = "y translation"
    Default = 0.5
    Hint = "This is the y translation of every function from the previosus. It doesn't \
            work with only one function"
  Endparam
  
  Param g_rot
    Caption = "Group rot."
    Default = 0.0
    Hint = "This is the angle of every function respect to the previous plotted function"
  Endparam
  
  Param maxite
    Caption = "Max. Iter."
    Default = 30
    Hint = "This is the max number of iteration of newton algorithm. It's useful \
            when the the formnula can0t find a minimun, and exit from a probable \
            infinite loop."
    Min = 1
  Endparam
  
  Param para
    Caption = "a"
    Default = 1.0
    Hint = "a parameter of functions"
  Endparam
  
  Param parb
    Caption = "b"
    Default = 0.0
    Hint = "b parameter of funcions"
  Endparam
  
  Param parc
    Caption = "c"
    Default = 1.0
    Hint = "c parameter of functions"
  Endparam
  
  Param pard
    Caption = "d"
    Default = 1.0
    Hint = "d parameter of functions"
  Endparam
  
  Param parf
    Caption = "f"
    Default = 1.0
    Hint = "f parameter of funcions"
  Endparam
  
  Param parg
    Caption = "g"
    Default = 1.0
    Hint = "g parameter of functions"
  Endparam
  
  Param tol
    Caption = "Tolerance"
    Default = 1e-2
    Hint = "This is the tolerance of numerical newton algorithm"
  Endparam
    
  Param seed
    Caption = "manual"
    Default = 0.0
    Hint = "This works only if it's selected in seed selection; with this, you can \
            manually add the seed of the newton iteration."
  Endparam
  
  Func init_func
    Caption = "Q()"
    Hint = "This is the function that appear in options of seed. If this does not \
            appear in the seed, doesn't work"
    Default = cos()
  Endfunc
  
    
      
  
  
}



dan-ucl-Waves(BOTH) {
;
; Polar trap based on polar functions.
; Functions are customizable in many aspects.
; You can choose different types of polar
; functions, but they are very flexible.
; you can change aspect ration and center
; of this trap.
; there's commented code. It correspond
; to another trap that I'm working on. It's
; not avaiable at the monent, because
; it doesn not work like I want.
; I'd like to improve this in future, but for
; now, this is the result! :-)
; 
; Any suggestion for improvement of this formula
; is accepted, of course! Let me know what and how
; I can add something.
;
; 18 Sep 2003
;
; Morph Parameters added by Toby Marshall
; Thanks Toby. These parameters are wonderful!
; 
; Created by Daniele Lupo
;
; danwolf@ciaoweb.it


global:

; From parameters to normal names 
  float ray_of_shape = @r
  float amplitude = @am
  float frequency = @fr
  int type_trap = @t_t
  float iterations_control = @i_c
  float thickness = 1/@th
  float aspect_ratio = @a_r
  int trapped_value_type = @t_v_t
  complex _center = @cn
  float scale = @sc
  bool ray_distance = @r_d
  
  If ray_distance == true
    bool solid_type = @s_t
  Else
    bool solid_type = false
  Endif
  
;  float slice_angle = 2*#pi/@s_a
;  float n_of_slices = @s_a
  int wave_function = @w_f
  float amplitude_1_1 = @am_1_1
  float amplitude_1_2 = @am_1_2
  float frequency_1_1 = @fr_1_1
  float frequency_1_2 = @fr_1_2
  float phase_1_1 = @ph_1_1*#pi/180
  float amplitude_2_1 = @am_2_1
  float amplitude_2_2 = @am_2_2
  float frequency_2_1 = @fr_2_1
  float frequency_2_2 = @fr_2_2
  float phase_2_1 = @ph_2_1*#pi/180
  float amplitude_3_1 = @am_3_1
  float amplitude_3_2 = @am_3_2
  float frequency_3_1 = @fr_3_1
  float phase_3_1 = @ph_3_1*#pi/180
  float rotation = @rot*#pi/180
  float amplitude_4_1 = @am_4_1
  float amplitude_4_2 = @am_4_2
  float frequency_4_1 = @fr_4_1
  float frequency_4_2 = @fr_4_2
  float phase_4_1 = @ph_4_1*#pi/180
;  float ray_1 = @r_1
;  float ray_2 = @r_2


init:
  
  float minimum = 1e20
  float maximum = 0.0
  float sum = 0.0
  float product = 1.0
  float first = 0.0
  float last = 0.0
  float _angle = 0.0
  complex real_imag = (0.0,0.0)
  int control = 0
  int iterations = 0
  float t = 0.0
  float error_value = 0.0  ;This contains error_value that it's calculated in trap
  complex trapped_value = (0.0,0.0) ;this contains changes of the z value
  complex ray = 0
  int iter = 0
;  float slice = 0.0
;  float temp_angle = 0.0
;  int angle_iter = 0
;  float distance_from_border = 0.0
  

loop:
  complex z1 = #z
  if @change == "1"
    z1 = @nufunc(z1-@tw)^@exp1m
  elseif @change == "2"
    z1 = @nufunc(z1-@tw)+@nufunc2((z1-@tw2)^@exp1m)
  elseif @change == "3"
    z1 = @nufunc(z1-@tw)-@nufunc2((z1-@tw2)^@exp1m)
  elseif @change == "4"
    z1 = @nufunc(z1-@tw)*@nufunc2((z1-@tw2)^@exp1m)
  elseif @change == "5"
    z1 = @nufunc(z1-@tw)/@nufunc2((z1-@tw2)^@exp1m)        
  elseif @change == "6"
    z1 = @nufunc(z1-@tw)^@nufunc2((z1-@tw2)^@exp1m)
  elseif @change == "7"
   if @opm == "+"
     z1 = @nufunc(z1-@tw)+(@nufunc2((z1-@tw2)^@exp1m)+@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "-"
     z1 = @nufunc(z1-@tw)+(@nufunc2((z1-@tw2)^@exp1m)-@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "*"
     z1 = @nufunc(z1-@tw)+(@nufunc2((z1-@tw2)^@exp1m)*@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "/"
     z1 = @nufunc(z1-@tw)+(@nufunc2((z1-@tw2)^@exp1m)/@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "^"
     z1 = @nufunc(z1-@tw)+(@nufunc2((z1-@tw2)^@exp1m)^@nufunc3((z1-@tw3)^@exp2m))
   endif
  elseif @change == "8"
   if @opm == "+"
     z1 = @nufunc(z1-@tw)-(@nufunc2((z1-@tw2)^@exp1m)+@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "-"
     z1 = @nufunc(z1-@tw)-(@nufunc2((z1-@tw2)^@exp1m)-@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "*"
     z1 = @nufunc(z1-@tw)-(@nufunc2((z1-@tw2)^@exp1m)*@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "/"
     z1 = @nufunc(z1-@tw)-(@nufunc2((z1-@tw2)^@exp1m)/@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "^"
     z1 = @nufunc(z1-@tw)-(@nufunc2((z1-@tw2)^@exp1m)^@nufunc3((z1-@tw3)^@exp2m))
   endif
  elseif @change == "9"
   if @opm == "+"
     z1 = @nufunc(z1-@tw)*(@nufunc2((z1-@tw2)^@exp1m)+@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "-"
     z1 = @nufunc(z1-@tw)*(@nufunc2((z1-@tw2)^@exp1m)-@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "*"
     z1 = @nufunc(z1-@tw)*(@nufunc2((z1-@tw2)^@exp1m)*@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "/"
     z1 = @nufunc(z1-@tw)*(@nufunc2((z1-@tw2)^@exp1m)/@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "^"
     z1 = @nufunc(z1-@tw)*(@nufunc2((z1-@tw2)^@exp1m)^@nufunc3((z1-@tw3)^@exp2m))
   endif
  elseif @change == "10"
   if @opm == "+"
     z1 = @nufunc(z1-@tw)/(@nufunc2((z1-@tw2)^@exp1m)+@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "-"
     z1 = @nufunc(z1-@tw)/(@nufunc2((z1-@tw2)^@exp1m)-@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "*"
     z1 = @nufunc(z1-@tw)/(@nufunc2((z1-@tw2)^@exp1m)*@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "/"
     z1 = @nufunc(z1-@tw)/(@nufunc2((z1-@tw2)^@exp1m)/@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "^"
     z1 = @nufunc(z1-@tw)/(@nufunc2((z1-@tw2)^@exp1m)^@nufunc3((z1-@tw3)^@exp2m))
   endif
  elseif @change == "11"
   if @opm == "+"
     z1 = @nufunc(z1-@tw)^(@nufunc2((z1-@tw2)^@exp1m)+@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "-"
     z1 = @nufunc(z1-@tw)^(@nufunc2((z1-@tw2)^@exp1m)-@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "*"
     z1 = @nufunc(z1-@tw)^(@nufunc2((z1-@tw2)^@exp1m)*@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "/"
     z1 = @nufunc(z1-@tw)^(@nufunc2((z1-@tw2)^@exp1m)/@nufunc3((z1-@tw3)^@exp2m))
   elseif @opm == "^"
     z1 = @nufunc(z1-@tw)^(@nufunc2((z1-@tw2)^@exp1m)^@nufunc3((z1-@tw3)^@exp2m))
   endif   
  endif
  


iter = iter + 1

; These code is used to change z properties

  If trapped_value_type == 0      ;z
    trapped_value = real(aspect_ratio*z1)+flip(imag(z1))
  Elseif trapped_value_type == 1  ;|z|
    trapped_value = |real(aspect_ratio*z1)+flip(imag(z1))|
  Elseif trapped_value_type == 2  ;real(z)
    trapped_value = real(aspect_ratio*z1)
  Elseif trapped_value_type == 3  ;imag(z)
    trapped_value = imag(z1)
  Elseif trapped_value_type == 4  ;angle(z)
    trapped_value = atan2(real(aspect_ratio*z1)+flip(imag(z1)))
  Elseif trapped_value_type == 5  ;f1(f2(z))
    trapped_value = @value_func_1(@value_func_2(real(aspect_ratio*z1)+flip(imag(z1))))
  Elseif trapped_value_type == 6  ;f1(real)+f2(imag)
    trapped_value = @value_func_1(aspect_ratio*real(z1))+@value_func_2(imag(z1))
  Elseif trapped_value_type == 7  ;f1(real)-f2(imag)
    trapped_value = @value_func_1(aspect_ratio*real(z1))-@value_func_2(imag(z1))
  Elseif trapped_value_type == 8  ;f1(imag)-f2(real)
    trapped_value = @value_func_1(imag(z1))-@value_func_2(aspect_ratio*real(z1))
  Elseif trapped_value_type == 9  ;f2(imag)*f1(real)
    trapped_value = @value_func_1(aspect_ratio*real(z1))*@value_func_2(imag(z1))
  Elseif trapped_value_type == 10  ;f1(real)/f2(imag)
    trapped_value = @value_func_1(aspect_ratio*real(z1))/@value_func_2(imag(z1))
  Elseif trapped_value_type == 11  ;f1(imag)/f2(real)
    trapped_value = @value_func_1(imag(z1))/@value_func_2(real(z1)*aspect_ratio)
  Elseif trapped_value_type == 12  ;f1(real)^f2(imag)
    trapped_value = @value_func_1(aspect_ratio*real(z1))^@value_func_2(imag(z1))
  Elseif trapped_value_type == 13  ;f1(imag)^f2(real)
    trapped_value = @value_func_1(imag(z1))^@value_func_2(aspect_ratio*real(z1))
  Endif

;Now we change center and rotation of our shape

  trapped_value = trapped_value - _center
  trapped_value = cabs(trapped_value)*exp(flip(rotation+atan2(trapped_value)))


; These lines determine amplitude of the
; polar function with the z angle

  t = atan2(trapped_value)
  If wave_function == 0
    ray = ray_of_shape+amplitude*@wave(frequency*t)
  Elseif wave_function == 1
    ray = amplitude_1_1*@wave_1_1(frequency_1_1*t)+amplitude_1_2*@wave_1_2(frequency_1_2*t+phase_1_1)
  Elseif wave_function == 2
    ray = amplitude_2_1*@wave_2_1(frequency_2_1*t)*amplitude_2_2*@wave_2_2(frequency_2_2*t+phase_2_1)
  Elseif wave_function == 3
    ray = amplitude_3_1*@wave_3_1(amplitude_3_2*@wave_3_2(frequency_3_1*t+phase_3_1))
  Elseif wave_function == 4
    ray = amplitude_4_1*@wave_4_1(frequency_4_1*t)*amplitude_4_2*@wave_4_2(frequency_4_2*t+phase_4_1)
;  Elseif wave_function == 5
;    temp_angle = atan2(trapped_value - _center)
;    If (temp_angle < 0)
;        temp_angle = temp_angle + #pi * 2
;      Endif

;    While temp_angle > slice_angle
;      temp_angle = temp_angle - slice_angle
;     angle_iter = angle_iter + 1
;    Endwhile
;    If ((angle_iter%2) == 0)
;      ray = ray_1
;    Else
;      ray = ray_2
;    Endif
;    
;    distance_from_border = atan2 (slice_angle - temp_angle)
;    If distance_from_border > (-atan2(temp_angle-slice_angle))
;      distance_from_border = -atan2(temp_angle-slice_angle)
;    Endif
;    If (cabs(ray_1-ray_2) > distance_from_border)
;      ray = distance_from_border
;    Endif
    
  Endif
  

  actual = cabs(trapped_value)*scale
  
  
; Use ray Option for first trap shape
  
  If wave_function == 0
    If ray_distance == true
      error_value = cabs(actual - ray) * thickness
    Else
      error_value = cabs(ray)*thickness*cabs(trapped_value)
    Endif
  Else
    error_value = cabs(actual - ray) * thickness
  Endif
  
  
; Type of Trap

  If ((type_trap == 0) || (type_trap == 9) || (type_trap == 8) || (type_trap == 10) || (type_trap == 7))      ;minimum & iterations & real min & max min & angle
    If error_value < minimum
      minimum = real(@trap_func(error_value))
      real_imag = @trap_func(trapped_value)
      iterations = iter
      _angle = real(@trap_func(atan2(trapped_value)))
      If (_angle < 0)
        _angle = real(@trap_func(_angle + #pi * 2))
      Endif
    Endif
  Elseif type_trap == 1  ;maximum
    If error_value > maximum
      maximum = real(@trap_func(error_value))
    Endif
  Elseif ((type_trap == 2) || (type_trap == 3))  ;maximum & average
    sum = sum + real(@trap_func(error_value))
  Elseif type_trap == 4  ;product
    product = product * real(@trap_func(error_value))
  Elseif type_trap == 5 ;first
    If control == 0
      first = real(@trap_func(error_value))
      control = 1
    Endif
  Elseif type_trap == 6 ;last
    last = real(@trap_func(error_value))
  Endif
  
    
final:

;Case of average
  If type_trap == 3
    sum = sum/#numiter
  Endif
  
  If solid_type == true
    If minimum*thickness > 1
      #solid = true
    Endif
  Endif
  
    
  If type_trap == 0      ;minimum
    #index = minimum * thickness
  Elseif type_trap == 1  ;maximum
    #index = maximum * thickness
  Elseif ((type_trap == 2) || (type_trap == 3))  ;sum & average
    #index = sum
  Elseif type_trap == 4  ;product
    #index = product
  Elseif type_trap == 5  ;first
    #index = first * thickness
  Elseif type_trap == 6  ;last
    #index = last * thickness
  Elseif type_trap == 7 ;iterations
    #index = iterations/iterations_control
  Elseif type_trap == 8  ;real minimum
    #index = abs(real(real_imag))
  Elseif type_trap == 9  ;imag minimum
    #index = abs(imag(real_imag))
  Elseif type_trap == 10  ;angle
    #index = _angle
  Endif
  
  
default:
  title = "Waves"
  
  Param w_f
    Caption = "Wave Funcion"
    Default = 0
    Enum = "Radial Func" "Sum of 2 Func" "Product of 2 Func" "Func of Func" \
           "Division of Two Func"
    Hint = "This parameter choose the type of polar function that you want \
            to use. Parameters are all editable"
  Endparam  
  
  Param t_t
    Caption = "Type of Trap"
    Default = 0
    Enum = "Minimum" "Maximum" "Average" "Sum" "Product" \
           "First value" "Last value" "Iterations" "Real" \
           "Imaginary" "Angle"
    Hint = "These are the kinds of trap of this formula"
  Endparam
  
  Func trap_func
    Caption = "Funciton Trap"
    Default = ident()
    Hint = "This function is applied to the trap value at every iteration."
  Endfunc
  
  Param t_v_t
    Caption = "Trapped Value"
    Default = 0
    Enum = "z" "|z|" "Real" "Imag" "Angle" "f1(f2(z))" "f1(Real)+f2(Imag)" \
           "f1(Real)-f2(Imag)" "f1(Imag)-f2(Real)" "f1(Real)*f2(Imag)" "f1(Real)/f2(Imag)" "f1(Imag)/f2(Real)" \
           "f1(Real)^f2(Imag)" "f1(Imag)^f2(Real)"
    Hint = "This is the value that it's trapped in interations."
  Endparam
  
  Func value_func_1
    Caption = "f1"
    Default = ident()
    Visible = (@t_v_t >= 5)
    Hint = "First function selected for the Trapped Value"
  Endfunc
  
  Func value_func_2
    Caption = "f2"
    Default = ident()
    Hint = "Second function selected for the Trapped Value"
    Visible = (@t_v_t >= 5)
  Endfunc
  
  Heading
    Caption = "---Parameters for Radial Func---"
    Visible = (@w_f == 0)
  Endheading
  
  Func wave
    Caption = "Wave Function"
    Default = cos()
    Hint = "This is the function that you use as wave function."
    Visible = (@w_f == 0)
  Endfunc
  
  Param r
    Caption = "Ray"
    Default = 1.3

    Min = 0.0
    Hint = "This is the ray of the function, without considering amplitude."
    Visible = (@w_f == 0)
  Endparam
  
  Param am
    Caption = "Amplitude"
    Default = 0.2
    Min = 0.0
    Hint = "This is the amplitude of waves of the function"
    Visible = (@w_f == 0)
  Endparam
  
  Param fr
    Caption = "Frequency"
    Default = 10.0
    Hint = "This is frequency of waves of the function."
    Visible = (@w_f == 0)
  Endparam
  
  Param r_d
    Caption = "Use shape ray?"
    Default = True
    Hint = "Click on in tif you want to create a shape with a determined \
            ray. Don't click to have a central shape trap."
    Visible = (@w_f == 0)
  Endparam
  


  Heading
    Caption = "-------------------------------------"
    Visible = (@w_f == 0)
  Endheading
  
  Heading
    Caption = "--- Parameters for Sum of 2 Func ---"


    Visible = (@w_f == 1)

  Endheading
  
  Param am_1_1
    Caption = "Amplitude 1"
    Default = 1.0
    Hint = "First amplitude value of this polar function."
    Visible = (@w_f == 1)
  Endparam
  
  Param fr_1_1
    Caption = "Frequency 1"
    Default = 1.0
    Hint = "This is the first frequency value of this polar function."
    Visible = (@w_f == 1)
  Endparam
  
  Param am_1_2
    Caption = "Amplitude 2"
    Default = 0.2
    Hint = "This is the second amplitude value of this polar function."
    Visible = (@w_f == 1)
  Endparam
  
  Param fr_1_2
    Caption = "Frequency 2"
    Default = 5
    Hint = "I'd like to know if you can find the meaning of this parameter :-)"
    Visible = (@w_f == 1)
  Endparam
  
  Param ph_1_1
    Caption = "Phase"
    Default = 0.0
    Hint = "This parameter specifies phase between two function of this shape."
    Visible = (@w_f == 1)
  Endparam
  
  Func wave_1_1
    Caption = "Wave Function 1"
    Default = sin()
    Hint = "This is the first function of this wave trap"
    Visible = (@w_f == 1)
  Endfunc
  
  Func wave_1_2
    Caption = "Wave Function 2"
    Default = sin()
    Hint = "This is the second function of this polar function."
    Visible = (@w_f == 1)
  Endfunc
  
  Heading
    Caption = "----------------------------------------"
    Visible = (@w_f == 1)
  Endheading
  
  Heading
    Caption = "--- Parameter for Product of 2 Func ---"
    Visible = (@w_f == 2)
  Endheading
  
  Param am_2_1
    Caption = "Amplitude 1"
    Default = 1.0
    Hint = "First amplitude value of this polar function."
    Visible = (@w_f == 2)
  Endparam
  
  Param fr_2_1
    Caption = "Frequency 1"
    Default = 4.0
    Hint = "This is the first frequency value of this polar function."
    Visible = (@w_f == 2)
  Endparam
  

  Param am_2_2
    Caption = "Amplitude 2"
    Default = 1.0
    Hint = "This is the second amplitude value of this polar function."
    Visible = (@w_f == 2)
  Endparam
  
  Param fr_2_2
    Caption = "Frequency 2"
    Default = 10.0
    Hint = "I'd like to know if you can find the meaning of this parameter :-)"
    Visible = (@w_f == 2)
  Endparam
  
  Param ph_2_1
    Caption = "Phase"
    Default = 0.0
    Hint = "This parameter specifies phase between two function of this shape."
    Visible = (@w_f == 2)
  Endparam
  
  Func wave_2_1
    Caption = "Wave Function 1"
    Default = sin()
    Hint = "This is the first function of this wave trap"
    Visible = (@w_f == 2)
  Endfunc


  
  Func wave_2_2
    Caption = "Wave Function 2"
    Default = cos()
    Hint = "This is the second function of this polar function."
    Visible = (@w_f == 2)
  Endfunc
  
  Heading
    Caption = "-------------------------------------------"
    Visible = (@w_f == 2)
  Endheading
 
 
  Heading
    Caption = "---Parameters for Func of Func---"
    Visible = (@w_f == 3)
  Endheading
  
  Param am_3_1
    Caption = "Amplitude 1"
    Default = 0.5
    Hint = "Amplitude of the external function of this polar shape."
    Visible = (@w_f == 3)
  Endparam
  
  Param am_3_2
    Caption = "Amplitude 2"
    Default = 2.0
    Hint = "Amplitude of the internal function of this polar shape."
    Visible = (@w_f == 3)
  Endparam
  
  Param fr_3_1
    Caption = "Frequency 1"
    Default = 4.0
    Hint = "This is the frequancy of the polar shape."
    Visible = (@w_f == 3)
  Endparam
  
  Param ph_3_1
    Caption = "Phase"
    Default = 0.0
    Hint = "Phase of this polar shape"
    Visible = (@w_f == 3)
  Endparam
  
  Func @wave_3_1
    Caption = "External Function"
    Default = sqrt()
    Hint = "External function of this shape"
    Visible = (@w_f == 3)
  Endfunc
  
  Func @wave_3_2
    Caption = "Internal Function"
    Default = cabs()
    Hint = "Internal fucntion of this shape"
    Visible = (@w_f == 3)
  Endfunc
  
  Heading
    Caption = "--------------------------------------"
    Visible = (@w_f == 3)
  Endheading
  
  Heading
    Caption = "--- Parameter for Division of 2 Func ---"
    Visible = (@w_f == 4)
  Endheading
  
  Param am_4_1
    Caption = "Amplitude 1"
    Default = 1.0
    Hint = "First amplitude value of this polar function."
    Visible = (@w_f == 4)
  Endparam
  
  Param fr_4_1
    Caption = "Frequency 1"
    Default = 4.0
    Hint = "This is the first frequency value of this polar function."
    Visible = (@w_f == 4)
  Endparam
  
  Param am_4_2
    Caption = "Amplitude 2"
    Default = 1.0
    Hint = "This is the second amplitude value of this polar function."
    Visible = (@w_f == 4)
  Endparam
  
  Param fr_4_2
    Caption = "Frequency 2"
    Default = 10.0
    Hint = "I'd like to know if you can find the meaning of this parameter :-)"
    Visible = (@w_f == 4)
  Endparam
  
  Param ph_4_1
    Caption = "Phase"
    Default = 0.0
    Hint = "This parameter specifies phase between two function of this shape."
    Visible = (@w_f == 4)
  Endparam
  
  Func wave_4_1
    Caption = "Wave Function 1"
    Default = recip()
    Hint = "This is the first function of this wave trap"
    Visible = (@w_f == 4)
  Endfunc
  
  Func wave_4_2
    Caption = "Wave Function 2"
    Default = sin()
    Hint = "This is the second function of this polar function."
    Visible = (@w_f == 4)
  Endfunc
  
  Heading
    Caption = "-------------------------------------------"
    Visible = (@w_f == 4)
  Endheading
  
  Heading
    Caption = "---Slices---"
    Visible = (@w_f == 5)
  Endheading
  
  Param s_a
    Caption = "N° of slices"
    Default = 8.0
    Hint = "N° of the slices of this shape."
    Visible = (@w_f == 5)
  Endparam 
  
  Param r_1
    Caption = "Ray 1"
    Default = 1.0
    Hint = "First ray of the shape"
    Visible = (@w_f == 5)
  Endparam
  
  Param r_2
    Caption = "Ray 2"
    Default = 0.5
    Hint = "Second ray of the shape"
    Visible = (@w_f == 5)
  Endparam
  
  Heading
    Caption = "---------------"
    Visible = (@w_f == 5)
  Endheading

  
  Param cn
    Caption = "Center"
    Default = (0.0,0.0)
    Hint = "This is the center of the Trap Shape"
  Endparam
  
  Param sc
    Caption = "Scale"
    Default = 1.0
    Hint = "Scale of the polar function"
  Endparam
  
  Param rot
    Caption = "Rotation"
    Default = 0.0
    Hint = "This is the angle rotation (in degrees) of the trap shape."
  Endparam
  
  Param a_r
    Caption = "Aspect Ratio"
    Default = 1.0
    Hint = "This parameter streches your trap shape"
  Endparam
       
  Param th
    Caption = "Thickness"
    Default = 1.0
    Hint = "Changing it you can change thickness if this shape"
  Endparam
  
  Param i_c
    Caption = "Iter. Control"
    Default = 21.5
    Hint = "This is the control of iterations. Use it to change \
            the range of values in which iterations work to change colors."
    Visible = (@t_t == 7)
  Endparam
  
  Param s_t
    Caption = "Solid?"
    Default = False
    Hint = "This value specifies what kind of solid color coloring \
            methid you want. In this way you can easily create frames \
            or other things."
    Visible = (@r_d == True)
  Endparam
  
; -------------------------------------------------------------------
; Advanced option Parameters
; ------------------------------------------------------------------- 

param adv
caption = "Advanced options"
hint = "Additional shaping parameters and functions."
default = false
endparam

param change
caption = "Morph"
enum = "1""2""3""4""5""6""7""8""9""10""11"
default = 0
visible = @adv
endparam

complex param exp1m
caption = "Bias 1"
default = (1,0)
visible = @adv
endparam

complex param exp3m
caption = "Bias 2"
default = (1,0)
visible = @adv && @change > 0
endparam

complex param exp2m
caption = "Bias 3"
default = (1,0)
visible = @adv && @change > 5
endparam

complex param tw
caption = "Twist 1"
default = (0,0)
visible = @adv
endparam

complex param tw2
caption = "Twist 2"
default = (0,0)
visible = @adv && @change > 0
endparam

complex param tw3
caption = "Twist 3"
default = (0,0)
visible = @adv && @change > 5
endparam

param opm
caption = "Operator"
enum = "+""-""*""/""^"
default = 0
visible = @adv && @change > 5
endparam

heading
caption = "Trap Functions"
visible = @adv
endheading

func nufunc
caption = "Morph function 1"
default = ident()
visible = @adv
endfunc

func nufunc2
caption = "Morph function 2"
default = ident ()
visible = @adv && @change > 0
endfunc



func nufunc3
caption = "Morph function 3"
default = ident ()
visible = @adv && @change > 5
endfunc

  
}