comment {
This file is the concatenation of

gwfa.ufm of October 21, 2000
gwfa-attractor.ufm of February 3, 2001
gwfa-matrix.ufm of June 1, 2001


into one file as gwfa.ufm  for purposes of Ultra Fractal Formula Database


}

;----------------------------------------------------------------------------------------------------------

comment {
  gwfa.ufm  Vers. 1.00.10a
  Fractal types for Ultra Fractal 2
  by G.W.F. Albrecht
  October 21, 2000

  For more information about this formula collection,
  please visit its home page:
  
      http://home.wtal.de/spiriteye/fractal/index.html

or

     http://www.gwfa.de


  All formulae are based on the function Phi(z):=(exp(z)-1)/z for |z| > 0
  or its reciprocal; in some cases - with transformations applied to z -
  the point z=0 has to be handled individually which has not been
  implemented consequently into all variants

  All algebraic operations take place in the algebras of complex numbers or Quaternions or complex 2x2 matrices
  equipped with the usual algebra norms. Quaternions are visualized by a 2dimensional cut through the origin,
  complex matrices by a simple entry. (Here is where your fantasy may come in!)

  While most of the formulae I used individually are now formally  included in this parametrizable set I made
  no attempt to reproduce any  image (if I have to render it I'll use the old formula).
  Some variants still have to be added.

  For compatibility reasons some older versions of formulas have been added at the end of this file (with
 titles starting with z-).
  However for future use the equivalent form of the newer versions should be used.
  Hence for new images do not use formulas below the line ...WwWwWw....


All formulas usually have a Mandelbroit and a Julia variant to be selected using the switch feature of UF.
Note that -  as not all parameters are propagated by switching - switching twice does not yield identity.

}



gwfa-Newton_1 {
; That is the simple Newton algorithm for Phi - 1 around 0
init:
complex  z = pixel 
loop:
complex u = exp(z), complex v=z*u , complex w = v - u + 1
  z = z - ( v - 2 * z ) / w
bailout:
  |z| <= @bailout
default:
  title = "Newton 1"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam


switch:
  type = "gwfa-Newton_1_Jul"
  seed = #pixel
  
  bailout = @bailout

}


gwfa-Newton_1_Jul {
; 
init:
complex  z = pixel 
loop:
complex u = exp(z), complex v=z*u , complex w = v - u + 1
  z = z - ( v - 2 * z ) / w + @seed
bailout:
  |z| <= @bailout
default:
  title = "Newton 1 Julia Variant"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

 param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam

switch:
  type = "gwfa-Newton_1"

  
  bailout = @bailout


}


gwfa-Newton_1_relaxed {
; That is the relaxed  Newton algorithm for Phi  around 0
init:
complex  z = pixel 
complex s2 = @sfact*@sfact
loop:
complex u = exp(z)
complex u2 = u * u
complex z1 = z - 1
complex z12 = z1 * z1
complex v= z * ( u2 * z1 * ( z1 - @rfact) + u * ( z * ( @rfact + 2) - 2*( @rfact + 1) ) + s2 * z^2 + @rfact + 1 ) 
complex w =  u2 *  z12 + 2* u*z1 + s2*z^4+1
  z =  v / w

bailout:
  |z| <= @bailout

default:
  title = "Newton 1 Relaxed  Phi=1"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param rfact
    caption = "Relaxation Factor"
    default = (1,0)
  endparam

  param sfact
    caption = "Pert.Relaxation Factor "
    default = (0,0)
  endparam


switch:
  type = "gwfa-Newton_1_relaxed_Jul"
  seed = #pixel
  rfact = @rfact
  sfact = @sfact
  bailout = @bailout

}

gwfa-Newton_1_relaxed_Jul {
; That is the relaxed  Newton algorithm for Phi  around 0 - Julia variant
init:
complex  z = pixel 
complex s2 = @sfact*@sfact
loop:
complex u = exp(z)
complex u2 = u * u
complex z1 = z - 1
complex z12 = z1 * z1
complex v= z * ( u2 * z1 * ( z1 - @rfact) + u * ( z * ( @rfact + 2) - 2*( @rfact + 1) ) + s2 * z^2 + @rfact + 1 ) 
complex w =  u2 *  z12 + 2* u*z1 + s2*z^4+1
  z =  v / w + @seed
bailout:
  |z| <= @bailout
default:
  title = "Newton 1Relaxed Julia Phi=1"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param rfact
    caption = "Relaxation Factor"
    default = (1,0)
  endparam

  param sfact
    caption = "Pert.Relaxation Factor "
    default = (0,0)
  endparam

param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam


switch:
  type = "gwfa-Newton_1_relaxed"
    rfact = @rfact
  sfact = @sfact
  bailout = @bailout

}


gwfa-Newton_2_relaxed {
; That is the relaxed  Newton algorithm for Phi  around 0
init:
complex  z = pixel 
complex s2 = @sfact*z*z*(0.0,1.0)
loop:
complex u = exp(z)
complex v = z *( u* (z -@rfact -1) + @rfact * z + @rfact + 1 +s2 )
complex w = u * ( z - 1) + 1 + s2
  z =  v / w

bailout:
  |z| <= @bailout

default:
  title = "Newton 2 Relaxed Phi=0"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param rfact
    caption = "Relaxation Factor"
    default = (1,0)
  endparam

  param sfact
    caption = "Pert.Relaxation Factor "
    default = (0,0)
  endparam


switch:
  type = "gwfa-Newton_2_relaxed_Jul"
  seed = #pixel
  rfact = @rfact
  sfact = @sfact
  bailout = @bailout

}

gwfa-Newton_2_relaxed_Jul {
; That is the relaxed  Newton algorithm for Phi - 1 around 0 - Julia variant
init:
complex  z = pixel 
complex s2 = @sfact*z*z*(0.0,1.0)
loop:
complex u = exp(z)
complex v = z *( u* (z -@rfact -1) + @rfact * z + @rfact + 1 +s2 )
complex w = u * ( z - 1) + 1 + s2
  z =  v / w + @seed

bailout:
  |z| <= @bailout
default:
  title = "Newton 2 Relaxed Phi=0 Julia variant"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param rfact
    caption = "Relaxation Factor"
    default = (1,0)
  endparam

  param sfact
    caption = "Pert.Relaxation Factor "
    default = (0,0)
  endparam

param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam


switch:
  type = "gwfa-Newton_2_relaxed"
    rfact = @rfact
  sfact = @sfact
  bailout = @bailout

}







gwfa-Newton_3 {
; another variant
init:
complex  z = pixel 
loop:

complex  u = exp(z)
complex v=z*u 
complex w = v - u + 1 
  z =  z * ( v + 2 * ( 1 + z - u ) ) / w
bailout:
  |z| <= @bailout
default:
  title = "Newton 3"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam


switch:
  type = "gwfa-Newton_3_Jul"
  seed = #pixel
  
  bailout = @bailout

}


gwfa-Newton_3_Jul {
; 
init:
complex  z = pixel 
loop:
complex  u = exp(z)
complex v=z*u 
complex w = v - u + 1 
  z =  z * ( v + 2 * ( 1 + z - u ) ) / w +@seed

bailout:
  |z| <= @bailout
default:
  title = "Newton 3 Julia Variant"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

 param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam

switch:
  type = "gwfa-Newton_3"

  
  bailout = @bailout


}




gwfa-Newton_4_relaxed {
; That is the relaxed  Newton algorithm for the first derivative of Phi  around 0
init:
complex  z = pixel 
complex s2 = @sfact*z*z*z*(0.0,1.0)
complex s3 = 2*(@rfact +2)  
loop:
complex u = exp(z)
complex v = z *( 2*u* (z^2 -z*(@rfact+2) + @rfact  +2 )  + @rfact*z*z - s3  +s2 )
complex w = 2*(u * ( z*z - 2*z +2) - s2 )
  z =  v / w

bailout:
  |z| <= @bailout

default:
  title = "Newton 4 Relaxed Phi'=1/2"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param rfact
    caption = "Relaxation Factor"
    default = (1,0)
  endparam

  param sfact
    caption = "Pert.Relaxation Factor "
    default = (0,0)
  endparam


switch:
  type = "gwfa-Newton_4_relaxed_Jul"
  seed = #pixel
  rfact = @rfact
  sfact = @sfact
  bailout = @bailout

}

gwfa-Newton_4_relaxed_Jul {
; That is the relaxed  Newton algorithm for the first derivative of  Phi  around 0 - Julia variant

init:
complex  z = pixel 
complex s2 = @sfact*z*z*z*(0.0,1.0)
complex s3 = 2*(@rfact +2)  
loop:
complex u = exp(z)
complex v = z *( 2*u* (z^2 -z*(@rfact+2) + @rfact  +2 )  + @rfact*z*z - s3  +s2 )
complex w = 2*(u * ( z*z - 2*z +2) - s2 )
  z =  v / w + @seed
bailout:
  |z| <= @bailout
default:
  title = "Newton 4 Relaxed Phi'=1/2 Julia variant"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param rfact
    caption = "Relaxation Factor"
    default = (1,0)
  endparam

  param sfact
    caption = "Pert.Relaxation Factor "
    default = (0,0)
  endparam

param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam


switch:
  type = "gwfa-Newton_4_relaxed"
    rfact = @rfact
  sfact = @sfact
  bailout = @bailout

}


gwfa-Newton_5_relaxed {
; That is the relaxed  Newton algorithm for the reciprocal of Phi  around 0
init:
complex  z = pixel 
complex s2 = @sfact*(0.0,1.0)
loop:
complex u = exp(z)
complex uu = u^2
complex f = (u - 1)^2
complex v = @rfact*( uu - u*(z+2) + z + 1)
complex w = u * ( 1- z ) - 1 + s2*f
  z = z +  v / w

bailout:
  |z| <= @bailout

default:
  title = "Newton 5 Relaxed 1/Phi=1"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param rfact
    caption = "Relaxation Factor"
    default = (1,0)
  endparam

  param sfact
    caption = "Pert.Relaxation Factor "
    default = (0,0)
  endparam


switch:
  type = "gwfa-Newton_5_relaxed_Jul"
  seed = #pixel
  rfact = @rfact
  sfact = @sfact
  bailout = @bailout

}

gwfa-Newton_5_relaxed_Jul {
; That is the relaxed  Newton algorithm for the reciprocal  around 0 - Julia variant

init:
complex  z = pixel 
complex s2 = @sfact*(0.0,1.0)
loop:
complex u = exp(z)
complex uu = u^2
complex f = (u - 1)^2
complex v = @rfact*( uu - u*(z+2) + z + 1)
complex w = u * ( 1- z ) - 1 + s2*f
  z = z +  v / w + @seed


bailout:
  |z| <= @bailout
default:
  title = "Newton 5 Relaxed 1/Phi=1 Julia variant"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param rfact
    caption = "Relaxation Factor"
    default = (1,0)
  endparam

  param sfact
    caption = "Pert.Relaxation Factor "
    default = (0,0)
  endparam

param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam


switch:
  type = "gwfa-Newton_5_relaxed"
    rfact = @rfact
  sfact = @sfact
  bailout = @bailout

}



gwfa-Newton_6_relaxed {
; That is the relaxed  Newton algorithm for a convex combination of  Phi and its reciprocal around 0

init:
complex  z = pixel 
complex r1 = @rfact + 1
complex r11 = @rfact - 1
complex l1 = @lambda - 1
complex s2 = @sfact*(0.0,1.0)
complex f1 =  6*@lambda*r1
complex f2 = 4*@lambda*r1

loop:

complex u = exp(z)
complex uu = u^2
complex uuu = u*uu
complex u1 = (u - 1) ^3

complex v =  @lambda*uu*uu*( r1 -z )
v = v + @lambda * uuu * ( 3*z - 4*r1)
v = v - uu*( z^3* l1 + z^2*r11*l1 + 3*@lambda*z - f1)
v = v + u *( z^3 * l1 + 2* z^2*r11*l1 + @lambda * z - f2)
v = v  - z^2* r11*l1 * @lambda*r1 - s2*z^2*u1
v = z * v 
complex w = @lambda*uu*uu*(1-z)
w = w + @lambda*uuu*( 3*z - 4)
w = w - uu*(z^3*l1 - z^2*l1+ 3* @lambda*z - 6* @lambda)
w = w + u * (z^3*l1 - 2*z^2*l1+  @lambda*z - 4* @lambda)
w = w + z^2 *l1 + @lambda - s2*z^2*u1

  z =   v / w

bailout:
  |z| <= @bailout

default:
  title = "Newton 6 Relaxed Convex"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param rfact
    caption = "Relaxation Factor"
    default = (1,0)
  endparam

  param sfact
    caption = "Pert.Relaxation Factor "
    default = (0,0)
  endparam

  param lambda
    caption = "convexity "
    default = (0.5,0)
  endparam


switch:
  type = "gwfa-Newton_6_relaxed_Jul"
  seed = #pixel
  rfact = @rfact
  sfact = @sfact
  lambda = @lambda
  bailout = @bailout

}

gwfa-Newton_6_relaxed_Jul {
; That is the relaxed  Newton algorithm for a convex combination of  Phi and its reciprocal around 0 - Julia variant


init:
complex  z = pixel 
complex r1 = @rfact + 1
complex r11 = @rfact - 1
complex l1 = @lambda - 1
complex s2 = @sfact*(0.0,1.0)
complex f1 =  6*@lambda*r1
complex f2 = 4*@lambda*r1

loop:

complex u = exp(z)
complex uu = u^2
complex uuu = u*uu
complex u1 = (u - 1) ^3

complex v =  @lambda*uu*uu*( r1 -z )
v = v + @lambda * uuu * ( 3*z - 4*r1)
v = v - uu*( z^3* l1 + z^2*r11*l1 + 3*@lambda*z - f1)
v = v + u *( z^3 * l1 + 2* z^2*r11*l1 + @lambda * z - f2)
v = v  - z^2* r11*l1 * @lambda*r1 - s2*z^2*u1
v = z * v 
complex w = @lambda*uu*uu*(1-z)
w = w + @lambda*uuu*( 3*z - 4)
w = w - uu*(z^3*l1 - z^2*l1+ 3* @lambda*z - 6* @lambda)
w = w + u * (z^3*l1 - 2*z^2*l1+  @lambda*z - 4* @lambda)
w = w + z^2 *l1 + @lambda - s2*z^2*u1

  z =   v / w + @seed


bailout:
  |z| <= @bailout
default:
  title = "Newton 6 Relaxed Convex - Julia variant"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param rfact
    caption = "Relaxation Factor"
    default = (1,0)
  endparam

  param sfact
    caption = "Pert.Relaxation Factor "
    default = (0,0)
  endparam

  param lambda
    caption = "convexity "
    default = (0.5,0)
  endparam


param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam


switch:
  type = "gwfa-Newton_6_relaxed"
    rfact = @rfact
  sfact = @sfact
  lambda = @lambda
  bailout = @bailout

}






gwfa-Phi_1 {
; 
init:
complex   z = Pixel
loop:
  z = z + @shift1
complex  a = fn1(z)
 a = a + @shift2
IF ( |a| > 0 )
  z = (exp(a)-1)/a + Pixel*@PixelWeight
ELSE
  z = 1 + Pixel * @PixelWeight
ENDIF
 z = fn2(z)
bailout:
    |z| <= @bailout

default:
  title = "Phi Version 1"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0,0)
  endparam


  param shift2
    caption = "Pixel Shift 2"
    default = (0,0)
  endparam

  param PixelWeight
    caption = "Pixel Weight"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function"
    default = ident()
  endfunc  

 func fn2
    caption = "Post function"
    default = ident()
  endfunc  
  

switch:
  type = "gwfa-Phi_1_Jul"
  seed = #pixel
  fn1 = fn1
  fn2 = fn2
  shift1 = @shift1
  shift2 = @shift2
  
  bailout = @bailout



  }


gwfa-Phi_1_Jul {
; This is the Julia variant for the previous formula
init:
complex   z = Pixel
loop:
  z = z + @shift1
complex  a = fn1(z)
 a = a + @shift2
IF ( |a| > 0 )
  z = (exp(a)-1)/a + @seed
ELSE
  z = 1 + @seed
ENDIF
 z = fn2(z)
bailout:
    |z| <= @bailout

default:
  title = "Phi Version 1Julia Variant"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0,0)
  endparam


  param shift2
    caption = "Pixel Shift 2"
    default = (0,0)
  endparam

  param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function"
    default = ident()
  endfunc  

 func fn2
    caption = "Post function"
    default = ident()
  endfunc  


switch:
  type = "gwfa-Phi_1"

  fn1 = fn1
  fn2 = fn2
  shift1 = @shift1
  shift2 = @shift2
  
  bailout = @bailout

  
  }








gwfa-Phi_2 {
; 
init:
complex z = Pixel
loop:
  z = z + @shift1
complex  a = fn1(z)
 a = a + @shift2

complex b = (exp(a)-1)

IF ( |b| > 0 )
  z = fn2(a)/b + Pixel*@PixelWeight
ELSE
  z = 1 + Pixel * @PixelWeight
; this branch can be dubious !!
ENDIF
 z = fn3(z)
bailout:
    |z| <= @bailout

default:
  title = "Phi Version 2 recip"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0,0)
  endparam


  param shift2
    caption = "Pixel Shift 2"
    default = (0,0)
  endparam

  param PixelWeight
    caption = "Pixel Weight"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc  

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc  

 func fn3
    caption = "Post function"
    default = ident()
  endfunc  
  

switch:
  type = "gwfa-Phi_2_Jul"
  seed = #pixel
  fn1 = fn1
  fn2 = fn2
 fn3 = fn3
 PixelWeight = @PixelWeight
  shift1 = @shift1
  shift2 = @shift2
  
  bailout = @bailout



  }


gwfa-Phi_2_Jul {
; 
init:
complex z = Pixel
loop:
  z = z + @shift1
complex  a = fn1(z)
 a = a + @shift2

complex b = (exp(a)-1)

IF ( |b| > 0 )
  z = fn2(a)/b + @seed*@PixelWeight
ELSE
  z = 1 + @seed * @PixelWeight
; this branch can be dubious !!
ENDIF
 z = fn3(z)
bailout:
    |z| <= @bailout

default:
  title = "Phi Version 2 recip Julia"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0,0)
  endparam


  param shift2
    caption = "Pixel Shift 2"
    default = (0,0)
  endparam

  param PixelWeight
    caption = "Pixel Weight"
    default = (1,0)
  endparam

  param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam



 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc  

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc  

 func fn3
    caption = "Post function"
    default = ident()
  endfunc  

switch:
  type = "gwfa-Phi_2"

  fn1 = fn1
  fn2 = fn2
 fn3 = fn3
 PixelWeight = @PixelWeight
  shift1 = @shift1
  shift2 = @shift2
  
  bailout = @bailout


  

  }



gwfa-Phi_3 {
; 
init:
complex  z = Pixel
loop:
  z = z + @shift1
complex a = fn1(z)
 a = a + @shift2

complex b = exp(a)


  z = fn2(a)/(b -1) + Pixel*@PixelWeight
; this line can be dubious !!

 z = fn3(z)
bailout:
    |z| <= @bailout

default:
  title = "Phi Version 3 recip"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0,0)
  endparam


  param shift2
    caption = "Pixel Shift 2"
    default = (0,0)
  endparam

  param PixelWeight
    caption = "Pixel Weight"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc  

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc  

 func fn3
    caption = "Post function"
    default = ident()
  endfunc  
  

switch:
  type = "gwfa-Phi_3_Jul"
  seed = #pixel
  fn1 = fn1
  fn2 = fn2
 fn3 = fn3
 PixelWeight = @PixelWeight
  shift1 = @shift1
  shift2 = @shift2
  
  bailout = @bailout




  }



gwfa-Phi_3_Jul {
; 
init:
complex  z = Pixel
loop:
  z = z + @shift1
complex a = fn1(z)
 a = a + @shift2

complex b = exp(a)


  z = fn2(a)/(b -1) + @seed*@PixelWeight
; this line can be dubious !!

 z = fn3(z)
bailout:
    |z| <= @bailout

default:
  title = "Phi Version 3 recip Julia"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0,0)
  endparam


  param shift2
    caption = "Pixel Shift 2"
    default = (0,0)
  endparam

  param PixelWeight
    caption = "Pixel Weight"
    default = (1,0)
  endparam

  param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam


 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc  

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc  

 func fn3
    caption = "Post function"
    default = ident()
  endfunc  
  


switch:
  type = "gwfa-Phi_3"

  fn1 = fn1
  fn2 = fn2
 fn3 = fn3
 PixelWeight = @PixelWeight
  shift1 = @shift1
  shift2 = @shift2
  
  bailout = @bailout



  }



;The following formulae have a  more incidental character and emerged while I looked for the
lost original Phi_4 (FractInt version)



gwfa-Phi_4_trial {; GA
init:
complex  z = Pixel 
loop:
complex  a = fn1(z)
  z = z + @shift1 
complex   b = exp(fn2(z))
  z = a/(b-1)   +  Pixel*@PixelWeight
bailout:
    |z| <= @bailout  

default:
  title = "Phi Version 4 recip"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0.1,0)
  endparam


  param PixelWeight
    caption = "Pixel Weight"
    default = (1,0)
  endparam




 func fn1
    caption = "Pre function1"
    default = sqr()
  endfunc  

 func fn2
    caption = "Pre function2"
    default = sqr()
  endfunc  


switch:
  type = "gwfa-Phi_4_trial_Jul"
  seed = #pixel
  shift1 = @shift1
  fn1 = fn1
  fn2 = fn2
  
  bailout = @bailout

  

  }




gwfa-Phi_4_trial_Jul {; GA
init:
complex  z = Pixel 
loop:
complex  a = fn1(z)
  z = z + @shift1 
complex   b = exp(fn2(z))
  z = a/(b-1)   +  @seed
bailout:
    |z| <= @bailout  

default:
  title = "Phi Version 4 recip Julia"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0.1,0)
  endparam

 param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam


 func fn1
    caption = "Pre function1"
    default = sqr()
  endfunc  

 func fn2
    caption = "Pre function2"
    default = sqr()
  endfunc  

switch:
  type = "gwfa-Phi_4_trial"

  shift1 = @shift1
  fn1 = fn1
  fn2 = fn2
  
  bailout = @bailout


  

  }




gwfa-Der1 {
; first derivative of Phi
init:
complex   z = Pixel 
complex c = Pixel 
loop:
  z = fn1(z)
  a = exp(z)
  b = z * z
IF ( |z| > 0 )
  z = ( (z-1)*a + 1 ) / b - 0.5 + @PixelWeight * c
ELSE
  z = @PixelWeight * c
ENDIF

z = fn2(z)

bailout:
    |z| <= @bailout  

default:
  title = "1st derivative of Phi"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam


  param PixelWeight
    caption = "Pixel Weight"
    default = (1,0)
  endparam


 func fn1
    caption = "Pre function"
    default = ident()
  endfunc  

 func fn2
    caption = "Post function"
    default = ident()
  endfunc  

switch:
  type = "gwfa-Der1_Jul"
  seed = #pixel
  shift1 = @shift1
  fn1 = fn1
  fn2 = fn2
 PixelWeight = @PixelWeight
  
  bailout = @bailout

}


gwfa-Der1_Jul {
; first derivative of Phi /Julia variant
complex   z = Pixel 

loop:
  z = fn1(z)
  a = exp(z)
  b = z * z
IF ( |z| > 0 )
  z = ( (z-1)*a + 1 ) / b - 0.5 +  @PixelWeight *@seed
ELSE
  z = @PixelWeight * @seed
ENDIF

z = fn2(z)

bailout:
    |z| <= @bailout  

default:
  title = "1st derivative of Phi / Julia"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam


  param PixelWeight
    caption = "Pixel Weight"
    default = (1,0)
  endparam

 param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam



 func fn1
    caption = "Pre function"
    default = ident()
  endfunc  

 func fn2
    caption = "Post function"
    default = ident()
  endfunc  


switch:
  type = "gwfa-Der1"

  shift1 = @shift1
  fn1 = fn1
  fn2 = fn2
 PixelWeight = @PixelWeight
  
  bailout = @bailout



}



gwfa-Phi_4_trial2 {; GA
init:
complex  z = Pixel 
loop:
complex  a = fn1(z)
  z = z + @shift1 
complex   b = exp(fn2(z))
  z = (b-1)/a   +  Pixel*@PixelWeight
bailout:
    |z| <= @bailout  

default:
  title = "Phi Version 4 "
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0.1,0)
  endparam


  param PixelWeight
    caption = "Pixel Weight"
    default = (1,0)
  endparam




 func fn1
    caption = "Pre function1"
    default = sqr()
  endfunc  

 func fn2
    caption = "Pre function2"
    default = sqr()
  endfunc  


switch:
  type = "gwfa-Phi_4_trial2_Jul"
  seed = #pixel
  shift1 = @shift1
  fn1 = fn1
  fn2 = fn2
  
  bailout = @bailout

  

  }




gwfa-Phi_4_trial2_Jul {; GA
init:
complex  z = Pixel 
loop:
complex  a = fn1(z)
  z = z + @shift1 
complex   b = exp(fn2(z))
  z = (b-1)/a   +  @seed
bailout:
    |z| <= @bailout  

default:
  title = "Phi Version 4  Julia"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0.1,0)
  endparam

 param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam


 func fn1
    caption = "Pre function1"
    default = sqr()
  endfunc  

 func fn2
    caption = "Pre function2"
    default = sqr()
  endfunc  

switch:
  type = "gwfa-Phi_4_trial2"
  shift1 = @shift1
  fn1 = fn1
  fn2 = fn2
  
  bailout = @bailout

  

  }

; following are some arbitrary tries for new effects

gwfa-Quot_1 {; GA
init:
complex   z = Pixel
loop:
  z = z + @shift1
  a = fn1(z)
  b = fn2(z)
  c = exp(a) - 1
  d = exp(b) - 1
  z = c / d - 1
bailout:
    |z| <= @bailout  

default:
  title = "Phi Quotients"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0.1,0)
  endparam

 func fn1
    caption = "Pre function1"
    default = sqr()
  endfunc  

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc  

switch:
  type = "gwfa-Quot_1_Jul"
  seed = #pixel
  shift1 = @shift1
  fn1 = fn1
  fn2 = fn2
  
  bailout = @bailout


 
  }

gwfa-Quot_1_Jul {; GA
init:
complex   z = Pixel
loop:
  z = z + @shift1
  a = fn1(z)
  b = fn2(z)
  c = exp(a) - 1
  d = exp(b) - 1
  z = c / d - 1 +@seed
bailout:
    |z| <= @bailout  

default:
  title = "Phi Quotients Julia"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0.1,0)
  endparam

 param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam



 func fn1
    caption = "Pre function1"
    default = sqr()
  endfunc  

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc  

switch:
  type = "gwfa-Quot_1"

  shift1 = @shift1
  fn1 = fn1
  fn2 = fn2
  
  bailout = @bailout

 
  }



gwfa-Phi_9_split {; GA
init:
 complex  z = pixel
loop:
 z  = z + @shift1
  x = real(z)
 y = imag(z)
 complex a = fn1(x) 
 complex  b = fn2(y)
 
IF ( |a| > 0 )
complex  c = (exp(a)-1)/a 
ELSE
  c = 1 
ENDIF

IF ( |b| > 0 )
complex  d = (exp(b)-1)/b 
ELSE
  d = 1 
ENDIF

z = fn3(  ( c + flip(d)    ) * fn4(z) ) + pixel
bailout:
    |z| <= @bailout 


default:
  title = "Real / Im Split "
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0.1,0)
  endparam


 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc  

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc  

 func fn3
    caption = "Iter function 1"
    default = exp()
  endfunc  

 func fn4
    caption = "Iter function2"
    default = log()
  endfunc  

switch:
  type = "gwfa-Phi_9_split_Jul"
  seed = #pixel
  shift1 = @shift1
  fn1 = fn1
  fn2 = fn2
  fn3 = fn3
  fn4 = fn4

  
  bailout = @bailout


 
  }


  
gwfa-Phi_9_split_Jul {; GA
init:
 complex  z = pixel
loop:
 z  = z + @shift1
  x = real(z)
 y = imag(z)
complex a = fn1(x) 
complex  b = fn2(y)
 
IF ( |a| > 0 )
complex  c = (exp(a)-1)/a 
ELSE
  c = 1 
ENDIF

IF ( |b| > 0 )
complex  d = (exp(b)-1)/b 
ELSE
  d = 1 
ENDIF

z = fn3(  ( c + flip(d)    ) * fn4(z) ) + @seed
bailout:
    |z| <= @bailout 


default:
  title = "Real / Im Split Julia"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param shift1
    caption = "Pixel Shift 1"
    default = (0.1,0)
  endparam


 param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam




 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc  

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc  

 func fn3
    caption = "Iter function 1"
    default = exp()
  endfunc  

 func fn4
    caption = "Iter function2"
    default = log()
  endfunc  


switch:
  type = "gwfa-Phi_9_split"

  shift1 = @shift1
  fn1 = fn1
  fn2 = fn2
  fn3 = fn3
  fn4 = fn4

  
  bailout = @bailout



 
  }


; the following formulae implement Phi1 in the real Quaternion algebra
; the first examples use
; no closed expression for Phi but partial sums (ie n=3 or n=7) 
; whereas quatphi implements (exp(z)-1)/z in the sense of the quaternions


;

gwfa-c-Quat_1 {
; Quaternionic version of Phi1 (with partial sum for n=3)
init:
complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(z))
float d = imag(fn1(z))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(z))
 d = imag(fn1(z))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(z))
 b = imag(fn1(z))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(z))
 c = imag(fn1(z))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(z))
 d = imag(fn1(z))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(z))
 c = imag(fn1(z))


endif



loop:

float t1 = (3*a^2+8*a-b^2-c^2-d^2 +12)/24


float anew = (a^3+4*a^2-3*a*(b^2+c^2+d^2-4)-4*(b^2+c^2+d^2-6))/24 
float bnew = b*t1
float cnew = c*t1
float dnew = d*t1

a = anew
b = bnew
c = cnew
d = dnew



if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  =  b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

bailout:
  |z| <= @bailout
default:
  title = "Quaternionic Phi  1 - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-c-Quat_1_Jul"
  seed = #pixel
  plane = @plane  
  bailout = @bailout
  fn1 = fn1
}

gwfa-c-Quat_1_Jul {
; Quaternionic version of Phi1 (with partial sum for n=3) - Julia variant
init:

complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(z))
float d = imag(fn1(z))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(z))
 d = imag(fn1(z))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(z))
 b = imag(fn1(z))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(z))
 c = imag(fn1(z))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(z))
 d = imag(fn1(z))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(z))
 c = imag(fn1(z))


endif



float u = real(@seed)
float v = imag(@seed)

loop:


float t1 = (3*a^2+8*a-b^2-c^2-d^2 +12)/24

float anew = (a^3+4*a^2-3*a*(b^2+c^2+d^2-4)-4*(b^2+c^2+d^2-6))/24 + u
float bnew = b*t1 + v
float cnew = c*t1 + u 
float dnew = d*t1 + v

a = anew
b = bnew
c = cnew
d = dnew

if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  = b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif


bailout:

|z| <= @bailout

default:
  title = "Quaternionic Phi  1 Julia  - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 

  endparam
 



  param seed
    caption = "fixed perturbation"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-c-Quat_1"
    plane = @plane
  bailout = @bailout
  fn1 = fn1

}

gwfa-c-Quat_1a {
; Quaternionic version of Phi1 (with partial sum for n=3)
; bailout following the quaternion norm
init:
complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 b = imag(fn1(@dim34))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(@dim34))
 c = imag(fn1(@dim34))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 c = imag(fn1(@dim34))


endif


float r = real(@bailout^2)

loop:

float t1 = (3*a^2+8*a-b^2-c^2-d^2 +12)/24


float anew = (a^3+4*a^2-3*a*(b^2+c^2+d^2-4)-4*(b^2+c^2+d^2-6))/24 
float bnew = b*t1
float cnew = c*t1
float dnew = d*t1

a = anew
b = bnew
c = cnew
d = dnew



if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  =  b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float  rho = a^2 +b^2+c^2+d^2

bailout:
rho <= r

default:
  title = "Quaternionic Phi  1a - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-c-Quat_1a_Jul"
  seed = #pixel
  plane = @plane  
  bailout = @bailout
  dim34 = @dim34
  fn1 = fn1
}

gwfa-c-Quat_1a_Jul {
; Quaternionic version of Phi1 (with partial sum for n=3) - Julia variant
; bailout following the quaternion norm
init:
complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 b = imag(fn1(@dim34))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(@dim34))
 c = imag(fn1(@dim34))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 c = imag(fn1(@dim34))


endif


float r = real(@bailout^2)

float u = real(@seed)
float v = imag(@seed)

loop:


float t1 = (3*a^2+8*a-b^2-c^2-d^2 +12)/24

float anew = (a^3+4*a^2-3*a*(b^2+c^2+d^2-4)-4*(b^2+c^2+d^2-6))/24 + u
float bnew = b*t1 + v
float cnew = c*t1 + u 
float dnew = d*t1 + v

a = anew
b = bnew
c = cnew
d = dnew

if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  = b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float rho = a^2 + b^2 + c^2 + d^2

bailout:
rho <= r
default:
  title = "Quaternionic Phi  1a Julia  - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 

  endparam
 



  param seed
    caption = "fixed perturbation"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-c-Quat_1a"
    plane = @plane
  bailout = @bailout
  dim34 = @dim34
  fn1 = fn1

}



gwfa-c-Quat_2 {
; Quaternionic version of Phi1 (with partial sum for n=7)
init:

complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(z))
float d = imag(fn1(z))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(z))
 d = imag(fn1(z))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(z))
 b = imag(fn1(z))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(z))
 c = imag(fn1(z))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(z))
 d = imag(fn1(z))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(z))
 c = imag(fn1(z))


endif


loop:


float t1 = (3*b^2+3*c^2+3*d^2-8)
float t2 = (5*b^2+5*c^2+5*d^2-14)
float t3 = (b^4+2*b^2*(c^2+d^2-8)+c^4+2*c^2*(d^2-8)+d^4-16*d^2+48)
float t4 = (5*b^4+2*b^2*(5*c^2+5*d^2-42)+5*c^4+2*c^2*(5*d^2-42)+5*d^4-84*d^2+280)
float t5 = (c^4+2*c^2*(d^2-14)+d^4-28*d^2+280)

float t6 =(b^2+c^2+d^2-8)
float t7 = (5*b^2+5*c^2+5*d^2-42)
float t8 = (3*b^4+2*b^2*(3*c^2+3*d^2-40)+3*c^4+2*c^2*(3*d^2-40)+3*d^4-80*d^2+720)
float t9 = (b^4+2*b^2*(c^2+d^2-14)+c^4+2*c^2*(d^2-14)+d^4-28*d^2+280)
float t10 = (3*c^4+2*c^2*(3*d^2-56)+3*d^4-112*d^2+1680)
float t11 = (3*d^4-112*d^2+1680)
float t13 = (3*c^2+3*d^2-40)
float t14 = (3*c^4+2*c^2*(3*d^2-40)+3*d^4-80*d^2+720)
float t15 = (3*d^4-80*d^2+720)
float t16 = (b^6+3*b^4*(c^2+d^2-14)+3*b^2*t5 +c^6+3*c^4*(d^2-14)+3*c^2*(d^4-28*d^2+280)+d^6-42*d^4+840*d^2-5040)

float t17=(7*a^6+48*a^5-35*a^4*t6-32*a^3*t7+7*a^2*t8+48*a*t9-b^6-b^4*(3*c^2+3*d^2-56)-b^2*t10-c^6+c^4*(56-3*d^2)-c^2*t11-d^6+56*d^4-1680*d^2+20160)/40320



float anew = (a^7+8*a^6-7*a^5*t1-24*a^4*t2+35*a^3*t3+24*a^2*t4-7*a*(b^6+b^4*t13+b^2*t14+c^6+c^4*(3*d^2-40)+c^2*t15+d^6-40*d^4+720*d^2-2880)-8*t16)/40320


float bnew  = b*t17

float cnew  = c*t17

float dnew  = d*t17

a = anew
b = bnew
c = cnew
d = dnew


if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  =  b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

bailout:
  |z| <= @bailout
default:
  title = "Quaternionic Phi  2 - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 1
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc



switch:
  type = "gwfa-c-Quat_2_Jul"
  seed = #pixel
  plane = @plane  
  bailout = @bailout
  fn1 = fn1

}

gwfa-c-Quat_2_Jul {
; Quaternionic version of Phi1 (with partial sum for n=7) - Julia variant
init:

complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(z))
float d = imag(fn1(z))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(z))
 d = imag(fn1(z))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(z))
 b = imag(fn1(z))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(z))
 c = imag(fn1(z))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(z))
 d = imag(fn1(z))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(z))
 c = imag(fn1(z))


endif




float u = real(@seed)
float v = imag(@seed)
float  w = real(@seed)
float s = imag(@seed)


loop:



float t1 = (3*b^2+3*c^2+3*d^2-8)
float t2 = (5*b^2+5*c^2+5*d^2-14)
float t3 = (b^4+2*b^2*(c^2+d^2-8)+c^4+2*c^2*(d^2-8)+d^4-16*d^2+48)
float t4 = (5*b^4+2*b^2*(5*c^2+5*d^2-42)+5*c^4+2*c^2*(5*d^2-42)+5*d^4-84*d^2+280)
float t5 = (c^4+2*c^2*(d^2-14)+d^4-28*d^2+280)

float t6 =(b^2+c^2+d^2-8)
float t7 = (5*b^2+5*c^2+5*d^2-42)
float t8 = (3*b^4+2*b^2*(3*c^2+3*d^2-40)+3*c^4+2*c^2*(3*d^2-40)+3*d^4-80*d^2+720)
float t9 = (b^4+2*b^2*(c^2+d^2-14)+c^4+2*c^2*(d^2-14)+d^4-28*d^2+280)
float t10 = (3*c^4+2*c^2*(3*d^2-56)+3*d^4-112*d^2+1680)
float t11 = (3*d^4-112*d^2+1680)
float t13 = (3*c^2+3*d^2-40)
float t14 = (3*c^4+2*c^2*(3*d^2-40)+3*d^4-80*d^2+720)
float t15 = (3*d^4-80*d^2+720)
float t16 = (b^6+3*b^4*(c^2+d^2-14)+3*b^2*t5 +c^6+3*c^4*(d^2-14)+3*c^2*(d^4-28*d^2+280)+d^6-42*d^4+840*d^2-5040)

float t17=(7*a^6+48*a^5-35*a^4*t6-32*a^3*t7+7*a^2*t8+48*a*t9-b^6-b^4*(3*c^2+3*d^2-56)-b^2*t10-c^6+c^4*(56-3*d^2)-c^2*t11-d^6+56*d^4-1680*d^2+20160)/40320


float anew = (a^7+8*a^6-7*a^5*t1-24*a^4*t2+35*a^3*t3+24*a^2*t4-7*a*(b^6+b^4*t13+b^2*t14+c^6+c^4*(3*d^2-40)+c^2*t15+d^6-40*d^4+720*d^2-2880)-8*t16)/40320

float bnew  = b*t17
float cnew = c*t17
float dnew  = d*t17



a = anew + u
b = bnew + v
c = cnew + w
d = dnew + s


if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  = b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif



bailout:
  |z| <= @bailout
default:
  title = "Quaternionic Phi  2 Julia  - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 1
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 

  endparam
 



  param seed
    caption = "fixed perturbation"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-c-Quat_2"
    plane = @plane
  bailout = @bailout
  fn1 = fn1

}





gwfa-c-QuatPhi1 {
; full Quaternionic version of Phi1 
; bailout following the quaternion norm
init:

complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 b = imag(fn1(@dim34))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(@dim34))
 c = imag(fn1(@dim34))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 c = imag(fn1(@dim34))


endif



float r = real(@bailout)^2
loop:



float norm1 = a^2+b^2+c^2+d^2
float norm2 = SQRT(b^2 + c^2 + d^2)
float s1 = SIN(norm2)
float c1 = COS(norm2)
float e1 = exp(a)
float f1 = (e1*(a*s1/(norm1*norm2)-c1/norm1)+ 1/norm1 )



float anew = (e1*(a*c1+norm2*s1)-a)/norm1
float bnew = b* f1
float cnew = c* f1
float dnew = d * f1




a = anew
b = bnew
c = cnew
d = dnew



if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  =  b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float  rho = a^2 +b^2+c^2+d^2

bailout:
rho <= r

default:
  title = "Quaternionic Phi 1Full - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-c-QuatPhi1_Jul"
  seed = #pixel
  plane = @plane  
  bailout = @bailout
  dim34 = @dim34
  fn1 = fn1
}

gwfa-c-QuatPhi1_Jul {
; full Quaternionic version of Phi1 - Julia variant
; bailout following the quaternion norm
init:

complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 b = imag(fn1(@dim34))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(@dim34))
 c = imag(fn1(@dim34))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 c = imag(fn1(@dim34))


endif



float r = real(@bailout)^2

float u = real(@seed)
float v = imag(@seed)

loop:


float norm1 = a^2+b^2+c^2+d^2
float norm2 = SQRT(b^2 + c^2 + d^2)
float s1 = SIN(norm2)
float c1 = COS(norm2)
float e1 = exp(a)
float f1 = (e1*(a*s1/(norm1*norm2)-c1/norm1)+ 1/norm1 )



float anew = (e1*(a*c1+norm2*s1)-a)/norm1 
float bnew = b* f1 
float cnew = c* f1 
float dnew = d * f1 




a = anew + u 
b = bnew + v 
c = cnew + u 
d = dnew + v

if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  = b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float rho = a^2 + b^2 + c^2 + d^2

bailout:
rho <= r
default:
  title = "Quaternionic Phi 1Full Jul - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 

  endparam
 



  param seed
    caption = "fixed perturbation"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-c-QuatPhi1"
    plane = @plane
  bailout = @bailout
  dim34 = @dim34
  fn1 = fn1

}


; added October 19,1999:  Reciprocal Quaternionic Phi:


gwfa-c-QuatPhi1Recip {
; full reciprocal Quaternionic version of Phi1 
; bailout following the quaternion norm
init:

complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 b = imag(fn1(@dim34))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(@dim34))
 c = imag(fn1(@dim34))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 c = imag(fn1(@dim34))


endif



float r = real(@bailout)^2
loop:



float norm1 = a^2+b^2+c^2+d^2
float norm2 = SQRT(b^2 + c^2 + d^2)
float s1 = SIN(norm2)
float c1 = COS(norm2)
float e1 = exp(a)
float g1= (e1^2 - 2*e1*c1+1)
float h1= (e1*(norm2*c1-a*s1)-norm2)/(norm2*g1)


 float anew = (e1*(a*c1+norm2*s1)-a)/g1
float bnew = b*h1
float cnew = c*h1
float dnew = d*h1



a = anew
b = bnew
c = cnew
d = dnew



if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  =  b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float  rho = a^2 +b^2+c^2+d^2

bailout:
rho <= r

default:
  title = "Quaternionic Reciprocal Phi 1Full - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.125
  
  param bailout
    caption = "Bailout Value"
    default = (10,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (0.1,0.1)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-c-QuatPhi1Recip_Jul"
  seed = #pixel
  plane = @plane  
  bailout = @bailout
  dim34 = @dim34
  fn1 = fn1
}

gwfa-c-QuatPhi1Recip_Jul {
; full Quaternionic version of Phi1Reciprocal - Julia variant
; bailout following the quaternion norm
init:

complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 b = imag(fn1(@dim34))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(@dim34))
 c = imag(fn1(@dim34))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 c = imag(fn1(@dim34))


endif



float r = real(@bailout)^2

float u = real(@seed)
float v = imag(@seed)

loop:


float norm1 = a^2+b^2+c^2+d^2
float norm2 = SQRT(b^2 + c^2 + d^2)
float s1 = SIN(norm2)
float c1 = COS(norm2)
float e1 = exp(a)
float g1= (e1^2 - 2*e1*c1+1)
float h1= (e1*(norm2*c1-a*s1)-norm2)/(norm2*g1)





 float anew = (e1*(a*c1+norm2*s1)-a)/g1
float bnew = b*h1
float cnew = c*h1
float dnew = d*h1





a = anew + u 
b = bnew + v 
c = cnew + u 
d = dnew + v

if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  = b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float rho = a^2 + b^2 + c^2 + d^2

bailout:
rho <= r
default:
  title = "Quaternionic Reciprocal Phi 1Full Jul - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.05
  
  param bailout
    caption = "Bailout Value"
    default = (10,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 

  endparam
 



  param seed
    caption = "fixed perturbation"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-c-QuatPhi1Recip"
    plane = @plane
  bailout = @bailout
  dim34 = @dim34
  fn1 = fn1

}





; added January 11, 2000:   Quaternionic Newton Algorithm for Phi = 1:




gwfa-c-QuatPhiNewt1 {
; full Quaternionic Newton Algorithm for the equation  Phi = 1 around zero 
; bailout following the quaternion norm
init:

complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 b = imag(fn1(@dim34))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(@dim34))
 c = imag(fn1(@dim34))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 c = imag(fn1(@dim34))


endif



float r = real(@bailout)^2
loop:




float norm1 = a^2+b^2+c^2+d^2
float  norm2quadrat = norm1 - a^2
float  norm2 = SQRT(norm2quadrat)
float  e1 = EXP(a)
float  s1 = SIN(norm2)
float  c1 = COS(norm2)

float  ausdr1 = a * (a^2-a + norm2quadrat) + norm2quadrat
float  ausdr2 = a^2 - 2*a + norm2quadrat
float  ausdr3 = e1 *  (ausdr1 * s1-ausdr2 * c1 * norm2 ) - 2 * a * norm2
float  ausdr4 = e1 * (a^3 * (c1 * norm2  - s1)+a^2 * norm2quadrat * s1+a * norm2quadrat * (c1 * norm2 -3 * s1) + norm2^3 * (2 * c1+ norm2 * s1)) -2 * norm2^3
float  ausdr5 = e1* (a^3 * (b^2* c1 *norm2 + s1 * (c^2+d^2)) +a^2 * norm2quadrat * (b^2 * s1-c1 * norm2quadrat )+a * norm2quadrat * (b^2 * (c1 * norm2quadrat -2 * s1) + s1 *  (c^2+d^2)) + norm2^3 *  (b^2 * (c1+ norm2 * s1) -c1 * (c^2+d^2)))
float  ausdr6 = e1* (ausdr1 * c1+ ausdr2 * norm2 * s1) + a^2-b^2-c^2-d^2
float  ausdr7 = e1*  (a^3 * (b^2 * s1+c^2 * c1 * norm2 +d^2 * s1) + a^2 * norm2quadrat * (c^2 * s1-c1 *norm2 )+a * norm2quadrat * (b^2 * s1+c^2 * ( c1 * norm2 - 2 * s1 ) + d^2 * s1 ) - norm2^3 * (b^2 * c1-c^2 *  (c1+ norm2 * s1) + c1 * d^2 )) + norm2^3  * (a^2+b^2-c^2+d^2)
float  ausdr8 = e1* (a^3 *  (b^2 * s1+c^2 * s1+c1 * d^2 * norm2) + a^2 * norm2quadrat * (d^2 * s1-c1 *  norm2 ) + a * norm2quadrat * (b^2 * s1+c^2 * s1 + d^2 * (c1 * norm2 - 2 * s1)) - norm2^3  * (b^2 * c1+c^2 * c1 - d^2 * (c1+ norm2 * s1))) + norm2^3 * (a^2+b^2+c^2-d^2)
float  ausdr9 = e1* (ausdr1 * s1 - ausdr2 * c1 *  norm2)
float  ausdr10 = norm2 * (a^2 * norm2^3 + ausdr5 - norm2^3 * (b^2-c^2-d^2))

float  ausdr11 = (ausdr10-ausdr4*b^2*norm2)*(ausdr4*d^2-ausdr8)
float  ausdr12 = (2*ausdr4*c^2*d^2-ausdr7*d^2-ausdr8*c^2)
float  ausdr13 = (ausdr7*d^2+ausdr8*c^2)
float  ausdr14 = (ausdr9*c1*d^2*norm2-ausdr8*s1)
float  ausdr15 = (ausdr10-ausdr4*b^2*norm2)*(ausdr4*c^2-ausdr7)
float  ausdr16 = (3*ausdr4^2*c^2*d^2-2*ausdr4*ausdr13+ausdr7*ausdr8)
float  ausdr17 = (ausdr9*s1+2*c1*norm2quadrat)
float  ausdr18 = (ausdr8*c1-3*d^2*norm2*ausdr17)
float  ausdr19 = ausdr4^2*(ausdr7*c1*d^2+c^2*ausdr18)
float  ausdr20 = (2*ausdr4^3*c^2*c1*d^2-ausdr19-2*ausdr4*norm2*ausdr13*ausdr17+ausdr7*ausdr8*norm2*ausdr17)
float  ausdr21 = (ausdr4^2*c^2*d^2-ausdr7*ausdr8)
float  ausdr22 = (ausdr4^2*c^2*d^2+4*ausdr4*c^2*d^2*norm2^3-ausdr7*(ausdr8+2*d^2*norm2^3)-2*ausdr8*c^2*norm2^3)
float  ausdr23 = (ausdr10*ausdr21-ausdr4^2*b^2*norm2*ausdr12)
float  ausdr24 = (2*ausdr4^3*c^2*d^2*norm2quadrat-ausdr4^2*(ausdr7*d^2*norm2quadrat+ausdr8*c^2*norm2quadrat+3*ausdr9*c^2*d^2*norm2^2)+2*ausdr4*ausdr9*norm2quadrat*ausdr13-ausdr7*ausdr8*ausdr9*norm2^2)
float  ausdr25 = (ausdr4^2*c^2*d^2-ausdr4*ausdr13+ausdr7*ausdr8)
float  ausdr26 = (ausdr8-ausdr4*d^2)*(ausdr4*c^2-ausdr7)
float  ausdr27 = (ausdr8-ausdr4*d^2)*(ausdr3*s1+ausdr6*c1)
float  ausdr28 = (ausdr10*(ausdr3*s1*(ausdr4*d^2-ausdr8)+ausdr4*ausdr6*c1*d^2-ausdr6*ausdr8*c1)+ausdr4*b^2*norm2*ausdr27)
float  ausdr29 = (ausdr10*(ausdr3*s1*(ausdr4*c^2-ausdr7)+ausdr4*ausdr6*c^2*c1-ausdr6*ausdr7*c1)+ausdr4*b^2*norm2*(ausdr7-ausdr4*c^2)*(ausdr3*s1+ausdr6*c1))

float  ausdr30 = ausdr10 * (2 * ausdr4 * c^2 * d^2- ausdr7 * d^2 - ausdr8 * c^2 ) - b^2 * norm2 * (3 * ausdr4^2 * c^2 * d^2 - 2 * ausdr4 * (ausdr7 * d^2 + ausdr8 * c^2 )+ausdr7 * ausdr8 )
float  ausdr31 = 2 * a * ausdr3 * norm2quadrat * ausdr30-ausdr10 * (ausdr3 * ausdr9 * norm2 * (2 * ausdr4 * c^2 * d^2-ausdr7 * d^2-ausdr8 * c^2)+ausdr6 * (ausdr4^2 * c^2 * d^2-ausdr7 * ausdr8))+b^2 * norm2 * (ausdr3 * ausdr9 * norm2 * (3 * ausdr4^2 * c^2 * d^2-2 * ausdr4 * (ausdr7 * d^2+ausdr8 * c^2)+ausdr7 * ausdr8)+ausdr4^2 * ausdr6 * (2 * ausdr4 * c^2 * d^2-ausdr7 * d^2-ausdr8 * c^2))

float ausdr32 = (ausdr10*ausdr12-b^2*norm2*ausdr16)
float ausdr33 = (ausdr4^2*c^2*c1*d^2+2*ausdr4*c^2*d^2*norm2*ausdr17-ausdr7*(ausdr8*c1+d^2*norm2*ausdr17)-ausdr8*c^2*norm2*ausdr17)
float ausdr34 = (ausdr4^2*c^2*d^2*s1-2*ausdr4*ausdr9*c^2*c1*d^2*norm2+ausdr7*ausdr14+ausdr8*ausdr9*c^2*c1*norm2)

float relax = real (@lambda)


float anew = (2*a^2*norm2quadrat*s1*ausdr32+a*(b^2*norm2*ausdr20-ausdr10*ausdr33)-norm2*(ausdr10*-b^2*norm2*(2*ausdr4^3*c^2*d^2*s1-ausdr4^2*(ausdr7*d^2*s1+c^2*(ausdr8*s1+3*ausdr9*c1*d^2*norm2))+2*ausdr4*ausdr9*c1*norm2*ausdr13-ausdr7*ausdr8*ausdr9*c1*norm2)))

anew = -relax*norm1*e1* anew

float anew = anew - (a^2*relax*norm1*ausdr23+a*(ausdr10*relax*norm1*ausdr22-2*ausdr4^3*b^2*c^2*d^2*relax*norm1*norm2+ausdr4^2*b^2*relax*norm1*norm2*(ausdr7*d^2+ausdr8*c^2-6*c^2*d^2*norm2^3)+4*ausdr4*b^2*relax*norm1*norm2quadrat^2*ausdr13-2*ausdr7*ausdr8*b^2*relax*norm1*norm2quadrat^2-ausdr31)+relax*norm1*(ausdr10*(ausdr4^2*c^2*d^2*(norm2quadrat)-2*ausdr4*ausdr9*c^2*d^2*norm2quadrat+ausdr7*(ausdr9*d^2*norm2quadrat-ausdr8*(norm2quadrat))+ausdr8*ausdr9*c^2*norm2quadrat)-b^2*norm2*(2*ausdr4^3*c^2*d^2*norm2quadrat-ausdr4^2*ausdr24)))

anew = anew / ausdr31

float bnew = b*relax*norm1*norm2^3*e1*(a*(ausdr3*c1*ausdr25+ausdr6*s1*ausdr26)+norm2*(ausdr3*s1*ausdr25+ausdr6*c1*(ausdr4*c^2-ausdr7)*(ausdr4*d^2-ausdr8)))/ausdr31-b*(a^2*ausdr3*relax*norm1*norm2^3*(ausdr4*c^2-ausdr7)*(ausdr4*d^2-ausdr8)+a*ausdr3*relax*norm1*norm2^3*ausdr25+ausdr3*relax*norm1*norm2^3*(norm2quadrat)*ausdr25+ausdr4^2*ausdr6*c^2*d^2*relax*norm1*norm2quadrat^2-ausdr4*ausdr6*relax*norm1*norm2quadrat^2*ausdr13+ausdr6*ausdr7*ausdr8*relax*norm1*norm2quadrat^2-ausdr31)/ausdr31

float cnew = c*(a^2*ausdr3*relax*norm1*norm2quadrat*ausdr11+a*ausdr3*relax*norm1*norm2quadrat*ausdr11+ausdr10*relax*norm1*norm2quadrat*(ausdr3*(b^2+c^2+d^2)+ausdr6*norm2)*(ausdr4*d^2-ausdr8)+ausdr3*ausdr4*b^2*relax*norm1*norm2^3*(norm2quadrat)*(ausdr8-ausdr4*d^2)-ausdr4^2*ausdr6*b^2*d^2*relax*norm1*norm2^4+ausdr4*ausdr6*ausdr8*b^2*relax*norm1*norm2quadrat^2+ausdr31)/ausdr31-c*relax*norm1*norm2^2*e1*(a*(ausdr10*(ausdr3*c1*(ausdr4*d^2-ausdr8)-ausdr4*ausdr6*d^2*s1+ausdr6*ausdr8*s1)+ausdr4*b^2*norm2*(ausdr8-ausdr4*d^2)*(ausdr3*c1-ausdr6*s1))+norm2*ausdr28)/ausdr31

float dnew = d*(a^2*ausdr3*relax*norm1*norm2quadrat*ausdr15+a*ausdr3*relax*norm1*norm2quadrat*ausdr15+ausdr10*relax*norm1*norm2quadrat*(ausdr3*(norm2quadrat)+ausdr6*norm2)*(ausdr4*c^2-ausdr7)+ausdr3*ausdr4*b^2*relax*norm1*norm2^3*(norm2quadrat)*(ausdr7-ausdr4*c^2)-ausdr4^2*ausdr6*b^2*c^2*relax*norm1*norm2quadrat^2+ausdr4*ausdr6*ausdr7*b^2*relax*norm1*norm2quadrat^2+ausdr31)/ausdr31-d*relax*norm1*norm2quadrat*e1*(a*(ausdr10*(ausdr3*c1*(ausdr4*c^2-ausdr7)-ausdr4*ausdr6*c^2*s1+ausdr6*ausdr7*s1)+ausdr4*b^2*norm2*(ausdr7-ausdr4*c^2)*(ausdr3*c1-ausdr6*s1))+norm2*ausdr29)/ausdr31




a = anew
b = bnew
c = cnew
d = dnew



if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  =  b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float  rho = a^2 +norm2quadrat

bailout:
rho <= r

default:
  title = "Quaternionic Newton 1 - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam

  param lambda
    caption = "Relaxation Factor"
    default = (1,0)
  endparam



  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-c-QuatPhiNewt1_Jul"
  seed = #pixel
  plane = @plane  
  lambda = @lambda
  bailout = @bailout
  dim34 = @dim34
  fn1 = fn1
}

gwfa-c-QuatPhiNewt1_Jul {
; full Quaternionic Newton Algorithm for the equation  Phi = 1 around zero - Julia Variant
; bailout following the quaternion norm
init:

complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 b = imag(fn1(@dim34))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(@dim34))
 c = imag(fn1(@dim34))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 c = imag(fn1(@dim34))


endif



float r = real(@bailout)^2

float u = real(fn2(@seed))
float v = imag(fn2(@seed))
float u1 = real(fn2(@seed2))
float v1 = imag(fn2(@seed2))



loop:





float norm1 = a^2+b^2+c^2+d^2
float  norm2quadrat = norm1 - a^2
float  norm2 = SQRT(norm2quadrat)
float  e1 = EXP(a)
float  s1 = SIN(norm2)
float  c1 = COS(norm2)

float  ausdr1 = a * (a^2-a + norm2quadrat) + norm2quadrat
float  ausdr2 = a^2 - 2*a + norm2quadrat
float  ausdr3 = e1 *  (ausdr1 * s1-ausdr2 * c1 * norm2 ) - 2 * a * norm2
float  ausdr4 = e1 * (a^3 * (c1 * norm2  - s1)+a^2 * norm2quadrat * s1+a * norm2quadrat * (c1 * norm2 -3 * s1) + norm2^3 * (2 * c1+ norm2 * s1)) -2 * norm2^3
float  ausdr5 = e1* (a^3 * (b^2* c1 *norm2 + s1 * (c^2+d^2)) +a^2 * norm2quadrat * (b^2 * s1-c1 * norm2quadrat )+a * norm2quadrat * (b^2 * (c1 * norm2quadrat -2 * s1) + s1 *  (c^2+d^2)) + norm2^3 *  (b^2 * (c1+ norm2 * s1) -c1 * (c^2+d^2)))
float  ausdr6 = e1* (ausdr1 * c1+ ausdr2 * norm2 * s1) + a^2-b^2-c^2-d^2
float  ausdr7 = e1*  (a^3 * (b^2 * s1+c^2 * c1 * norm2 +d^2 * s1) + a^2 * norm2quadrat * (c^2 * s1-c1 *norm2 )+a * norm2quadrat * (b^2 * s1+c^2 * ( c1 * norm2 - 2 * s1 ) + d^2 * s1 ) - norm2^3 * (b^2 * c1-c^2 *  (c1+ norm2 * s1) + c1 * d^2 )) + norm2^3  * (a^2+b^2-c^2+d^2)
float  ausdr8 = e1* (a^3 *  (b^2 * s1+c^2 * s1+c1 * d^2 * norm2) + a^2 * norm2quadrat * (d^2 * s1-c1 *  norm2 ) + a * norm2quadrat * (b^2 * s1+c^2 * s1 + d^2 * (c1 * norm2 - 2 * s1)) - norm2^3  * (b^2 * c1+c^2 * c1 - d^2 * (c1+ norm2 * s1))) + norm2^3 * (a^2+b^2+c^2-d^2)
float  ausdr9 = e1* (ausdr1 * s1 - ausdr2 * c1 *  norm2)
float  ausdr10 = norm2 * (a^2 * norm2^3 + ausdr5 - norm2^3 * (b^2-c^2-d^2))

float  ausdr11 = (ausdr10-ausdr4*b^2*norm2)*(ausdr4*d^2-ausdr8)
float  ausdr12 = (2*ausdr4*c^2*d^2-ausdr7*d^2-ausdr8*c^2)
float  ausdr13 = (ausdr7*d^2+ausdr8*c^2)
float  ausdr14 = (ausdr9*c1*d^2*norm2-ausdr8*s1)
float  ausdr15 = (ausdr10-ausdr4*b^2*norm2)*(ausdr4*c^2-ausdr7)
float  ausdr16 = (3*ausdr4^2*c^2*d^2-2*ausdr4*ausdr13+ausdr7*ausdr8)
float  ausdr17 = (ausdr9*s1+2*c1*norm2quadrat)
float  ausdr18 = (ausdr8*c1-3*d^2*norm2*ausdr17)
float  ausdr19 = ausdr4^2*(ausdr7*c1*d^2+c^2*ausdr18)
float  ausdr20 = (2*ausdr4^3*c^2*c1*d^2-ausdr19-2*ausdr4*norm2*ausdr13*ausdr17+ausdr7*ausdr8*norm2*ausdr17)
float  ausdr21 = (ausdr4^2*c^2*d^2-ausdr7*ausdr8)
float  ausdr22 = (ausdr4^2*c^2*d^2+4*ausdr4*c^2*d^2*norm2^3-ausdr7*(ausdr8+2*d^2*norm2^3)-2*ausdr8*c^2*norm2^3)
float  ausdr23 = (ausdr10*ausdr21-ausdr4^2*b^2*norm2*ausdr12)
float  ausdr24 = (2*ausdr4^3*c^2*d^2*norm2quadrat-ausdr4^2*(ausdr7*d^2*norm2quadrat+ausdr8*c^2*norm2quadrat+3*ausdr9*c^2*d^2*norm2^2)+2*ausdr4*ausdr9*norm2quadrat*ausdr13-ausdr7*ausdr8*ausdr9*norm2^2)
float  ausdr25 = (ausdr4^2*c^2*d^2-ausdr4*ausdr13+ausdr7*ausdr8)
float  ausdr26 = (ausdr8-ausdr4*d^2)*(ausdr4*c^2-ausdr7)
float  ausdr27 = (ausdr8-ausdr4*d^2)*(ausdr3*s1+ausdr6*c1)
float  ausdr28 = (ausdr10*(ausdr3*s1*(ausdr4*d^2-ausdr8)+ausdr4*ausdr6*c1*d^2-ausdr6*ausdr8*c1)+ausdr4*b^2*norm2*ausdr27)
float  ausdr29 = (ausdr10*(ausdr3*s1*(ausdr4*c^2-ausdr7)+ausdr4*ausdr6*c^2*c1-ausdr6*ausdr7*c1)+ausdr4*b^2*norm2*(ausdr7-ausdr4*c^2)*(ausdr3*s1+ausdr6*c1))

float  ausdr30 = ausdr10 * (2 * ausdr4 * c^2 * d^2- ausdr7 * d^2 - ausdr8 * c^2 ) - b^2 * norm2 * (3 * ausdr4^2 * c^2 * d^2 - 2 * ausdr4 * (ausdr7 * d^2 + ausdr8 * c^2 )+ausdr7 * ausdr8 )
float  ausdr31 = 2 * a * ausdr3 * norm2quadrat * ausdr30-ausdr10 * (ausdr3 * ausdr9 * norm2 * (2 * ausdr4 * c^2 * d^2-ausdr7 * d^2-ausdr8 * c^2)+ausdr6 * (ausdr4^2 * c^2 * d^2-ausdr7 * ausdr8))+b^2 * norm2 * (ausdr3 * ausdr9 * norm2 * (3 * ausdr4^2 * c^2 * d^2-2 * ausdr4 * (ausdr7 * d^2+ausdr8 * c^2)+ausdr7 * ausdr8)+ausdr4^2 * ausdr6 * (2 * ausdr4 * c^2 * d^2-ausdr7 * d^2-ausdr8 * c^2))

float ausdr32 = (ausdr10*ausdr12-b^2*norm2*ausdr16)
float ausdr33 = (ausdr4^2*c^2*c1*d^2+2*ausdr4*c^2*d^2*norm2*ausdr17-ausdr7*(ausdr8*c1+d^2*norm2*ausdr17)-ausdr8*c^2*norm2*ausdr17)
float ausdr34 = (ausdr4^2*c^2*d^2*s1-2*ausdr4*ausdr9*c^2*c1*d^2*norm2+ausdr7*ausdr14+ausdr8*ausdr9*c^2*c1*norm2)

float relax = real (@lambda)


float anew = (2*a^2*norm2quadrat*s1*ausdr32+a*(b^2*norm2*ausdr20-ausdr10*ausdr33)-norm2*(ausdr10*-b^2*norm2*(2*ausdr4^3*c^2*d^2*s1-ausdr4^2*(ausdr7*d^2*s1+c^2*(ausdr8*s1+3*ausdr9*c1*d^2*norm2))+2*ausdr4*ausdr9*c1*norm2*ausdr13-ausdr7*ausdr8*ausdr9*c1*norm2)))

anew = -relax*norm1*e1* anew

float anew = anew - (a^2*relax*norm1*ausdr23+a*(ausdr10*relax*norm1*ausdr22-2*ausdr4^3*b^2*c^2*d^2*relax*norm1*norm2+ausdr4^2*b^2*relax*norm1*norm2*(ausdr7*d^2+ausdr8*c^2-6*c^2*d^2*norm2^3)+4*ausdr4*b^2*relax*norm1*norm2quadrat^2*ausdr13-2*ausdr7*ausdr8*b^2*relax*norm1*norm2quadrat^2-ausdr31)+relax*norm1*(ausdr10*(ausdr4^2*c^2*d^2*(norm2quadrat)-2*ausdr4*ausdr9*c^2*d^2*norm2quadrat+ausdr7*(ausdr9*d^2*norm2quadrat-ausdr8*(norm2quadrat))+ausdr8*ausdr9*c^2*norm2quadrat)-b^2*norm2*(2*ausdr4^3*c^2*d^2*norm2quadrat-ausdr4^2*ausdr24)))

anew = anew / ausdr31

float bnew = b*relax*norm1*norm2^3*e1*(a*(ausdr3*c1*ausdr25+ausdr6*s1*ausdr26)+norm2*(ausdr3*s1*ausdr25+ausdr6*c1*(ausdr4*c^2-ausdr7)*(ausdr4*d^2-ausdr8)))/ausdr31-b*(a^2*ausdr3*relax*norm1*norm2^3*(ausdr4*c^2-ausdr7)*(ausdr4*d^2-ausdr8)+a*ausdr3*relax*norm1*norm2^3*ausdr25+ausdr3*relax*norm1*norm2^3*(norm2quadrat)*ausdr25+ausdr4^2*ausdr6*c^2*d^2*relax*norm1*norm2quadrat^2-ausdr4*ausdr6*relax*norm1*norm2quadrat^2*ausdr13+ausdr6*ausdr7*ausdr8*relax*norm1*norm2quadrat^2-ausdr31)/ausdr31

float cnew = c*(a^2*ausdr3*relax*norm1*norm2quadrat*ausdr11+a*ausdr3*relax*norm1*norm2quadrat*ausdr11+ausdr10*relax*norm1*norm2quadrat*(ausdr3*(b^2+c^2+d^2)+ausdr6*norm2)*(ausdr4*d^2-ausdr8)+ausdr3*ausdr4*b^2*relax*norm1*norm2^3*(norm2quadrat)*(ausdr8-ausdr4*d^2)-ausdr4^2*ausdr6*b^2*d^2*relax*norm1*norm2^4+ausdr4*ausdr6*ausdr8*b^2*relax*norm1*norm2quadrat^2+ausdr31)/ausdr31-c*relax*norm1*norm2^2*e1*(a*(ausdr10*(ausdr3*c1*(ausdr4*d^2-ausdr8)-ausdr4*ausdr6*d^2*s1+ausdr6*ausdr8*s1)+ausdr4*b^2*norm2*(ausdr8-ausdr4*d^2)*(ausdr3*c1-ausdr6*s1))+norm2*ausdr28)/ausdr31

float dnew = d*(a^2*ausdr3*relax*norm1*norm2quadrat*ausdr15+a*ausdr3*relax*norm1*norm2quadrat*ausdr15+ausdr10*relax*norm1*norm2quadrat*(ausdr3*(norm2quadrat)+ausdr6*norm2)*(ausdr4*c^2-ausdr7)+ausdr3*ausdr4*b^2*relax*norm1*norm2^3*(norm2quadrat)*(ausdr7-ausdr4*c^2)-ausdr4^2*ausdr6*b^2*c^2*relax*norm1*norm2quadrat^2+ausdr4*ausdr6*ausdr7*b^2*relax*norm1*norm2quadrat^2+ausdr31)/ausdr31-d*relax*norm1*norm2quadrat*e1*(a*(ausdr10*(ausdr3*c1*(ausdr4*c^2-ausdr7)-ausdr4*ausdr6*c^2*s1+ausdr6*ausdr7*s1)+ausdr4*b^2*norm2*(ausdr7-ausdr4*c^2)*(ausdr3*c1-ausdr6*s1))+norm2*ausdr29)/ausdr31







a = anew + u 
b = bnew + v 
c = cnew + u1 
d = dnew + v1

if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  = b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float rho = a^2 + b^2 + c^2 + d^2

bailout:
rho <= r
default:
  title = "Quaternionic Newton Julia 1 - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 

  endparam
 
  param lambda
    caption = "Relaxation Factor"
    default = (1,0)
  endparam




  param seed
    caption = "fixed perturbation"
    default = (1,0)
  endparam

  param seed2
    caption = "2nd fixed perturbation"
    default = (0,0)
  endparam



 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc

 func fn2
    caption = "Seed function"
    default = ident()
  endfunc



switch:
  type = "gwfa-c-QuatPhiNewt1"
    plane = @plane
  bailout = @bailout
  lambda = @lambda
  dim34 = @dim34
  fn1 = fn1

}





gwfa-c-QuatPhiNewt2 {
; full Quaternionic Newton Algorithm for the equation  Phi = 1 around zero 
; bailout following the quaternion norm
init:

complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 b = imag(fn1(@dim34))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(@dim34))
 c = imag(fn1(@dim34))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 c = imag(fn1(@dim34))


endif



float r = real(@bailout)^2
float relax = real (@lambda)

loop:



float norm1 = a^2+b^2+c^2+d^2
float  norm2quadrat = norm1 - a^2
float  norm2 = SQRT(norm2quadrat)
float  e1 = EXP(a)
float  s1 = SIN(norm2)
float  c1 = COS(norm2)

float  ausdr1 = a * (a^2-a + norm2quadrat) + norm2quadrat
float  ausdr2 = a^2 - 2*a + norm2quadrat
float  ausdr3 = e1 *  (ausdr1 * s1-ausdr2 * c1 * norm2 ) - 2 * a * norm2
float  ausdr4 = e1 * (a^3 * (c1 * norm2  - s1)+a^2 * norm2quadrat * s1+a * norm2quadrat * (c1 * norm2 -3 * s1) + norm2^3 * (2 * c1+ norm2 * s1)) -2 * norm2^3
float  ausdr5 = e1* (a^3 * (b^2* c1 *norm2 + s1 * (c^2+d^2)) +a^2 * norm2quadrat * (b^2 * s1-c1 * norm2quadrat )+a * norm2quadrat * (b^2 * (c1 * norm2quadrat -2 * s1) + s1 *  (c^2+d^2)) + norm2^3 *  (b^2 * (c1+ norm2 * s1) -c1 * (c^2+d^2)))
float  ausdr6 = e1* (ausdr1 * c1+ ausdr2 * norm2 * s1) + a^2-b^2-c^2-d^2
float  ausdr7 = e1*  (a^3 * (b^2 * s1+c^2 * c1 * norm2 +d^2 * s1) + a^2 * norm2quadrat * (c^2 * s1-c1 *norm2 )+a * norm2quadrat * (b^2 * s1+c^2 * ( c1 * norm2 - 2 * s1 ) + d^2 * s1 ) - norm2^3 * (b^2 * c1-c^2 *  (c1+ norm2 * s1) + c1 * d^2 )) + norm2^3  * (a^2+b^2-c^2+d^2)
float  ausdr8 = e1* (a^3 *  (b^2 * s1+c^2 * s1+c1 * d^2 * norm2) + a^2 * norm2quadrat * (d^2 * s1-c1 *  norm2 ) + a * norm2quadrat * (b^2 * s1+c^2 * s1 + d^2 * (c1 * norm2 - 2 * s1)) - norm2^3  * (b^2 * c1+c^2 * c1 - d^2 * (c1+ norm2 * s1))) + norm2^3 * (a^2+b^2+c^2-d^2)
float  ausdr9 = e1* (ausdr1 * s1 - ausdr2 * c1 *  norm2)
float  ausdr10 = norm2 * (a^2 * norm2^3 + ausdr5 - norm2^3 * (b^2-c^2-d^2))

float  ausdr11 = (ausdr10-ausdr4*b^2*norm2)*(ausdr4*d^2-ausdr8)
float  ausdr12 = (2*ausdr4*c^2*d^2-ausdr7*d^2-ausdr8*c^2)
float  ausdr13 = (ausdr7*d^2+ausdr8*c^2)
float  ausdr14 = (ausdr9*c1*d^2*norm2-ausdr8*s1)
float  ausdr15 = (ausdr10-ausdr4*b^2*norm2)*(ausdr4*c^2-ausdr7)
float  ausdr16 = (3*ausdr4^2*c^2*d^2-2*ausdr4*ausdr13+ausdr7*ausdr8)

float  ausdr21 = (ausdr4^2*c^2*d^2-ausdr7*ausdr8)
float  ausdr22 = (ausdr4^2*c^2*d^2+4*ausdr4*c^2*d^2*norm2^3-ausdr7*(ausdr8+2*d^2*norm2^3)-2*ausdr8*c^2*norm2^3)
float  ausdr23 = (ausdr10*ausdr21-ausdr4^2*b^2*norm2*ausdr12)

float  ausdr25 = (ausdr4^2*c^2*d^2-ausdr4*ausdr13+ausdr7*ausdr8)
float  ausdr26 = (ausdr8-ausdr4*d^2)*(ausdr4*c^2-ausdr7)
float  ausdr27 = (ausdr8-ausdr4*d^2)*(ausdr3*s1+ausdr6*c1)
float  ausdr28 = (ausdr10*(ausdr3*s1*(ausdr4*d^2-ausdr8)+ausdr4*ausdr6*c1*d^2-ausdr6*ausdr8*c1)+ausdr4*b^2*norm2*ausdr27)



float  ausdr30 = ausdr10 * (2 * ausdr4 * c^2 * d^2- ausdr7 * d^2 - ausdr8 * c^2 ) - b^2 * norm2 * (3 * ausdr4^2 * c^2 * d^2 - 2 * ausdr4 * (ausdr7 * d^2 + ausdr8 * c^2 )+ausdr7 * ausdr8 )
float  ausdr31 = 2 * a * ausdr3 * norm2quadrat * ausdr30-ausdr10 * (ausdr3 * ausdr9 * norm2 * (2 * ausdr4 * c^2 * d^2-ausdr7 * d^2-ausdr8 * c^2)+ausdr6 * (ausdr4^2 * c^2 * d^2-ausdr7 * ausdr8))+b^2 * norm2 * (ausdr3 * ausdr9 * norm2 * (3 * ausdr4^2 * c^2 * d^2-2 * ausdr4 * (ausdr7 * d^2+ausdr8 * c^2)+ausdr7 * ausdr8)+ausdr4^2 * ausdr6 * (2 * ausdr4 * c^2 * d^2-ausdr7 * d^2-ausdr8 * c^2))

float ausdr32 = (ausdr10*ausdr12-b^2*norm2*ausdr16)

float ausdr34 = (ausdr4^2*c^2*d^2*s1-2*ausdr4*ausdr9*c^2*c1*d^2*norm2+ausdr7*ausdr14+ausdr8*ausdr9*c^2*c1*norm2)


float ausdr35 = (2*ausdr4^3*c^2*d^2*(b^2+c^2+d^2)-ausdr4^2*(ausdr7*d^2*(b^2+c^2+d^2)+ausdr8*c^2*(b^2+c^2+d^2)+3*ausdr9*c^2*d^2*norm2quadrat)+2*ausdr4*ausdr9*norm2quadrat*ausdr13-ausdr7*ausdr8*ausdr9*norm2quadrat)

float ausdr36 = (ausdr10*(ausdr4^2*c^2*d^2*(b^2+c^2+d^2)-2*ausdr4*ausdr9*c^2*d^2*norm2quadrat+ausdr7*(ausdr9*d^2*norm2quadrat-ausdr8*(b^2+c^2+d^2))+ausdr8*ausdr9*c^2*norm2quadrat)-b^2*norm2*ausdr35)


float ausdr37 = (ausdr10*relax*norm1*ausdr22-2*ausdr4^3*b^2*c^2*d^2*relax*norm1*norm2+ausdr4^2*b^2*relax*norm1*norm2*(ausdr7*d^2+ausdr8*c^2-6*c^2*d^2*norm2^3)+4*ausdr4*b^2*relax*norm1*norm2^4*ausdr13-2*ausdr7*ausdr8*b^2*relax*norm1*norm2^4-ausdr31)

float ausdr38 = (2*ausdr4^3*c^2*d^2*s1-ausdr4^2*(ausdr7*d^2*s1+c^2*(ausdr8*s1+3*ausdr9*c1*d^2*norm2))+2*ausdr4*ausdr9*c1*norm2*ausdr13-ausdr7*ausdr8*ausdr9*c1*norm2)

float ausdr39 = (ausdr4^2*c^2*c1*d^2+2*ausdr4*c^2*d^2*norm2*(ausdr9*s1+2*c1*norm2quadrat)-ausdr7*(ausdr8*c1+d^2*norm2*(ausdr9*s1+2*c1*norm2quadrat))-ausdr8*c^2*norm2*(ausdr9*s1+2*c1*norm2quadrat))

float ausdr40 = ausdr4*b^2*norm2*(ausdr7-ausdr4*c^2)*(ausdr3*s1+ausdr6*c1)

float ausdr41 = norm2*(ausdr10*(ausdr3*s1*(ausdr4*c^2-ausdr7)+ausdr4*ausdr6*c^2*c1-ausdr6*ausdr7*c1)+ausdr40)

float anew = -relax*norm1*e1*(2*a^2*norm2quadrat*s1*ausdr32+a*(b^2*norm2*(2*ausdr4^3*c^2*c1*d^2-ausdr4^2*(ausdr7*c1*d^2+c^2*(ausdr8*c1-3*d^2*norm2*(ausdr9*s1+2*c1*norm2quadrat)))-2*ausdr4*norm2*ausdr13*(ausdr9*s1+2*c1*norm2quadrat)+ausdr7*ausdr8*norm2*(ausdr9*s1+2*c1*norm2quadrat))-ausdr10*ausdr39)-norm2*(ausdr10*ausdr34-b^2*norm2*ausdr38))/ausdr31

float anew = anew - (a^2*relax*norm1*ausdr23+a*ausdr37+relax*norm1*ausdr36)/ausdr31



float bnew = b*relax*norm1*norm2^3*e1*(a*(ausdr3*c1*ausdr25+ausdr6*s1*ausdr26)+norm2*(ausdr3*s1*ausdr25+ausdr6*c1*(ausdr4*c^2-ausdr7)*(ausdr4*d^2-ausdr8)))/ausdr31-b*(a^2*ausdr3*relax*norm1*norm2^3*(ausdr4*c^2-ausdr7)*(ausdr4*d^2-ausdr8)+a*ausdr3*relax*norm1*norm2^3*ausdr25+ausdr3*relax*norm1*norm2^3*(b^2+c^2+d^2)*ausdr25+ausdr4^2*ausdr6*c^2*d^2*relax*norm1*norm2^4-ausdr4*ausdr6*relax*norm1*norm2^4*ausdr13+ausdr6*ausdr7*ausdr8*relax*norm1*norm2^4-ausdr31)/ausdr31


float cnew = c*(a^2*ausdr3*relax*norm1*norm2quadrat*ausdr11+a*ausdr3*relax*norm1*norm2quadrat*ausdr11+ausdr10*relax*norm1*norm2quadrat*(ausdr3*(b^2+c^2+d^2)+ausdr6*norm2)*(ausdr4*d^2-ausdr8)+ausdr3*ausdr4*b^2*relax*norm1*norm2^3*(b^2+c^2+d^2)*(ausdr8-ausdr4*d^2)-ausdr4^2*ausdr6*b^2*d^2*relax*norm1*norm2^4+ausdr4*ausdr6*ausdr8*b^2*relax*norm1*norm2^4+ausdr31)/ausdr31-c*relax*norm1*norm2quadrat*e1*(a*(ausdr10*(ausdr3*c1*(ausdr4*d^2-ausdr8)-ausdr4*ausdr6*d^2*s1+ausdr6*ausdr8*s1)+ausdr4*b^2*norm2*(ausdr8-ausdr4*d^2)*(ausdr3*c1-ausdr6*s1))+norm2*ausdr28)/ausdr31


float dnew = d*(a^2*ausdr3*relax*norm1*norm2quadrat*ausdr15+a*ausdr3*relax*norm1*norm2quadrat*ausdr15+ausdr10*relax*norm1*norm2quadrat*(ausdr3*(b^2+c^2+d^2)+ausdr6*norm2)*(ausdr4*c^2-ausdr7)+ausdr3*ausdr4*b^2*relax*norm1*norm2^3*(b^2+c^2+d^2)*(ausdr7-ausdr4*c^2)-ausdr4^2*ausdr6*b^2*c^2*relax*norm1*norm2^4+ausdr4*ausdr6*ausdr7*b^2*relax*norm1*norm2^4+ausdr31)/ausdr31-d*relax*norm1*norm2quadrat*e1*(a*(ausdr10*(ausdr3*c1*(ausdr4*c^2-ausdr7)-ausdr4*ausdr6*c^2*s1+ausdr6*ausdr7*s1)+ausdr4*b^2*norm2*(ausdr7-ausdr4*c^2)*(ausdr3*c1-ausdr6*s1))+ausdr41)/ausdr31


a = anew
b = bnew
c = cnew
d = dnew



if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  =  b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float  rho = a^2 +norm2quadrat

bailout:
rho <= r

default:
  title = "Quaternionic Newton 2 - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam

  param lambda
    caption = "Relaxation Factor"
    default = (1,0)
  endparam



  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-c-QuatPhiNewt2_Jul"
  seed = #pixel
  plane = @plane  
  lambda = @lambda
  bailout = @bailout
  dim34 = @dim34
  fn1 = fn1
}

gwfa-c-QuatPhiNewt2_Jul {
; full Quaternionic Newton Algorithm for the equation  Phi = 1 around zero - Julia Variant
; bailout following the quaternion norm
init:

complex  z = pixel 

if @plane == 0
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))

elseif @plane == 1
 b = real(z)
 c = imag(z)
 a = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 2
 c = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 b = imag(fn1(@dim34))

elseif @plane == 3
 d = real(z)
 a = imag(z)
 b = real(fn1(@dim34))
 c = imag(fn1(@dim34))

elseif @plane == 4
 a = real(z)
 c = imag(z)
 b = real(fn1(@dim34))
 d = imag(fn1(@dim34))

elseif @plane == 5
 b = real(z)
 d = imag(z)
 a = real(fn1(@dim34))
 c = imag(fn1(@dim34))


endif



float r = real(@bailout)^2

complex jul1 = @PixelWeight*fn2(@seed)
complex jul2 = @PixelWeight*fn2(@seed2)


float u = real(jul1)
float v = imag(jul1)
float u1 = real(jul2)
float v1 = imag(jul2)

float relax = real (@lambda)

loop:



float norm1 = a^2+b^2+c^2+d^2
float  norm2quadrat = norm1 - a^2
float  norm2 = SQRT(norm2quadrat)
float  e1 = EXP(a)
float  s1 = SIN(norm2)
float  c1 = COS(norm2)

float  ausdr1 = a * (a^2-a + norm2quadrat) + norm2quadrat
float  ausdr2 = a^2 - 2*a + norm2quadrat
float  ausdr3 = e1 *  (ausdr1 * s1-ausdr2 * c1 * norm2 ) - 2 * a * norm2
float  ausdr4 = e1 * (a^3 * (c1 * norm2  - s1)+a^2 * norm2quadrat * s1+a * norm2quadrat * (c1 * norm2 -3 * s1) + norm2^3 * (2 * c1+ norm2 * s1)) -2 * norm2^3
float  ausdr5 = e1* (a^3 * (b^2* c1 *norm2 + s1 * (c^2+d^2)) +a^2 * norm2quadrat * (b^2 * s1-c1 * norm2quadrat )+a * norm2quadrat * (b^2 * (c1 * norm2quadrat -2 * s1) + s1 *  (c^2+d^2)) + norm2^3 *  (b^2 * (c1+ norm2 * s1) -c1 * (c^2+d^2)))
float  ausdr6 = e1* (ausdr1 * c1+ ausdr2 * norm2 * s1) + a^2-b^2-c^2-d^2
float  ausdr7 = e1*  (a^3 * (b^2 * s1+c^2 * c1 * norm2 +d^2 * s1) + a^2 * norm2quadrat * (c^2 * s1-c1 *norm2 )+a * norm2quadrat * (b^2 * s1+c^2 * ( c1 * norm2 - 2 * s1 ) + d^2 * s1 ) - norm2^3 * (b^2 * c1-c^2 *  (c1+ norm2 * s1) + c1 * d^2 )) + norm2^3  * (a^2+b^2-c^2+d^2)
float  ausdr8 = e1* (a^3 *  (b^2 * s1+c^2 * s1+c1 * d^2 * norm2) + a^2 * norm2quadrat * (d^2 * s1-c1 *  norm2 ) + a * norm2quadrat * (b^2 * s1+c^2 * s1 + d^2 * (c1 * norm2 - 2 * s1)) - norm2^3  * (b^2 * c1+c^2 * c1 - d^2 * (c1+ norm2 * s1))) + norm2^3 * (a^2+b^2+c^2-d^2)
float  ausdr9 = e1* (ausdr1 * s1 - ausdr2 * c1 *  norm2)
float  ausdr10 = norm2 * (a^2 * norm2^3 + ausdr5 - norm2^3 * (b^2-c^2-d^2))

float  ausdr11 = (ausdr10-ausdr4*b^2*norm2)*(ausdr4*d^2-ausdr8)
float  ausdr12 = (2*ausdr4*c^2*d^2-ausdr7*d^2-ausdr8*c^2)
float  ausdr13 = (ausdr7*d^2+ausdr8*c^2)
float  ausdr14 = (ausdr9*c1*d^2*norm2-ausdr8*s1)
float  ausdr15 = (ausdr10-ausdr4*b^2*norm2)*(ausdr4*c^2-ausdr7)
float  ausdr16 = (3*ausdr4^2*c^2*d^2-2*ausdr4*ausdr13+ausdr7*ausdr8)

float  ausdr21 = (ausdr4^2*c^2*d^2-ausdr7*ausdr8)
float  ausdr22 = (ausdr4^2*c^2*d^2+4*ausdr4*c^2*d^2*norm2^3-ausdr7*(ausdr8+2*d^2*norm2^3)-2*ausdr8*c^2*norm2^3)
float  ausdr23 = (ausdr10*ausdr21-ausdr4^2*b^2*norm2*ausdr12)

float  ausdr25 = (ausdr4^2*c^2*d^2-ausdr4*ausdr13+ausdr7*ausdr8)
float  ausdr26 = (ausdr8-ausdr4*d^2)*(ausdr4*c^2-ausdr7)
float  ausdr27 = (ausdr8-ausdr4*d^2)*(ausdr3*s1+ausdr6*c1)
float  ausdr28 = (ausdr10*(ausdr3*s1*(ausdr4*d^2-ausdr8)+ausdr4*ausdr6*c1*d^2-ausdr6*ausdr8*c1)+ausdr4*b^2*norm2*ausdr27)



float  ausdr30 = ausdr10 * (2 * ausdr4 * c^2 * d^2- ausdr7 * d^2 - ausdr8 * c^2 ) - b^2 * norm2 * (3 * ausdr4^2 * c^2 * d^2 - 2 * ausdr4 * (ausdr7 * d^2 + ausdr8 * c^2 )+ausdr7 * ausdr8 )
float  ausdr31 = 2 * a * ausdr3 * norm2quadrat * ausdr30-ausdr10 * (ausdr3 * ausdr9 * norm2 * (2 * ausdr4 * c^2 * d^2-ausdr7 * d^2-ausdr8 * c^2)+ausdr6 * (ausdr4^2 * c^2 * d^2-ausdr7 * ausdr8))+b^2 * norm2 * (ausdr3 * ausdr9 * norm2 * (3 * ausdr4^2 * c^2 * d^2-2 * ausdr4 * (ausdr7 * d^2+ausdr8 * c^2)+ausdr7 * ausdr8)+ausdr4^2 * ausdr6 * (2 * ausdr4 * c^2 * d^2-ausdr7 * d^2-ausdr8 * c^2))

float ausdr32 = (ausdr10*ausdr12-b^2*norm2*ausdr16)

float ausdr34 = (ausdr4^2*c^2*d^2*s1-2*ausdr4*ausdr9*c^2*c1*d^2*norm2+ausdr7*ausdr14+ausdr8*ausdr9*c^2*c1*norm2)


float ausdr35 = (2*ausdr4^3*c^2*d^2*(b^2+c^2+d^2)-ausdr4^2*(ausdr7*d^2*(b^2+c^2+d^2)+ausdr8*c^2*(b^2+c^2+d^2)+3*ausdr9*c^2*d^2*norm2quadrat)+2*ausdr4*ausdr9*norm2quadrat*ausdr13-ausdr7*ausdr8*ausdr9*norm2quadrat)

float ausdr36 = (ausdr10*(ausdr4^2*c^2*d^2*(b^2+c^2+d^2)-2*ausdr4*ausdr9*c^2*d^2*norm2quadrat+ausdr7*(ausdr9*d^2*norm2quadrat-ausdr8*(b^2+c^2+d^2))+ausdr8*ausdr9*c^2*norm2quadrat)-b^2*norm2*ausdr35)


float ausdr37 = (ausdr10*relax*norm1*ausdr22-2*ausdr4^3*b^2*c^2*d^2*relax*norm1*norm2+ausdr4^2*b^2*relax*norm1*norm2*(ausdr7*d^2+ausdr8*c^2-6*c^2*d^2*norm2^3)+4*ausdr4*b^2*relax*norm1*norm2^4*ausdr13-2*ausdr7*ausdr8*b^2*relax*norm1*norm2^4-ausdr31)

float ausdr38 = (2*ausdr4^3*c^2*d^2*s1-ausdr4^2*(ausdr7*d^2*s1+c^2*(ausdr8*s1+3*ausdr9*c1*d^2*norm2))+2*ausdr4*ausdr9*c1*norm2*ausdr13-ausdr7*ausdr8*ausdr9*c1*norm2)

float ausdr39 = (ausdr4^2*c^2*c1*d^2+2*ausdr4*c^2*d^2*norm2*(ausdr9*s1+2*c1*norm2quadrat)-ausdr7*(ausdr8*c1+d^2*norm2*(ausdr9*s1+2*c1*norm2quadrat))-ausdr8*c^2*norm2*(ausdr9*s1+2*c1*norm2quadrat))

float ausdr40 = ausdr4*b^2*norm2*(ausdr7-ausdr4*c^2)*(ausdr3*s1+ausdr6*c1)

float ausdr41 = norm2*(ausdr10*(ausdr3*s1*(ausdr4*c^2-ausdr7)+ausdr4*ausdr6*c^2*c1-ausdr6*ausdr7*c1)+ausdr40)

float anew = -relax*norm1*e1*(2*a^2*norm2quadrat*s1*ausdr32+a*(b^2*norm2*(2*ausdr4^3*c^2*c1*d^2-ausdr4^2*(ausdr7*c1*d^2+c^2*(ausdr8*c1-3*d^2*norm2*(ausdr9*s1+2*c1*norm2quadrat)))-2*ausdr4*norm2*ausdr13*(ausdr9*s1+2*c1*norm2quadrat)+ausdr7*ausdr8*norm2*(ausdr9*s1+2*c1*norm2quadrat))-ausdr10*ausdr39)-norm2*(ausdr10*ausdr34-b^2*norm2*ausdr38))/ausdr31

float anew = anew - (a^2*relax*norm1*ausdr23+a*ausdr37+relax*norm1*ausdr36)/ausdr31



float bnew = b*relax*norm1*norm2^3*e1*(a*(ausdr3*c1*ausdr25+ausdr6*s1*ausdr26)+norm2*(ausdr3*s1*ausdr25+ausdr6*c1*(ausdr4*c^2-ausdr7)*(ausdr4*d^2-ausdr8)))/ausdr31-b*(a^2*ausdr3*relax*norm1*norm2^3*(ausdr4*c^2-ausdr7)*(ausdr4*d^2-ausdr8)+a*ausdr3*relax*norm1*norm2^3*ausdr25+ausdr3*relax*norm1*norm2^3*(b^2+c^2+d^2)*ausdr25+ausdr4^2*ausdr6*c^2*d^2*relax*norm1*norm2^4-ausdr4*ausdr6*relax*norm1*norm2^4*ausdr13+ausdr6*ausdr7*ausdr8*relax*norm1*norm2^4-ausdr31)/ausdr31


float cnew = c*(a^2*ausdr3*relax*norm1*norm2quadrat*ausdr11+a*ausdr3*relax*norm1*norm2quadrat*ausdr11+ausdr10*relax*norm1*norm2quadrat*(ausdr3*(b^2+c^2+d^2)+ausdr6*norm2)*(ausdr4*d^2-ausdr8)+ausdr3*ausdr4*b^2*relax*norm1*norm2^3*(b^2+c^2+d^2)*(ausdr8-ausdr4*d^2)-ausdr4^2*ausdr6*b^2*d^2*relax*norm1*norm2^4+ausdr4*ausdr6*ausdr8*b^2*relax*norm1*norm2^4+ausdr31)/ausdr31-c*relax*norm1*norm2quadrat*e1*(a*(ausdr10*(ausdr3*c1*(ausdr4*d^2-ausdr8)-ausdr4*ausdr6*d^2*s1+ausdr6*ausdr8*s1)+ausdr4*b^2*norm2*(ausdr8-ausdr4*d^2)*(ausdr3*c1-ausdr6*s1))+norm2*ausdr28)/ausdr31


float dnew = d*(a^2*ausdr3*relax*norm1*norm2quadrat*ausdr15+a*ausdr3*relax*norm1*norm2quadrat*ausdr15+ausdr10*relax*norm1*norm2quadrat*(ausdr3*(b^2+c^2+d^2)+ausdr6*norm2)*(ausdr4*c^2-ausdr7)+ausdr3*ausdr4*b^2*relax*norm1*norm2^3*(b^2+c^2+d^2)*(ausdr7-ausdr4*c^2)-ausdr4^2*ausdr6*b^2*c^2*relax*norm1*norm2^4+ausdr4*ausdr6*ausdr7*b^2*relax*norm1*norm2^4+ausdr31)/ausdr31-d*relax*norm1*norm2quadrat*e1*(a*(ausdr10*(ausdr3*c1*(ausdr4*c^2-ausdr7)-ausdr4*ausdr6*c^2*s1+ausdr6*ausdr7*s1)+ausdr4*b^2*norm2*(ausdr7-ausdr4*c^2)*(ausdr3*c1-ausdr6*s1))+ausdr41)/ausdr31






a = anew + u 
b = bnew + v 
c = cnew + u1 
d = dnew + v1

if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  = b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float rho = a^2 + b^2 + c^2 + d^2

bailout:
rho <= r
default:
  title = "Quaternionic Newton Julia 2 - c"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 

  endparam
 
  param lambda
    caption = "Relaxation Factor"
    default = (1,0)
  endparam




  param seed
    caption = "fixed perturbation"
    default = (1,0)
  endparam

  param seed2
    caption = "2nd fixed perturbation"
    default = (0,0)
  endparam


  param PixelWeight
    caption = "Weighing Factor"
    default = (1,0)
  endparam





 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc

 func fn2
    caption = "Seed function"
    default = ident()
  endfunc



switch:
  type = "gwfa-c-QuatPhiNewt2"
    plane = @plane
  bailout = @bailout
  lambda = @lambda
  dim34 = @dim34
  fn1 = fn1

}



; added on February 10,2000: The iteration of phi-1 in the complex algebra of 2x2-matrices:




gwfa-c-MatrixPhi_1 {
; Iteration of phi-1 in the algebra of complex 2x2 matrices
; bailout following the matrix norm Tr(aa')
; matrix notation like [[a,b],[c,d]]
init:

complex  z = pixel 

if @plane == 0
complex a = z
complex b = @bp1
complex c = @bp2
complex d = @bp3


elseif @plane == 1
complex b = z
complex c = @bp1
complex d = @bp2
complex a = @bp3


elseif @plane == 2
complex c = z
complex d = @bp1
complex a= @bp2
complex b = @bp3


elseif @plane == 3
complex d = z
complex a = @bp1
complex b = @bp2
complex c = @bp3



endif



float r = real(@bailout)^2

loop:

float  rho = cabs(a)^2+cabs(b)^2+cabs(c)^2+cabs(d)^2

if rho == 0

complex anew = 0
complex bnew =  0
complex cnew =  0
complex dnew = 0


else

complex bb1 = -SQRT(a^2-2*a*d+4*b*c+d^2)

complex bb2 = (SQRT(a^2-2*a*d+4*b*c+d^2)-a+d)

complex bb3 = (SQRT(a^2-2*a*d+4*b*c+d^2)/2+a/2+d/2)

complex bb4 = (SQRT(a^2-2*a*d+4*b*c+d^2)+a-d)

complex bb5 = SQRT(a^2-2*a*d+4*b*c+d^2)

complex eb1 = EXP(bb1/2+a/2+d/2)

complex eb2 = EXP(bb3)

complex eb3 = EXP(a/2+d/2)

complex eb4 =EXP(bb1/2)-EXP(bb5/2)

complex dd = bb5 * (a*d - b*c )

complex dd2 = (b*c*eb1-b*c*eb2 - dd )

anew = ( d*(eb1*bb2+eb2*bb4-2*bb5)/2 + dd2 ) /dd


bnew =  ( - b*(eb1*bb2+eb2*bb4-2*bb5)/ 2 - a*b*eb3*eb4 ) /dd

cnew =  ( - c*(eb1*bb4+eb2*bb2-2*bb5)/2 - c*d*eb3*eb4 ) /dd


dnew = ( a*(eb1*bb4+eb2*bb2-2*bb5)/2 +  dd2 ) /dd


endif




a = fn1(anew)
b = fn2(bnew)
c = fn3(cnew)
d = fn4(dnew)



if @plane == 0
z  = a 
elseif @plane == 1
z  =  b
elseif @plane == 2
z  = c 
elseif @plane == 3
z  = d 
endif

  rho = cabs(a)^2+cabs(b)^2+cabs(c)^2+cabs(d)^2

bailout:
rho <= r

default:
  title = "Matrix Phi 1"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param bp1
    caption = "Init Value 2nd Entry"
    default = (0.1,0.1)
  endparam

param bp2
    caption = "Init Value 3rd Entry"
    default = (0.1,0.1)
  endparam

param bp3
    caption = "Init Value 4th Entry"
    default = (0.1,0.1)
  endparam




  param plane
   caption = "Select visible matrix entry"
   default = 0
    enum = "a" "b" "c" "d" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc

 func fn3
    caption = "Pre function3"
    default = ident()
  endfunc

 func fn4
    caption = "Pre function4"
    default = ident()
  endfunc







switch:
  type = "gwfa-c-MatrixPhi_1_Jul"
  seed = #pixel
  plane = @plane  

  bailout = @bailout
  bp1 = @bp1
  bp2 = @bp2
  bp3 = @bp3
  fn1 = fn1
  fn2 = fn2
  fn3 = fn3
  fn4 = fn4
}




gwfa-c-MatrixPhi_1_Jul {
; Iteration of phi-1 in the algebra of complex 2x2 matrices
; Julia variant
; bailout following the matrix norm Tr(aa')
init:

complex  z = pixel 

if @plane == 0
complex a = z
complex b = @bp1
complex c = @bp2
complex d = @bp3


elseif @plane == 1
complex b = z
complex c = @bp1
complex d = @bp2
complex a = @bp3


elseif @plane == 2
complex c = z
complex d = @bp1
complex a= @bp2
complex b = @bp3


elseif @plane == 3
complex d = z
complex a = @bp1
complex b = @bp2
complex c = @bp3



endif


if @disturb == 0
complex aseed = @seed
complex bseed = 0
complex cseed = 0
complex aseed = @seed

elseif @disturb == 1
complex cseed = @seed
complex aseed = 0
complex dseed = 0
complex bseed = @seed

elseif @disturb == 2
complex aseed = @seed
complex bseed = @seed
complex cseed = @seed
complex aseed = @seed

endif




float r = real(@bailout)^2

loop:


float  rho = cabs(a)^2+cabs(b)^2+cabs(c)^2+cabs(d)^2

if rho == 0

complex anew = 0
complex bnew =  0
complex cnew =  0
complex dnew = 0


else


complex bb1 = -SQRT(a^2-2*a*d+4*b*c+d^2)

complex bb2 = (SQRT(a^2-2*a*d+4*b*c+d^2)-a+d)

complex bb3 = (SQRT(a^2-2*a*d+4*b*c+d^2)/2+a/2+d/2)

complex bb4 = (SQRT(a^2-2*a*d+4*b*c+d^2)+a-d)

complex bb5 = SQRT(a^2-2*a*d+4*b*c+d^2)

complex eb1 = EXP(bb1/2+a/2+d/2)

complex eb2 = EXP(bb3)

complex eb3 = EXP(a/2+d/2)

complex eb4 =EXP(bb1/2)-EXP(bb5/2)

complex dd = bb5 * (a*d - b*c )

complex dd2 = (b*c*eb1-b*c*eb2 - dd )

anew = ( d*(eb1*bb2+eb2*bb4-2*bb5)/2 + dd2 ) /dd


bnew =  ( - b*(eb1*bb2+eb2*bb4-2*bb5)/ 2 - a*b*eb3*eb4 ) /dd

cnew =  ( - c*(eb1*bb4+eb2*bb2-2*bb5)/2 - c*d*eb3*eb4 ) /dd


dnew = ( a*(eb1*bb4+eb2*bb2-2*bb5)/2 +  dd2 ) /dd


endif


if @variant == 0
a = fn1(anew) +aseed
b = fn2(bnew) + bseed
c = fn3(cnew) + cseed
d = fn4(dnew) + dseed


elseif @variant == 1
a = fn1(anew +aseed)
b = fn2(bnew + bseed)
c = fn3(cnew + cseed)
d = fn4(dnew + dseed)

endif


if @plane == 0
z  = a 
elseif @plane == 1
z  =  b
elseif @plane == 2
z  = c 
elseif @plane == 3
z  = d 
endif

float  rho = cabs(a)^2+cabs(b)^2+cabs(c)^2+cabs(d)^2

bailout:
rho <= r

default:
  title = "Matrix Phi 1Julia"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param bp1
    caption = "Init Value 2nd Entry"
    default = (0.1,0.1)
  endparam

param bp2
    caption = "Init Value 3rd Entry"
    default = (0.1,0.1)
  endparam

param bp3
    caption = "Init Value 4th Entry"
    default = (0.1,0.1)
  endparam

param seed
    caption = "Julia Parameter"
    default = (0.1,0.1)
  endparam

  param disturb
   caption = "Julia param effects diagonal:"
   default = 0
    enum = "Main" "Minor" "All"
  endparam

  param variant
   caption = "Julia function effects :"
   default = 0
    enum = "Variant 1" "Variant 2"
  endparam





  param plane
   caption = "Select visible matrix entry"
   default = 0
    enum = "a" "b" "c" "d" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc

 func fn3
    caption = "Pre function3"
    default = ident()
  endfunc

 func fn4
    caption = "Pre function4"
    default = ident()
  endfunc







switch:
  type = "gwfa-c-MatrixPhi_1"
   plane = @plane  

  bailout = @bailout
  bp1 = @bp1
  bp2 = @bp2
  bp3 = @bp3
  fn1 = fn1
  fn2 = fn2
  fn3 = fn3
  fn4 = fn4
}





;-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-..-
; added July 4, 2000 (old formulas from GA9906.frm refreshed)

gwfa-8ap_1 {
; 
init:
complex  z = pixel 
loop:

complex a = fn1(z)

if (cabs(a) > 0)
  z  = (exp(a) - 1) / a + pixel
else
  z = 1 + pixel
endif

z  = fn2(z)

bailout:
  |z| <= @bailout
default:
  title = "Phi 8ap1"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = 100
  endparam

 func fn1
    caption = "Pre function"
    default = ident()
  endfunc  


 func fn2
    caption = "Post function"
    default = ident()
  endfunc  
    


switch:
  type = "gwfa-8ap_1_Jul"
   seed = #pixel
  fn1 = fn1
  fn2 = fn2
  bailout = @bailout

}


gwfa-8ap_1_Jul {
; 
init:
complex  z = pixel 
loop:


complex a = fn1(z * @seed ) + @seed

if (cabs(a) > 0)
  z  = (exp(a) - 1) / a + pixel
else
  z = 1 + pixel
endif

z  = fn2(z)

bailout:
  |z| <= @bailout
default:
  title = "Phi 8ap1 Julia Variant"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
 param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam

  param bailout
    caption = "Bailout Value"
    default = 100
  endparam

 func fn1
    caption = "Pre function"
    default = ident()
  endfunc  


 func fn2
    caption = "Post function"
    default = ident()
  endfunc  
    



switch:
  type = "gwfa-8ap_1"
 fn1 = fn1
  fn2 = fn2
  bailout = @bailout


}




gwfa-8apn_1 {
; 
init:
complex  z = pixel 
loop:

complex a = fn1(z + @shift1)

if (cabs(a) > 0)
  z  = (exp(a) - 1) / a + pixel
else
  z = 1 + pixel
endif

z  = fn2(z) * @blow1

bailout:
  |z| <= @bailout
default:
  title = "Phi 8apn1"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = 100
  endparam

  param shift1
    caption = "Shift Value"
    default = (0.0,0.0)
  endparam

  param blow1
    caption = "Blow Up"
    default = (1.0,0.0)
  endparam



 func fn1
    caption = "Pre function"
    default = ident()
  endfunc  


 func fn2
    caption = "Post function"
    default = ident()
  endfunc  
    


switch:
  type = "gwfa-8apn_1_Jul"
   seed = #pixel
  fn1 = fn1
  fn2 = fn2
  shift1 = @shift1
  blow1 = @blow2
  bailout = @bailout

}


gwfa-8apn_1_Jul {
; 
init:
complex  z = pixel 
loop:


complex a = fn1(z * @seed  + @shift1 ) + @seed

if (cabs(a) > 0)
  z  = (exp(a) - 1) / a + pixel
else
  z = 1 + pixel
endif

z  = fn2(z) * @blow1

bailout:
  |z| <= @bailout
default:
  title = "Phi 8apn1 Julia Variant"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
 param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam

  param bailout
    caption = "Bailout Value"
    default = 100
  endparam


  param shift1
    caption = "Shift Value"
    default = (0.0,0.0)
  endparam

  param blow1
    caption = "Blow Up"
    default = (1.0,0.0)
  endparam


 func fn1
    caption = "Pre function"
    default = ident()
  endfunc  


 func fn2
    caption = "Post function"
    default = ident()
  endfunc  
    



switch:
  type = "gwfa-8apn_1"
 fn1 = fn1
  fn2 = fn2
  shift1 = @shift1
  blow1 = @blow2
  bailout = @bailout


}




gwfa-8apn_2 {
; 
init:
complex  z = pixel 
loop:

complex a = fn1(z )
 a =  a + @shift1

if (cabs(a) > 0)
  z  = (exp(a) - 1) / a + pixel
else
  z = 1 + pixel
endif

z  = fn2(z) * @blow1

bailout:
  |z| <= @bailout
default:
  title = "Phi 8apn2"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = 100
  endparam

  param shift1
    caption = "Shift Value"
    default = (0.0,0.0)
  endparam

  param blow1
    caption = "Blow Up"
    default = (1.0,0.0)
  endparam



 func fn1
    caption = "Pre function"
    default = ident()
  endfunc  


 func fn2
    caption = "Post function"
    default = ident()
  endfunc  
    


switch:
  type = "gwfa-8apn_2_Jul"
   seed = #pixel
  fn1 = fn1
  fn2 = fn2
  shift1 = @shift1
  blow1 = @blow1
  bailout = @bailout

}


gwfa-8apn_2_Jul {
; 
init:
complex  z = pixel 
loop:


complex a = fn1(z * @seed )
 a =  a + @shift1 + @seed


if (cabs(a) > 0)
  z  = (exp(a) - 1) / a + pixel
else
  z = 1 + pixel
endif

z  = fn2(z) * @blow1

bailout:
  |z| <= @bailout
default:
  title = "Phi 8apn2 Julia Variant"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
 param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam

  param bailout
    caption = "Bailout Value"
    default = 100
  endparam


  param shift1
    caption = "Shift Value"
    default = (0.0,0.0)
  endparam

  param blow1
    caption = "Blow Up"
    default = (1.0,0.0)
  endparam


 func fn1
    caption = "Pre function"
    default = ident()
  endfunc  


 func fn2
    caption = "Post function"
    default = ident()
  endfunc  
    



switch:
  type = "gwfa-8apn_2"
 fn1 = fn1
  fn2 = fn2
  shift1 = @shift1
  blow1 = @blow1
  bailout = @bailout


}





gwfa-4recon7_1 {
; 
init:
complex  z = pixel 
loop:

complex a = fn1(z ) + @shift1 

  z  =  a / (exp(a + Pixel ) - 1)  + pixel

z  = fn2(z) * @blow1

bailout:
  |z| <= @bailout
default:
  title = "Phi 4recon71"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = 100
  endparam

  param shift1
    caption = "Shift Value"
    default = (0.0,0.0)
  endparam

  param blow1
    caption = "Blow Up"
    default = (1.0,0.0)
  endparam



 func fn1
    caption = "Pre function"
    default = ident()
  endfunc  


 func fn2
    caption = "Post function"
    default = ident()
  endfunc  
    


switch:
  type = "gwfa-4recon7_1_Jul"
   seed = #pixel
  fn1 = fn1
  fn2 = fn2
  shift1 = @shift1
  blow1 = @blow2
  bailout = @bailout

}


gwfa-4recon7_1_Jul {
; 
init:
complex  z = pixel 
loop:




complex a = fn1(z  * @seed ) + @shift1 + @seed

  z  =  a / (exp(a + Pixel ) - 1)  + pixel

z  = fn2(z) * @blow1



bailout:
  |z| <= @bailout
default:
  title = "Phi 4recon71 Julia Variant"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
 param seed
    caption = "fix perturbation"
    default = (1,0)
  endparam

  param bailout
    caption = "Bailout Value"
    default = 100
  endparam


  param shift1
    caption = "Shift Value"
    default = (0.0,0.0)
  endparam

  param blow1
    caption = "Blow Up"
    default = (1.0,0.0)
  endparam


 func fn1
    caption = "Pre function"
    default = ident()
  endfunc  


 func fn2
    caption = "Post function"
    default = ident()
  endfunc  
    



switch:
  type = "gwfa-4recon7_1"
 fn1 = fn1
  fn2 = fn2
  shift1 = @shift1
  blow1 = @blow2
  bailout = @bailout


}

;-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-..-












;-------------------------------------------------------
;WwWwWwWwWwWwWwWwWwWwWwWwWwWwWwWwWwWwWwWwWwW
;-------------------------------------------------------
;The formulas below are only for compatibility purposes and should not be used any more
; In public they will be needed only for the UF parameter collection of June 1999

Phi_4_trial3 {; GA
init:
  z = Pixel , d = p3 * Pixel
loop:

  z = z + p1 
  z = fn1(z)/(exp(fn2(z))-1) + d 
bailout:
    |z| <= p2

default:
  title = "z-look at the parameter file"
  maxiter = 100
  periodicity = 2
  center = (0,0)
  magn = 1.0
  
  }


Phi_4_recon7 {; GA
init:
  z = Pixel  
loop:
  u = fn1(z)
  u = u + p1
  z = u / (exp(u + Pixel ) - 1 ) + Pixel
bailout:
    |z| <= p2
default:
  title = "z-look at the parameter file"
  maxiter = 100
  periodicity = 2
  center = (0,0)
  magn = 1.0
  
  }



Phi_4_recon14 {; GA
init:
  z = Pixel  
loop:

  z = z + p1
  u = fn1(z)
  z =  z/(exp(u) - 1  + Pixel)  
bailout:
    |z| <= p2
default:
  title = "z-look at the parameter file"
  maxiter = 100
  periodicity = 2
  center = (0,0)
  magn = 1.0
  
  }

Phi_8ap_n {; GA
init:
  z = Pixel  
loop:
  a = fn1(z)
  z = z + p1 
IF ( |a| > 0 )
  z = (exp(a)-1)/a + Pixel
ELSE
  z = 1 + Pixel 
ENDIF
  z = fn2(z)* p2
bailout:
    |z| <= p3
default:
  title = "z-look at the parameter file"
  maxiter = 100
  periodicity = 2
  center = (0,0)
  magn = 1.0
  
  }



Phi_4_recon051301 {; GA
init:
  z = Pixel  

loop:
;   z = z + p1
   z =  (exp(fn1(z ))  - 1  )/(z+p1)   + Pixel
bailout:
    |z| <= p2
default:
  title = "z-look at the parameter file"
  maxiter = 100
  periodicity = 2
  center = (0,0)
  magn = 1.0
  
  }


; following are earlier versions of the quaternion formulas which did not handle the plane choice appropriately but gave nonetheless
; nice images

gwfa-Quat_1 {
; Quaternionic version of Phi1 (with partial sum for n=3)
init:
complex  z = pixel 
float a = real(z)
float b = imag(z)
float c = real(fn1(z))
float d = imag(fn1(z))
loop:

float t1 = (3*a^2+8*a-b^2-c^2-d^2 +12)/24


float anew = (a^3+4*a^2-3*a*(b^2+c^2+d^2-4)-4*(b^2+c^2+d^2-6))/24 
float bnew = b*t1
float cnew = c*t1
float dnew = d*t1

a = anew
b = bnew
c = cnew
d = dnew



if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  =  b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

bailout:
  |z| <= @bailout
default:
  title = "z-Quaternionic Phi  1"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-Quat_1_Jul"
  seed = #pixel
  plane = @plane  
  bailout = @bailout
  fn1 = fn1
}

gwfa-Quat_1_Jul {
; Quaternionic version of Phi1 (with partial sum for n=3) - Julia variant
init:
complex  z = pixel 
float a = real(z)
float b = imag(z)
float c = real(fn1(z))
float d = imag(fn1(z))

float u = real(@seed)
float v = imag(@seed)

loop:


float t1 = (3*a^2+8*a-b^2-c^2-d^2 +12)/24

float anew = (a^3+4*a^2-3*a*(b^2+c^2+d^2-4)-4*(b^2+c^2+d^2-6))/24 + u
float bnew = b*t1 + v
float cnew = c*t1 + u 
float dnew = d*t1 + v

a = anew
b = bnew
c = cnew
d = dnew

if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  = b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif


bailout:

|z| <= @bailout

default:
  title = "z-Quaternionic Phi  1 Julia "
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 

  endparam
 



  param seed
    caption = "fixed perturbation"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-Quat_1"
    plane = @plane
  bailout = @bailout
  fn1 = fn1

}

gwfa-Quat_1a {
; Quaternionic version of Phi1 (with partial sum for n=3)
; bailout following the quaternion norm
init:
complex  z = pixel 
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))
float r = real(@bailout^2)
loop:

float t1 = (3*a^2+8*a-b^2-c^2-d^2 +12)/24


float anew = (a^3+4*a^2-3*a*(b^2+c^2+d^2-4)-4*(b^2+c^2+d^2-6))/24 
float bnew = b*t1
float cnew = c*t1
float dnew = d*t1

a = anew
b = bnew
c = cnew
d = dnew



if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  =  b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float  rho = a^2 +b^2+c^2+d^2

bailout:
rho <= r

default:
  title = "z-Quaternionic Phi  1a"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-Quat_1a_Jul"
  seed = #pixel
  plane = @plane  
  bailout = @bailout
  dim34 = @dim34
  fn1 = fn1
}

gwfa-Quat_1a_Jul {
; Quaternionic version of Phi1 (with partial sum for n=3) - Julia variant
; bailout following the quaternion norm
init:
complex  z = pixel 
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))

float r = real(@bailout^2)

float u = real(@seed)
float v = imag(@seed)

loop:


float t1 = (3*a^2+8*a-b^2-c^2-d^2 +12)/24

float anew = (a^3+4*a^2-3*a*(b^2+c^2+d^2-4)-4*(b^2+c^2+d^2-6))/24 + u
float bnew = b*t1 + v
float cnew = c*t1 + u 
float dnew = d*t1 + v

a = anew
b = bnew
c = cnew
d = dnew

if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  = b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float rho = a^2 + b^2 + c^2 + d^2

bailout:
rho <= r
default:
  title = "z-Quaternionic Phi  1a Julia "
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 

  endparam
 



  param seed
    caption = "fixed perturbation"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-Quat_1a"
    plane = @plane
  bailout = @bailout
  dim34 = @dim34
  fn1 = fn1

}



gwfa-Quat_2 {
; Quaternionic version of Phi1 (with partial sum for n=7)
init:
complex  z = pixel 
float a = real(z)
float b = imag(z)
float c = real(fn1(z))
float d = imag(fn1(z))
loop:


float t1 = (3*b^2+3*c^2+3*d^2-8)
float t2 = (5*b^2+5*c^2+5*d^2-14)
float t3 = (b^4+2*b^2*(c^2+d^2-8)+c^4+2*c^2*(d^2-8)+d^4-16*d^2+48)
float t4 = (5*b^4+2*b^2*(5*c^2+5*d^2-42)+5*c^4+2*c^2*(5*d^2-42)+5*d^4-84*d^2+280)
float t5 = (c^4+2*c^2*(d^2-14)+d^4-28*d^2+280)

float t6 =(b^2+c^2+d^2-8)
float t7 = (5*b^2+5*c^2+5*d^2-42)
float t8 = (3*b^4+2*b^2*(3*c^2+3*d^2-40)+3*c^4+2*c^2*(3*d^2-40)+3*d^4-80*d^2+720)
float t9 = (b^4+2*b^2*(c^2+d^2-14)+c^4+2*c^2*(d^2-14)+d^4-28*d^2+280)
float t10 = (3*c^4+2*c^2*(3*d^2-56)+3*d^4-112*d^2+1680)
float t11 = (3*d^4-112*d^2+1680)
float t13 = (3*c^2+3*d^2-40)
float t14 = (3*c^4+2*c^2*(3*d^2-40)+3*d^4-80*d^2+720)
float t15 = (3*d^4-80*d^2+720)
float t16 = (b^6+3*b^4*(c^2+d^2-14)+3*b^2*t5 +c^6+3*c^4*(d^2-14)+3*c^2*(d^4-28*d^2+280)+d^6-42*d^4+840*d^2-5040)

float t17=(7*a^6+48*a^5-35*a^4*t6-32*a^3*t7+7*a^2*t8+48*a*t9-b^6-b^4*(3*c^2+3*d^2-56)-b^2*t10-c^6+c^4*(56-3*d^2)-c^2*t11-d^6+56*d^4-1680*d^2+20160)/40320



float anew = (a^7+8*a^6-7*a^5*t1-24*a^4*t2+35*a^3*t3+24*a^2*t4-7*a*(b^6+b^4*t13+b^2*t14+c^6+c^4*(3*d^2-40)+c^2*t15+d^6-40*d^4+720*d^2-2880)-8*t16)/40320


float bnew  = b*t17

float cnew  = c*t17

float dnew  = d*t17

a = anew
b = bnew
c = cnew
d = dnew


if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  =  b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

bailout:
  |z| <= @bailout
default:
  title = "z-Quaternionic Phi  2"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 1
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc



switch:
  type = "gwfa-Quat_2_Jul"
  seed = #pixel
  plane = @plane  
  bailout = @bailout
  fn1 = fn1

}

gwfa-Quat_2_Jul {
; Quaternionic version of Phi1 (with partial sum for n=7) - Julia variant
init:
complex  z = pixel 
float a = real(z)
float b = imag(z)
float c = real(fn1(z))
float d = imag(fn1(z))

float u = real(@seed)
float v = imag(@seed)
float  w = real(@seed)
float s = imag(@seed)


loop:



float t1 = (3*b^2+3*c^2+3*d^2-8)
float t2 = (5*b^2+5*c^2+5*d^2-14)
float t3 = (b^4+2*b^2*(c^2+d^2-8)+c^4+2*c^2*(d^2-8)+d^4-16*d^2+48)
float t4 = (5*b^4+2*b^2*(5*c^2+5*d^2-42)+5*c^4+2*c^2*(5*d^2-42)+5*d^4-84*d^2+280)
float t5 = (c^4+2*c^2*(d^2-14)+d^4-28*d^2+280)

float t6 =(b^2+c^2+d^2-8)
float t7 = (5*b^2+5*c^2+5*d^2-42)
float t8 = (3*b^4+2*b^2*(3*c^2+3*d^2-40)+3*c^4+2*c^2*(3*d^2-40)+3*d^4-80*d^2+720)
float t9 = (b^4+2*b^2*(c^2+d^2-14)+c^4+2*c^2*(d^2-14)+d^4-28*d^2+280)
float t10 = (3*c^4+2*c^2*(3*d^2-56)+3*d^4-112*d^2+1680)
float t11 = (3*d^4-112*d^2+1680)
float t13 = (3*c^2+3*d^2-40)
float t14 = (3*c^4+2*c^2*(3*d^2-40)+3*d^4-80*d^2+720)
float t15 = (3*d^4-80*d^2+720)
float t16 = (b^6+3*b^4*(c^2+d^2-14)+3*b^2*t5 +c^6+3*c^4*(d^2-14)+3*c^2*(d^4-28*d^2+280)+d^6-42*d^4+840*d^2-5040)

float t17=(7*a^6+48*a^5-35*a^4*t6-32*a^3*t7+7*a^2*t8+48*a*t9-b^6-b^4*(3*c^2+3*d^2-56)-b^2*t10-c^6+c^4*(56-3*d^2)-c^2*t11-d^6+56*d^4-1680*d^2+20160)/40320


float anew = (a^7+8*a^6-7*a^5*t1-24*a^4*t2+35*a^3*t3+24*a^2*t4-7*a*(b^6+b^4*t13+b^2*t14+c^6+c^4*(3*d^2-40)+c^2*t15+d^6-40*d^4+720*d^2-2880)-8*t16)/40320

float bnew  = b*t17
float cnew = c*t17
float dnew  = d*t17



a = anew + u
b = bnew + v
c = cnew + w
d = dnew + s


if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  = b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif



bailout:
  |z| <= @bailout
default:
  title = "z-Quaternionic Phi  2 Julia "
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 1
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 

  endparam
 



  param seed
    caption = "fixed perturbation"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-Quat_2"
    plane = @plane
  bailout = @bailout
  fn1 = fn1

}





gwfa-QuatPhi1 {
; full Quaternionic version of Phi1 
; bailout following the quaternion norm
init:
complex  z = pixel 
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))
float r = real(@bailout)^2
loop:



float norm1 = a^2+b^2+c^2+d^2
float norm2 = SQRT(b^2 + c^2 + d^2)
float s1 = SIN(norm2)
float c1 = COS(norm2)
float e1 = exp(a)
float f1 = (e1*(a*s1/(norm1*norm2)-c1/norm1)+ 1/norm1 )



float anew = (e1*(a*c1+norm2*s1)-a)/norm1
float bnew = b* f1
float cnew = c* f1
float dnew = d * f1




a = anew
b = bnew
c = cnew
d = dnew



if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  =  b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float  rho = a^2 +b^2+c^2+d^2

bailout:
rho <= r

default:
  title = "z-Quaternionic Phi 1Full"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-QuatPhi1_Jul"
  seed = #pixel
  plane = @plane  
  bailout = @bailout
  dim34 = @dim34
  fn1 = fn1
}

gwfa-QuatPhi1_Jul {
; full Quaternionic version of Phi1 - Julia variant
; bailout following the quaternion norm
init:
complex  z = pixel 
float a = real(z)
float b = imag(z)
float c = real(fn1(@dim34))
float d = imag(fn1(@dim34))

float r = real(@bailout)^2

float u = real(@seed)
float v = imag(@seed)

loop:


float norm1 = a^2+b^2+c^2+d^2
float norm2 = SQRT(b^2 + c^2 + d^2)
float s1 = SIN(norm2)
float c1 = COS(norm2)
float e1 = exp(a)
float f1 = (e1*(a*s1/(norm1*norm2)-c1/norm1)+ 1/norm1 )



float anew = (e1*(a*c1+norm2*s1)-a)/norm1 
float bnew = b* f1 
float cnew = c* f1 
float dnew = d * f1 




a = anew + u 
b = bnew + v 
c = cnew + u 
d = dnew + v

if @plane == 0
z  = a + flip(b)
elseif @plane == 1
z  = b + flip(c)
elseif @plane == 2
z  = c + flip(d)
elseif @plane == 3
z  = d + flip(a)
elseif @plane == 4
z  = a + flip(c)
elseif @plane == 5
z  = b + flip(d)
endif

float rho = a^2 + b^2 + c^2 + d^2

bailout:
rho <= r
default:
  title = "z-Quaternionic Phi 1Full Jul"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param dim34
    caption = "Init Value c, d"
    default = (1,1)
  endparam


  param plane
   caption = "Select 2dim plane"
   default = 0
    enum = "ab" "bc" "cd" "da" "ac" "bd" 

  endparam
 



  param seed
    caption = "fixed perturbation"
    default = (1,0)
  endparam

 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc


switch:
  type = "gwfa-QuatPhi1"
    plane = @plane
  bailout = @bailout
  dim34 = @dim34
  fn1 = fn1

}












;-------------------------------------------------------------------------------------------------
comment {
  gwfa-attractor.ufm  Vers. 1.01.02a
  Fractal types for Ultra Fractal 2
  by G.W.F. Albrecht
  February 3, 2001

  For more information about this formula collection,
  please visit its home page:
  
      http://home.wtal.de/spiriteye/fractal/index.html

or

     http://www.gwfa.de

}


Cockatoo {; G.W.F. Albrecht, February 2001
; partly inspired from Ron Barnetts Strange Attractor Formulas
; 
init:
  float xx = 0
  float x = real(#pixel)
  float y = imag(#pixel)
  
  float xabs = cabs(x)
  bool done = false
  int touche_it = 0
loop:
  
xx   = x
x = y*(@a+SIN(@b*xx))-@c*SQRT(xabs)
y = @d-xx

#z = x + flip(y)
complex u = #z
IF cabs(#z-#pixel) <= @epsilon
   touche_it = touche_it + 1
   IF (touche_it >= @max_touch)
      done = true
      IF @variant == 0
       #z = cabs(#z-#pixel) + flip(cabs(#z))
      ELSEIF @variant == 1
       #z = #z
      ELSEIF @variant == 2
       IF cabs(#z) > 0
        #z = @fn1((EXP(u) - 1) / u)
       ELSE 
        #z = @fn1(1)
       ENDIF
      ELSEIF @variant == 3
       IF cabs(#z) > 0
        #z = @fn1(u /(EXP(u) - 1) )
       ELSE 
        #z = @fn1(1)
       ENDIF

      ENDIF
    ENDIF
ENDIF


bailout:
  done == false
default:
  title = "Cockatoo Attractor"
  center = (0.0, 0.0)
  maxiter = 1000
  periodicity = 0
  method = onepass
param a
  caption = "Parameter a"
  default = 1.0
endparam
param b
  caption = "Parameter b"
  default = 0.7
endparam
param c
  caption = "Parameter c"
  default = 1.2
endparam
param d
  caption = "Parameter d"
  default = 0.21
endparam
param epsilon
  caption = "epsilon"
  default = 0.5
endparam
param max_touch
  caption = "find pixel attempts"
  default = 2
endparam
param variant
  caption = "Variant"
  enum = "A" "B" "C" "D"
  default = 0
endparam
 func fn1
    caption = "Pixel function"
    default = ident()
  endfunc  



switch:
  type = "Cockatoo2"
  seed = #pixel
  a = @a
  b = @b
  c = @c
  d = @d
  epsilon = @epsilon
  max_touch = @max_touch
  variant = @variant
  fn1 = fn1

}




Cockatoo2 {; G.W.F. Albrecht, February 2001
; partly inspired from Ron Barnetts Strange Attractor Formulas
; 
init:
  float xx = 0
  float x = real(#pixel)
  float y = imag(#pixel)
  
  float xabs = cabs(x)
  bool done = false
  int touche_it = 0
loop:
  
xx   = x
x = y*(@a+SIN(@b*xx))-@c*SQRT(xabs)
y = @d-xx

#z = x + flip(y) + @fn2(@seed)
complex u = #z
IF cabs(#z-#pixel) <= @epsilon
   touche_it = touche_it + 1
   IF (touche_it >= @max_touch)
      done = true
      IF @variant == 0
       #z = cabs(#z-#pixel) + flip(cabs(#z))
      ELSEIF @variant == 1
       #z = #z
      ELSEIF @variant == 2
       IF cabs(#z) > 0
        #z = @fn1((EXP(u) - 1) / u)
       ELSE 
        #z = @fn1(1)
       ENDIF
      ELSEIF @variant == 3
       IF cabs(#z) > 0
        #z = @fn1(u /(EXP(u) - 1) )
       ELSE 
        #z = @fn1(1)
       ENDIF

      ENDIF
    ENDIF
ENDIF


bailout:
  done == false
default:
  title = "Cockatoo Attractor Julia Variant"
  center = (0.0, 0.0)
  maxiter = 1000
  periodicity = 0
  method = onepass
param a
  caption = "Parameter a"
  default = 1.0
endparam
param b
  caption = "Parameter b"
  default = 0.7
endparam
param c
  caption = "Parameter c"
  default = 1.2
endparam
param d
  caption = "Parameter d"
  default = 0.21
endparam
param epsilon
  caption = "epsilon"
  default = 0.5
endparam
param max_touch
  caption = "find pixel attempts"
  default = 2
endparam
param variant
  caption = "Variant"
  enum = "A" "B" "C" "D"
  default = 0
endparam

param seed
  caption = " Julia Perturbation"
  default =( 0,0)
endparam


 func fn1
    caption = "Pixel function"
    default = ident()
  endfunc  


 func fn2
    caption = "Seed function"
    default = ident()
  endfunc  




switch:
  type = "Cockatoo"

  a = @a
  b = @b
  c = @c
  d = @d
  epsilon = @epsilon
  max_touch = @max_touch
  variant = @variant
  fn1 = fn1

}










;-------------------------------------------------------------------------------------------------------------

comment {
  gwfa-matrix.ufm  Vers. 1.01.06a
  Fractal types for Ultra Fractal 2
  by G.W.F. Albrecht
  June 1, 2001

  For more information about this formula collection,
  please visit its home page:
  
      http://home.wtal.de/spiriteye/fractal/index.html

or

     http://www.gwfa.de

 This collection uses formulas derived from the transformation z -> z^n
 in the algebra of complex 2x2 or 3x3 matrices.



}




gwfa-c-MatrixPlain_1 {
; Iteration of z^2 in the algebra of complex 2x2 matrices
; bailout following the matrix norm Tr(aa')
; matrix notation like [[a,b],[c,d]]
init:

complex  z = pixel 

if @plane == 0
complex a = z
complex b = @bp1
complex c = @bp2
complex d = @bp3


elseif @plane == 1
complex b = z
complex c = @bp1
complex d = @bp2
complex a = @bp3


elseif @plane == 2
complex c = z
complex d = @bp1
complex a= @bp2
complex b = @bp3


elseif @plane == 3
complex d = z
complex a = @bp1
complex b = @bp2
complex c = @bp3



endif



float r = real(@bailout)^2

loop:


complex anew = a*a  + b *c
complex bnew = b*(a +  d )
complex cnew = c* (a + d) 
complex dnew = c*b + d*d







a = fn1(anew)
b = fn2(bnew)
c = fn3(cnew)
d = fn4(dnew)



if @plane == 0
z  = a 
elseif @plane == 1
z  =  b
elseif @plane == 2
z  = c 
elseif @plane == 3
z  = d 
endif

float   rho = cabs(a)^2+cabs(b)^2+cabs(c)^2+cabs(d)^2

bailout:
rho <= r

default:
  title = "Matrix Plain 1"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param bp1
    caption = "Init Value 2nd Entry"
    default = (0.1,0.1)
  endparam

param bp2
    caption = "Init Value 3rd Entry"
    default = (0.1,0.1)
  endparam

param bp3
    caption = "Init Value 4th Entry"
    default = (0.1,0.1)
  endparam




  param plane
   caption = "Select visible matrix entry"
   default = 0
    enum = "a" "b" "c" "d" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc

 func fn3
    caption = "Pre function3"
    default = ident()
  endfunc

 func fn4
    caption = "Pre function4"
    default = ident()
  endfunc







switch:
  type = "gwfa-c-MatrixPlain_1_Jul"
  seed = #pixel
  plane = @plane  

  bailout = @bailout
  bp1 = @bp1
  bp2 = @bp2
  bp3 = @bp3
  fn1 = fn1
  fn2 = fn2
  fn3 = fn3
  fn4 = fn4
}




gwfa-c-MatrixPlain_1_Jul {
; Iteration of z^2 in the algebra of complex 2x2 matrices
; Julia variant
; bailout following the matrix norm Tr(aa')
init:

complex  z = pixel 

if @plane == 0
complex a = z
complex b = @bp1
complex c = @bp2
complex d = @bp3


elseif @plane == 1
complex b = z
complex c = @bp1
complex d = @bp2
complex a = @bp3


elseif @plane == 2
complex c = z
complex d = @bp1
complex a= @bp2
complex b = @bp3


elseif @plane == 3
complex d = z
complex a = @bp1
complex b = @bp2
complex c = @bp3



endif


if @disturb == 0
complex aseed = @seed
complex bseed = 0
complex cseed = 0
complex aseed = @seed

elseif @disturb == 1
complex cseed = @seed
complex aseed = 0
complex dseed = 0
complex bseed = @seed

elseif @disturb == 2
complex aseed = @seed
complex bseed = @seed
complex cseed = @seed
complex aseed = @seed

endif




float r = real(@bailout)^2

loop:

complex anew = a*a  + b *c
complex bnew = b*(a +  d )
complex cnew = c* (a + d) 
complex dnew = c*b + d*d



if @variant == 0
a = fn1(anew) +aseed
b = fn2(bnew) + bseed
c = fn3(cnew) + cseed
d = fn4(dnew) + dseed


elseif @variant == 1
a = fn1(anew +aseed)
b = fn2(bnew + bseed)
c = fn3(cnew + cseed)
d = fn4(dnew + dseed)

endif


if @plane == 0
z  = a 
elseif @plane == 1
z  =  b
elseif @plane == 2
z  = c 
elseif @plane == 3
z  = d 
endif

float  rho = cabs(a)^2+cabs(b)^2+cabs(c)^2+cabs(d)^2

bailout:
rho <= r

default:
  title = "Matrix Plain 1Julia"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param bp1
    caption = "Init Value 2nd Entry"
    default = (0.1,0.1)
  endparam

param bp2
    caption = "Init Value 3rd Entry"
    default = (0.1,0.1)
  endparam

param bp3
    caption = "Init Value 4th Entry"
    default = (0.1,0.1)
  endparam

param seed
    caption = "Julia Parameter"
    default = (0.1,0.1)
  endparam

  param disturb
   caption = "Julia param effects diagonal:"
   default = 0
    enum = "Main" "Minor" "All"
  endparam

  param variant
   caption = "Julia function effects :"
   default = 0
    enum = "Variant 1" "Variant 2"
  endparam





  param plane
   caption = "Select visible matrix entry"
   default = 0
    enum = "a" "b" "c" "d" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc

 func fn3
    caption = "Pre function3"
    default = ident()
  endfunc

 func fn4
    caption = "Pre function4"
    default = ident()
  endfunc







switch:
  type = "gwfa-c-MatrixPlain_1"
   plane = @plane  

  bailout = @bailout
  bp1 = @bp1
  bp2 = @bp2
  bp3 = @bp3
  fn1 = fn1
  fn2 = fn2
  fn3 = fn3
  fn4 = fn4
}


;------------------------------------

gwfa-c-MatrixPower_1 {
; Iteration of z^n in the algebra of complex 2x2 matrices
; bailout following the matrix norm Tr(aa')
; matrix notation like [[a,b],[c,d]]
init:

complex  z = pixel 

if @plane == 0
complex a = z
complex b = @bp1
complex c = @bp2
complex d = @bp3


elseif @plane == 1
complex b = z
complex c = @bp1
complex d = @bp2
complex a = @bp3


elseif @plane == 2
complex c = z
complex d = @bp1
complex a= @bp2
complex b = @bp3


elseif @plane == 3
complex d = z
complex a = @bp1
complex b = @bp2
complex c = @bp3



endif



float r = real(@bailout)^2

loop:



complex expr1 = SQRT(a^2-2*a*d+4*b*c+d^2)
complex expr2 = 1 / ( 2^@mpower)
complex expr3 = expr2 / 2


complex anew = (expr3*(expr1-a+d)*(-expr1+a+d)^@mpower+expr3*(expr1+a-d)*(expr1+a+d)^@mpower) /expr1
complex bnew = expr2*b*((expr1+a+d)^@mpower-(-expr1+a+d)^@mpower)/expr1
complex cnew = expr2*c*((expr1+a+d)^@mpower-(-expr1+a+d)^@mpower)/expr1
complex dnew = (expr3*(expr1+a-d)*(-expr1+a+d)^@mpower + expr3*(expr1-a+d)*(expr1+a+d)^@mpower)/expr1








a = fn1(anew)
b = fn2(bnew)
c = fn3(cnew)
d = fn4(dnew)



if @plane == 0
z  = a 
elseif @plane == 1
z  =  b
elseif @plane == 2
z  = c 
elseif @plane == 3
z  = d 
endif

float   rho = cabs(a)^2+cabs(b)^2+cabs(c)^2+cabs(d)^2

bailout:
rho <= r

default:
  title = "Matrix Power 1"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param bp1
    caption = "Init Value 2nd Entry"
    default = (0.1,0.1)
  endparam

param bp2
    caption = "Init Value 3rd Entry"
    default = (0.1,0.1)
  endparam

param bp3
    caption = "Init Value 4th Entry"
    default = (0.1,0.1)
  endparam

param mpower
    caption = "Complex Power"
    default = (3,0)
  endparam




  param plane
   caption = "Select visible matrix entry"
   default = 0
    enum = "a" "b" "c" "d" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc

 func fn3
    caption = "Pre function3"
    default = ident()
  endfunc

 func fn4
    caption = "Pre function4"
    default = ident()
  endfunc







switch:
  type = "gwfa-c-MatrixPower_1_Jul"
  seed = #pixel
  plane = @plane  

  bailout = @bailout
  bp1 = @bp1
  bp2 = @bp2
  bp3 = @bp3
  mpower = @mpower
  fn1 = fn1
  fn2 = fn2
  fn3 = fn3
  fn4 = fn4
}




gwfa-c-MatrixPower_1_Jul {
; Iteration of z^n in the algebra of complex 2x2 matrices
; Julia variant
; bailout following the matrix norm Tr(aa')
init:

complex  z = pixel 

if @plane == 0
complex a = z
complex b = @bp1
complex c = @bp2
complex d = @bp3


elseif @plane == 1
complex b = z
complex c = @bp1
complex d = @bp2
complex a = @bp3


elseif @plane == 2
complex c = z
complex d = @bp1
complex a= @bp2
complex b = @bp3


elseif @plane == 3
complex d = z
complex a = @bp1
complex b = @bp2
complex c = @bp3



endif


if @disturb == 0
complex aseed = @seed
complex bseed = 0
complex cseed = 0
complex aseed = @seed

elseif @disturb == 1
complex cseed = @seed
complex aseed = 0
complex dseed = 0
complex bseed = @seed

elseif @disturb == 2
complex aseed = @seed
complex bseed = @seed
complex cseed = @seed
complex aseed = @seed

endif




float r = real(@bailout)^2

loop:


complex expr1 = SQRT(a^2-2*a*d+4*b*c+d^2)
complex expr2 = 1 / ( 2^@mpower)
complex expr3 = expr2 / 2


complex anew = (expr3*(expr1-a+d)*(-expr1+a+d)^@mpower+expr3*(expr1+a-d)*(expr1+a+d)^@mpower) /expr1
complex bnew = expr2*b*((expr1+a+d)^@mpower-(-expr1+a+d)^@mpower)/expr1
complex cnew = expr2*c*((expr1+a+d)^@mpower-(-expr1+a+d)^@mpower)/expr1
complex dnew = (expr3*(expr1+a-d)*(-expr1+a+d)^@mpower + expr3*(expr1-a+d)*(expr1+a+d)^@mpower)/expr1





if @variant == 0
a = fn1(anew) +aseed
b = fn2(bnew) + bseed
c = fn3(cnew) + cseed
d = fn4(dnew) + dseed


elseif @variant == 1
a = fn1(anew +aseed)
b = fn2(bnew + bseed)
c = fn3(cnew + cseed)
d = fn4(dnew + dseed)

endif


if @plane == 0
z  = a 
elseif @plane == 1
z  =  b
elseif @plane == 2
z  = c 
elseif @plane == 3
z  = d 
endif

float  rho = cabs(a)^2+cabs(b)^2+cabs(c)^2+cabs(d)^2

bailout:
rho <= r

default:
  title = "Matrix Power 1Julia"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = (100,0)
  endparam

  param bp1
    caption = "Init Value 2nd Entry"
    default = (0.1,0.1)
  endparam

param bp2
    caption = "Init Value 3rd Entry"
    default = (0.1,0.1)
  endparam

param bp3
    caption = "Init Value 4th Entry"
    default = (0.1,0.1)
  endparam

param mpower
    caption = "Complex Power"
    default = (3,0)
  endparam

param seed
    caption = "Julia Parameter"
    default = (0.1,0.1)
  endparam

  param disturb
   caption = "Julia param effects diagonal:"
   default = 0
    enum = "Main" "Minor" "All"
  endparam

  param variant
   caption = "Julia function effects :"
   default = 0
    enum = "Variant 1" "Variant 2"
  endparam





  param plane
   caption = "Select visible matrix entry"
   default = 0
    enum = "a" "b" "c" "d" 
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc

 func fn3
    caption = "Pre function3"
    default = ident()
  endfunc

 func fn4
    caption = "Pre function4"
    default = ident()
  endfunc







switch:
  type = "gwfa-c-MatrixPower_1"
   plane = @plane  

  bailout = @bailout
  bp1 = @bp1
  bp2 = @bp2
  bp3 = @bp3
  mpower = @mpower
  fn1 = fn1
  fn2 = fn2
  fn3 = fn3
  fn4 = fn4
}





gwfa-c-MatrixPlain_3 {
; Iteration of z^2 in the algebra of complex 3x3 matrices
; bailout following the matrix norm Tr(aa')
; matrix notation like [[a,b,c],[d,eh,f], [g,h,i]]
init:

complex  z = pixel 

if @plane == 0
complex a = z
complex b = @bp1
complex c = @bp2
complex d = @bp3
complex eh = @bp4
complex f = @bp5
complex g = @bp6
complex h = @bp7
complex i = @bp8




elseif @plane == 1
complex b = z
complex c = @bp1
complex d = @bp2
complex eh = @bp3
complex f = @bp4
complex g = @bp5
complex h = @bp6
complex i = @bp7
complex a = @bp8


elseif @plane == 2
complex c = z
complex d = @bp1
complex eh = @bp2
complex f = @bp3
complex g = @bp4
complex h = @bp5
complex i = @bp6
complex a = @bp7
complex b = @bp8



elseif @plane == 3
complex d = z
complex eh = @bp1
complex f = @bp2
complex g = @bp3
complex h = @bp4
complex i = @bp5
complex a = @bp6
complex b = @bp7
complex c = @bp8



elseif @plane == 4
complex eh = z
complex f = @bp1
complex g = @bp2
complex h = @bp3
complex i = @bp4
complex a = @bp5
complex b = @bp6
complex c = @bp7
complex d = @bp8



elseif @plane == 5
complex f = z
complex g = @bp1
complex h = @bp2
complex i = @bp3
complex a = @bp4
complex b = @bp5
complex c = @bp6
complex d = @bp7
complex eh = @bp8



elseif @plane == 6
complex g = z
complex h = @bp1
complex i = @bp2
complex a = @bp3
complex b = @bp4
complex c = @bp5
complex d = @bp6
complex eh = @bp7
complex f = @bp8



elseif @plane == 7
complex h = z
complex i = @bp1
complex a = @bp2
complex b = @bp3
complex c = @bp4
complex d = @bp5
complex eh = @bp6
complex f = @bp7
complex g = @bp8


elseif @plane == 8
complex i = z
complex a = @bp1
complex b = @bp2
complex c = @bp3
complex d = @bp4
complex eh = @bp5
complex f = @bp6
complex g = @bp7
complex h = @bp8



endif



float r = real(@bailout)^2

loop:




complex anew = a^2+b*d+c*g

complex bnew = a*b+b*eh+c*h

complex cnew = a*c+b*f+c*i

complex dnew = a*d+d*eh+f*g

complex enew = b*d+eh^2+f*h

complex fnew = c*d+eh*f+f*i

complex gnew = a*g+d*h+g*i

complex hnew = b*g+eh*h+h*i

complex inew = c*g+f*h+i^2









a = fn1(anew)
b = fn2(bnew)
c = fn3(cnew)
d = fn1(dnew)

eh = fn2(enew)
f = fn3(fnew)
g = fn1(gnew)
h = fn2(hnew)

i = fn3(inew)


if @plane == 0
z  = a 
elseif @plane == 1
z  =  b
elseif @plane == 2
z  = c 
elseif @plane == 3
z  = d 
elseif @plane == 4
z  =  eh
elseif @plane == 5
z  = f 
elseif @plane == 6
z  = g 
elseif @plane == 7
z  =  h
elseif @plane == 8
z  = i
endif

float   rho = cabs(a)^2+cabs(b)^2+cabs(c)^2+cabs(d)^2
            + cabs(eh)^2+cabs(f)^2+cabs(g)^2+cabs(h)^2 + cabs(i)^2

bailout:
rho <= r

default:
  title = "Matrix Plain 3x3"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = 100
  endparam

  param bp1
    caption = "Init Value 2nd Entry"
    default = (0,0)
  endparam

param bp2
    caption = "Init Value 3rd Entry"
    default = (0,0)
  endparam

param bp3
    caption = "Init Value 4th Entry"
    default = (0,0)
  endparam


param bp4
    caption = "Init Value 5th Entry"
    default = (0.1,0.1)
  endparam


param bp5
    caption = "Init Value 6th Entry"
    default = (0,0)
  endparam



param bp6
    caption = "Init Value 7th Entry"
    default = (0,0)
  endparam



param bp7
    caption = "Init Value 8th Entry"
    default = (0,0)
  endparam


param bp8
    caption = "Init Value 9th Entry"
    default = (0.1,0.1)
  endparam







  param plane
   caption = "Select visible matrix entry"
   default = 0
    enum = "a" "b" "c" "d" "e" "f" "g" "h" "i"
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc

 func fn3
    caption = "Pre function3"
    default = ident()
  endfunc




switch:
  type = "gwfa-c-MatrixPlain_3_Jul"
  seed = #pixel
  plane = @plane  
  bailout = @bailout
  bp1 = @bp1
  bp2 = @bp2
  bp3 = @bp3
  bp4 = @bp4
  bp5 = @bp5
  bp6 = @bp6
  bp7 = @bp7
  bp8 = @bp8
  fn1 = fn1
  fn2 = fn2
  fn3 = fn3
}



gwfa-c-MatrixPlain_3_Jul {
; Iteration of z^2 in the algebra of complex 3x3 matrices
; bailout following the matrix norm Tr(aa')
; matrix notation like [[a,b,c],[d,eh,f], [g,h,i]]
init:

complex  z = pixel 

if @plane == 0
complex a = z
complex b = @bp1
complex c = @bp2
complex d = @bp3
complex eh = @bp4
complex f = @bp5
complex g = @bp6
complex h = @bp7
complex i = @bp8




elseif @plane == 1
complex b = z
complex c = @bp1
complex d = @bp2
complex eh = @bp3
complex f = @bp4
complex g = @bp5
complex h = @bp6
complex i = @bp7
complex a = @bp8


elseif @plane == 2
complex c = z
complex d = @bp1
complex eh = @bp2
complex f = @bp3
complex g = @bp4
complex h = @bp5
complex i = @bp6
complex a = @bp7
complex b = @bp8



elseif @plane == 3
complex d = z
complex eh = @bp1
complex f = @bp2
complex g = @bp3
complex h = @bp4
complex i = @bp5
complex a = @bp6
complex b = @bp7
complex c = @bp8



elseif @plane == 4
complex eh = z
complex f = @bp1
complex g = @bp2
complex h = @bp3
complex i = @bp4
complex a = @bp5
complex b = @bp6
complex c = @bp7
complex d = @bp8



elseif @plane == 5
complex f = z
complex g = @bp1
complex h = @bp2
complex i = @bp3
complex a = @bp4
complex b = @bp5
complex c = @bp6
complex d = @bp7
complex eh = @bp8



elseif @plane == 6
complex g = z
complex h = @bp1
complex i = @bp2
complex a = @bp3
complex b = @bp4
complex c = @bp5
complex d = @bp6
complex eh = @bp7
complex f = @bp8



elseif @plane == 7
complex h = z
complex i = @bp1
complex a = @bp2
complex b = @bp3
complex c = @bp4
complex d = @bp5
complex eh = @bp6
complex f = @bp7
complex g = @bp8


elseif @plane == 8
complex i = z
complex a = @bp1
complex b = @bp2
complex c = @bp3
complex d = @bp4
complex eh = @bp5
complex f = @bp6
complex g = @bp7
complex h = @bp8



endif



float r = real(@bailout)^2

loop:




complex anew = a^2+b*d+c*g + @seed

complex bnew = a*b+b*eh+c*h + @seed

complex cnew = a*c+b*f+c*i + @seed

complex dnew = a*d+d*eh+f*g + @seed

complex enew = b*d+eh^2+f*h + @seed

complex fnew = c*d+eh*f+f*i + @seed

complex gnew = a*g+d*h+g*i + @seed

complex hnew = b*g+eh*h+h*i + @seed

complex inew = c*g+f*h+i^2 + @seed









a = fn1(anew)
b = fn2(bnew)
c = fn3(cnew)
d = fn1(dnew)

eh = fn2(enew)
f = fn3(fnew)
g = fn1(gnew)
h = fn2(hnew)

i = fn3(inew)


if @plane == 0
z  = a 
elseif @plane == 1
z  =  b
elseif @plane == 2
z  = c 
elseif @plane == 3
z  = d 
elseif @plane == 4
z  =  eh
elseif @plane == 5
z  = f 
elseif @plane == 6
z  = g 
elseif @plane == 7
z  =  h
elseif @plane == 8
z  = i
endif

float   rho = cabs(a)^2+cabs(b)^2+cabs(c)^2+cabs(d)^2
            + cabs(eh)^2+cabs(f)^2+cabs(g)^2+cabs(h)^2 + cabs(i)^2

bailout:
rho <= r

default:
  title = "Matrix Plain 3x3 Julia"
  maxiter = 1000
  periodicity = 2
  center = (0,0)
  magn = 0.5
  
  param bailout
    caption = "Bailout Value"
    default = 100
  endparam

  param bp1
    caption = "Init Value 2nd Entry"
    default = (0,0)
  endparam

param bp2
    caption = "Init Value 3rd Entry"
    default = (0,0)
  endparam

param bp3
    caption = "Init Value 4th Entry"
    default = (0,0)
  endparam


param bp4
    caption = "Init Value 5th Entry"
    default = (0.1,0.1)
  endparam


param bp5
    caption = "Init Value 6th Entry"
    default = (0,0)
  endparam



param bp6
    caption = "Init Value 7th Entry"
    default = (0,0)
  endparam



param bp7
    caption = "Init Value 8th Entry"
    default = (0,0)
  endparam


param bp8
    caption = "Init Value 9th Entry"
    default = (0.1,0.1)
  endparam



param seed
    caption = "Pertubation"
    default = (0,0)
  endparam






  param plane
   caption = "Select visible matrix entry"
   default = 0
    enum = "a" "b" "c" "d" "e" "f" "g" "h" "i"
  endparam
 
 func fn1
    caption = "Pre function1"
    default = ident()
  endfunc

 func fn2
    caption = "Pre function2"
    default = ident()
  endfunc

 func fn3
    caption = "Pre function3"
    default = ident()
  endfunc




switch:
  type = "gwfa-c-MatrixPlain_3"
  
  plane = @plane  
  bailout = @bailout
  bp1 = @bp1
  bp2 = @bp2
  bp3 = @bp3
  bp4 = @bp4
  bp5 = @bp5
  bp6 = @bp6
  bp7 = @bp7
  bp8 = @bp8
  fn1 = fn1
  fn2 = fn2
  fn3 = fn3
}