reb
Class InvertMobius

Object
  extended by common:Generic
      extended by reb:InvertSphere
          extended by reb:InvertMobius

class 
InvertSphere:InvertMobius

Derived from the public C code of Curtis McMullen and material from Indra's Pearls by Mumford, Series and Wright A new method for InitSpheres is needed, as both spheres and Kleinian Group matricies must be initalized.


Ultra Fractal Source

Toggle UF Source Code Display

 class InvertMobius(InvertSphere) {
 ; <p>
 ;  Derived from the public C code of Curtis McMullen and
 ;  material from Indra's Pearls by Mumford, Series and Wright
 ;  A new method for InitSpheres is needed, as both spheres
 ;  and Kleinian Group matricies must be initalized.
 ;
 public:
   import "common.ulb"
 
 ;  constructor for Mobius (Kleinian) sphere inversions
   func InvertMobius(Generic pparent)
     invertSphere.Invertsphere(0)
     PARAFUZZ  = 1e-3, LINEFUZZ  = 1e-5, im = (0,1),gi = 0
     fii = 4, fks = 1, fk = fks, scale = 0.7, sscale = 1, iend = 4
 
     ss = new SphereArray(50000000)
 
     InitSpheresK()
   endfunc
 
 ; Initialize spheres, inversion matricies and trace values
   func InitSpheresK()
 ; This function initializes the base sphere(s), the Kleinian group
 ; inversion matrices and the trace values
     complex Fc[4]
     float Fz[4]
     float Fr[4]
     int i = 0
     complex ta = @ta
     complex tb = @tb
     float cir = @cir
     complex cira = @cira
     Complex cirs = @cirs
     float rada = @rada
     float rads = @rads
     float radam = @radam
     complex m = @m
     complex tab = 0
     complex gz0 = 0
     complex gz = 0
     complex gQ = 0
     complex gR = 0
     complex tabAB = @tabAB
     Fz[0] = 0
     Fz[1] = 0
     Fz[2] = 0
     Fz[3] = 0
 
     if @cusptype == "Nearby group"
       ta =  (1.91,-0.05)
       tb =  (2,0)
       cir = 0.297
     elseif @cusptype == "Accident"
       ta = (1.9021130325903071442328786667588,0)
       tb = (1.9021130325903071442328786667588,0)
     elseif @cusptype == "Troel's Point"
       ta = (1.6168866453,-0.7056734968)
       tb = (2,0)
       cir = 0.13091
     elseif @cusptype == "Degenerate Spiral"
       ta = (1.9264340533,-0.0273817919)
       tb = (2,0)
       cir = 0.2419
     elseif @cusptype == "Double Cusp"
       tb =  (2,0)
       if @cusp == "0/1"
         ta =  (2,0)
         cir = 0.5
       elseif  @cusp == "1/100"
         ta =  (1.99599050389291,-0.000302826761310194)
         cir = 0.3091
       elseif  @cusp == "2/99"
         ta =  (1.99599050389291,-0.000302826761310194)
         cir = 0.2676
       elseif  @cusp == "3/100"
         ta =  (1.99118310579765,-0.000982081178197968)
         cir = 0.2427
       elseif  @cusp == "4/99"
         ta =  (1.98427416931522,-0.00237401043014847)
         cir = 0.2351
       elseif  @cusp == "1/15"
         ta =  (1.95859103011179,-0.0112785606117653)
         cir = 0.32870768
       elseif @cusp == "1/10"
         ta =  (1.91342329586682,-0.0362880775225429)
         cir = 0.32365169
       elseif @cusp == "2/19"
         ta =  (1.90377999779718,-0.0395799512688805)
         cir = 0.2604321
       elseif @cusp == "1/9"
         ta =  (1.89640725094921,-0.0487530128986506)
         cir = 0.32176507
       elseif @cusp == "7/43"
         ta =  (1.80785523999043,-0.136998687907137)
         cir = 0.229796
       elseif @cusp == "3/10"
         ta =  (1.65831239517773,-0.5)
         cir = 0.2316625
       elseif @cusp == "2/5"
         ta =  (1.64213876865348,-0.766588417465459)
         cir = 0.2665884
       elseif @cusp == "1/2"
         ta =  (1.73205080756888,-1)
         cir = 0.36602539
       elseif @cusp == "21/34"
         ta =  (1.6179907967521,-1.29170028664218)
         cir = 0.166461
       endif
     endif
 
     ; Grandma's Special algorithm
 
     if !(@cusptype == "Grandma's Special #2" || @cusptype == "Maskit")
       if @showcusp == "Both"
         Fc[0] = 1-cir
         Fr[0] = cir
         Fc[1] = -(1-cir)
         Fr[1] = cir
         Fc[2] = (0,0)
         Fr[2] = 1
         s.sph[0] = new Sphere(Fc[0],Fz[0],Fr[0],0,0)
         s.sph[1] = new Sphere(Fc[1],Fz[1],Fr[1],0,0)
         s.sph[2] = new Sphere(Fc[2],Fz[2],Fr[2],0,0)
         ss.sph[0] = s.sph[0]
         ss.sph[1] = s.sph[1]
         ss.sph[2] = s.sph[2]
       elseif @showcusp == "Alternate"
         Fc[0] = (1-cir)
         Fr[0] = cir
         Fc[1] = -(1-cir)
         Fr[1] = cir
         s.sph[0] = new Sphere(Fc[0],Fz[0],Fr[0],0,0)
         s.sph[1] = new Sphere(Fc[1],Fz[1],Fr[1],0,0)
         ss.sph[0] = s.sph[0]
         ss.sph[1] = s.sph[1]
       elseif @showcusp == "Standard"
         Fc[0] = (0,0)
         Fr[0] = 1
         s.sph[0] = new Sphere(Fc[0],Fz[0],Fr[0],0,0)
         ss.sph[0] = s.sph[0]
       endif
       complex tab = (ta*tb+sqrt(ta*ta*tb*tb-4*(ta*ta+tb*tb)))/2
       if @alttab
         tab = (ta*tb-sqrt(ta*ta*tb*tb-4*(ta*ta+tb*tb)))/2
       endif
       complex gz0 = (tab-2)*tb/(tb*tab-2*ta+2*(0,1)*tab)
       mob[0] = new Mobius((tb-(0,2))/2, tb/2, tb/2, conj((tb-(0,2))/2))
       mob[1] = new Mobius(mob[0].m_d, -mob[0].m_b, -mob[0].m_c, mob[0].m_a)
       mob[2] = new Mobius(ta/2, (ta*tab-2*tb-(0,4))*gz0/(2*tab-4), \
                 (ta*tab-2*tb+(0,4))/((2*tab+4)*gz0), ta/2)
       if @trans
         temp1 = mob[2].m_b
         mob[2].m_b = mob[2].m_c
         mob[2].m_c = temp1
       endif
       mob[3] = new Mobius(mob[2].m_d, -mob[2].m_b, -mob[2].m_c, mob[2].m_a)
     endif
 
     ; Figure 11.1 from the Indra's Pearls book
     if @cusptype == "Indra 11.1"
       tabAB = (0,0)
       tb = (2,0)
       sscale = 2.4
       ta = (1.924781, 0.047529)
       cirs = (0, 0.272591599)
       rads = 0.41421381
       cira = (0.168208, 0)
       rada = 0.2211131
     endif
 
     ; Grandma's Special #2 algorithm
 
     if !(@cusptype == "Grandma's Special" || @cusptype == "Maskit" ||\
          @cusptype == "Double Cusp" || @cusptype == "Nearby Group" ||\
          @cusptype == "Troel's Point"|| @cusptype == "Degenerate Spiral")
       if @showcusp == "Both"
         Fc[0] = cirs
         Fr[0] = rads-cabs(cirs)
         if @lines
           Fr[0] = -Fr[0]
         endif
         Fc[1] = -cirs
         Fr[1] = Fr[0]
         Fc[2] = cira
         Fr[2] = rada-cabs(cira)
         if @linea
           Fr[2] = -Fr[2]
         endif
         Fc[3] = -cira
         Fr[3] = Fr[2]
         s.sph[0] = new Sphere(Fc[0],Fz[0],Fr[0],0,0)
         s.sph[1] = new Sphere(Fc[1],Fz[1],Fr[1],0,0)
         s.sph[2] = new Sphere(Fc[2],Fz[2],Fr[2],0,0)
         s.sph[3] = new Sphere(Fc[3],Fz[3],Fr[3],0,0)
         ss.sph[0] = s.sph[0]
         ss.sph[1] = s.sph[1]
         ss.sph[2] = s.sph[2]
         ss.sph[3] = s.sph[3]
       elseif @showcusp == "Alternate"
         Fc[0] = cira
         Fr[0] = rada-cabs(cira)
         if @linea
           Fr[0] = -Fr[0]
         endif
         Fc[1] = -cira
         Fr[1] = Fr[0]
         s.sph[0] = new Sphere(Fc[0],Fz[0],Fr[0],0,0)
         s.sph[1] = new Sphere(Fc[1],Fz[1],Fr[1],0,0)
         ss.sph[0] = s.sph[0]
         ss.sph[1] = s.sph[1]
       else
         Fc[0] = cirs
         Fr[0] = rads-cabs(cirs)
         if @lines
           Fr[0] = -Fr[0]
         endif
         Fc[1] = -cirs
         Fr[1] = Fr[0]
         s.sph[0] = new Sphere(Fc[0],Fz[0],Fr[0],0,0)
         s.sph[1] = new Sphere(Fc[1],Fz[1],Fr[1],0,0)
         ss.sph[0] = s.sph[0]
         ss.sph[1] = s.sph[1]
       endif
       gz = 0.5*sqrt(ta^2*tb^2-4*ta^2-4*tb^2+4*tabAB+8)
       tab = 0.5*ta*tb-gz
       if @alttab
         tab = 0.5*ta*tb+gz
       endif
       gQ = sqrt(2-tabAB)
       gR = sqrt(2+tabAB)
       if cabs(tabAB+im*gQ*gR)<2
         gR = -gR
       endif
       gz0 = (tab-2)*(tb+gR)/(tb*tab-2*ta+im*gQ*tab)
       mob[0] = new Mobius((tb-im*gQ)/2, (tb*tab-2*ta-im*gQ*tab)/((2*tab+4)*gz0),\
                   (tb*tab-2*ta+im*gQ*tab)*gz0/(2*tab-4), (tb+im*gQ)/2)
       mob[1] = new Mobius(mob[0].m_d, -mob[0].m_b, -mob[0].m_c, mob[0].m_a)
       mob[2] = new Mobius(ta/2, (ta*tab-2*tb+2*im*gQ)/((2*tab+4)*gz0), \
                 (ta*tab-2*tb-2*im*gQ)*gz0/(2*tab-4), ta/2)
       if @trans
         temp1 = mob[2].m_b
         mob[2].m_b = mob[2].m_c
         mob[2].m_c = temp1
       endif
       mob[3] = new Mobius(mob[2].m_d, -mob[2].m_b, -mob[2].m_c, mob[2].m_a)
     endif
 
     if @cusptype == "Maskit"
       if @showcuspm == "Predefined"
         if @cuspm == "0/1"
           m =  (0,2)
           radam = 0.5
         elseif @cuspm == "1/15"
           m = (-0.0112785606117653, 1.95859103011179)
           radam = 0.489661311
         elseif @cuspm == "1/10"
           m =  (-0.0362880775225429,1.91342329586682)
           radam = 0.47852888
         elseif @cuspm == "2/19"
           m =  (-0.0395799512688805,1.90377999779718)
           radam = 0.35214005
         elseif @cuspm == "1/9"
           m =  (-0.0487530128986506,1.89640725094921)
           radam = 0.4744165
         elseif @cuspm == "7/43"
           m =  (-0.136998687907137,1.80785523999043)
           radam = 0.2983515
         elseif @cuspm == "3/10"
           m =  (-0.5,1.65831239517773)
           radam = 0.3015111
         elseif @cuspm == "2/5"
           m =  (-0.766588417465459,1.64213876865348)
           radam = 0.363491
         elseif @cuspm == "1/2"
           m =  (-1,1.73205080756888)
           radam = 0.57735015
         elseif @cuspm == "21/34"
           m =  (1.29170028664218,1.6179907967521)
           radam = 0.1997052
         endif
       endif
       if @showcuspm == "Both" || (@showcuspm == "Predefined" \
              && @showcuspm2 == "Both")
         Fc[0] = 1 + flip(radam)
         Fr[0] = radam
         Fc[1] = (1,0)
         Fr[1] = 0
         s.sph[0] = new Sphere(Fc[0],Fz[0],Fr[0],0,0)
         s.sph[1] = new Sphere(Fc[1],Fz[1],Fr[1],0,0)
         ss.sph[0] = s.sph[0]
         ss.sph[1] = s.sph[1]
       elseif @showcuspm == "Alternate" || (@showcuspm == "Predefined" \
              && @showcuspm2 == "Alternate")
         Fc[0] = 1 + flip(radam)
         Fr[0] = radam
         s.sph[0] = new Sphere(Fc[0],Fz[0],Fr[0],0,0)
         ss.sph[0] = s.sph[0]
       else
         Fc[0] = (1,0)
         Fr[0] = 0
         s.sph[0] = new Sphere(Fc[0],Fz[0],Fr[0],0,0)
         ss.sph[0] = s.sph[0]
       endif
 
       mob[0] = new Mobius(m, 1, 1, 0)
       mob[1] = new Mobius(mob[0].m_d, -mob[0].m_b, -mob[0].m_c, mob[0].m_a)
       mob[2] = new Mobius(1, 0, 2, 1)
       mob[3] = new Mobius(mob[2].m_d, -mob[2].m_b, -mob[2].m_c, mob[2].m_a)
     endif
 
 ; normalize the matrices
     i = 0
     while i < 4
       mob[i].Normalize(mob[i])
       i = i + 1
     endwhile
   endfunc
 
 ; Inverts a sphere - Mobius inversion The return value may be a sphere or a plane. <br>
 ; @param target = sphere to be inverted
 ; @param inverter = not used.
 ; @param level = recursion level
 ; @param gen = index of generator matrix
 ;------------------------------------------------
    Sphere func Inverse(Sphere target, Sphere inverter, int level, int gen)
      mob[5] = new Mobius(0,0,0,0)
      mob[5].CircToMob(mob[5],target)
      mob[gen].MConj(mob[5],mob[5])
      Sphere ns = new Sphere(0,0,0,0,0)
      mob[5].MobToCirc(level,gen,ns)
       return ns
    endfunc
 
  ;  Top level recursion function
  func Recurse()
    int l = 0
    int ir = 0
    int si = 0
    int sj = 0
    int j = 0
    bool done = false
    int c = 0
    float ddx = 0
    float ddy = 0
    float ddr = 0
    sphere rra
    bool continue = false
    int ifinal = 0
    int istack = 0
    int nold = 0
    int nfinal = 1
    int iii = 0
    int k = 0
    int lev = 0
    if @showcusp == "Both"
      j = 3
      iii = 2
    elseif @showcusp == "Standard"
      j = 1
      iii = 0
    else
      j = 2
      iii = 1
    endif
    if @cusptype == "Grandma's Special #2"
      if @showcusp == "Both"
        j = 4
        iii = 3
      else
        j = 2
        iii = 1
      endif
    endif
    if @cusptype == "Maskit"
      if @showcuspm == "Both" || (@showcuspm == "Predefined" && \
          @showcuspm2 == "Both")
        j = 2
        iii = 1
      else
        j = 1
        iii = 0
      endif
    endif
    iii = iii + 1
    lev = 1
    while lev <= mlevel && j < flength
      gi = 0
      while gi < iend  && j < flength
        k = 0
        while k <= iii  && j < flength
          ss.sph[j] = Inverse(ss.sph[k],ss.sph[0],lev,gi)
 
          if ss.sph[j].frad >= 0 && abs(ss.sph[j].frad) < fminsphere
            j = j - 1
          endif
          if j > flength
            lev = mlevel + 1
            k = iii + 1
          endif
          j = j+1
          k = k + 1
        endwhile
        if j > flength
          gi = 5
        endif
        gi = gi + 1
      endwhile
 
   ; sort array by circle identity
         l = round(j/2)+1
         ir = j
         continue = true
         repeat
           if l > 1
             l = l-1
             rra = ss.sph[l-1]
           else
             rra = ss.sph[ir-1]
             ss.sph[ir-1] = ss.sph[0]
             ir = ir-1
             if ir == 0
               ss.sph[0] = rra
               continue = false
             endif
           endif
 
           if continue == true
             si = l
             sj = 2*l
           endif
 
           while (sj <= ir) && (continue == true)
             if sj < ir
               ddx = real(ss.sph[sj-1].fcen)-real(ss.sph[sj].fcen)
               ddy = imag(ss.sph[sj-1].fcen)-imag(ss.sph[sj].fcen)
               ddr = ss.sph[sj-1].frad-ss.sph[sj].frad
               c = 0
               done = false
 
           ; Lines precede circles
             if ss.sph[sj-1].frad <= 0 && ss.sph[sj].frad > 0
                 c = -1
                 done = true
             elseif ss.sph[sj].frad <= 0 && ss.sph[sj-1].frad > 0
               c = 1
               done = true
 
               elseif ddx < -thresh
                 c = -1
                 done = true
               elseif ddx > thresh && !done
                 c = 1
                 done = true
               elseif ddy < -thresh && !done
                 c = -1
                 done = true
               elseif ddy > thresh && !done
                 c = 1
                 done = true
               elseif ss.sph[sj-1].frad <= 0 && !done     ; Lines
                 if ddr < -PARAFUZZ && !done
                   c = -1
                   done = true
                 elseif ddr > PARAFUZZ && !done
                   c = 1
                   done = true
                 endif
               elseif !done                       ; Circles
                 if ddr < -thresh && !done
                   c = -1
                   done = true
                 elseif ddr > thresh && !done
                   c = 1
                   done = true
                 endif
               endif
               if c < 0
                 sj = sj + 1
               endif
             endif
             ddx = real(rra.fcen)-real(ss.sph[sj-1].fcen)
             ddy = imag(rra.fcen)-imag(ss.sph[sj-1].fcen)
             ddr = rra.frad-ss.sph[sj-1].frad
             c = 0
             done = false
 
   ;        ; Lines precede circles
            if rra.frad <= 0 && ss.sph[sj-1].frad > 0
              c = -1
              done = true
            elseif ss.sph[sj-1].frad <= 0 && rra.frad > 0
              c = 1
              done = true
            elseif ddx < -thresh
              c = -1
              done = true
            elseif ddx > thresh && !done
              c = 1
              done = true
            elseif ddy < -thresh && !done
              c = -1
              done = true
            elseif ddy > thresh && !done
              c = 1
              done = true
            elseif rra.frad <= 0 && !done           ; Lines
              if ddr < -PARAFUZZ && !done
                c = -1
                done = true
              elseif ddr > PARAFUZZ && !done
                c = 1
                done = true
              endif
            elseif !done                       ; Circles
              if ddr < -thresh && !done
                c = -1
                done = true
              elseif ddr > thresh && !done
                c = 1
                done = true
              endif
            endif
            if c < 0
              ss.sph[si-1] = ss.sph[sj-1]
              si = sj
              sj = sj + sj
            else
              sj = ir + 1
            endif
          endwhile
          if (continue == true)
            ss.sph[si-1] = rra
          endif
        until continue == false
 
     ; eliminate duplicates in place
 
 ;          ir = oldj-iii+1
 ;          l = oldj-iii
           ir = j+1
           l = j
           while ir < j
             ddx = real(ss.sph[ir].fcen)-real(ss.sph[ir-1].fcen)
             ddy = imag(ss.sph[ir].fcen)-imag(ss.sph[ir-1].fcen)
             ddr = ss.sph[ir].frad-ss.sph[ir-1].frad
             c = 0
             done = false
 
           ; Lines precede circles
             if ss.sph[ir].frad <= 0 && ss.sph[ir-1].frad > 0
               c = -1
               done = true
             elseif ss.sph[ir-1].frad <= 0 && ss.sph[ir].frad > 0
               c = 1
               done = true
             elseif ddx < -thresh && !done
               c = -1
               done = true
             elseif ddx > thresh && !done
               c = 1
               done = true
             elseif ddy < -thresh && !done
               c = -1
               done = true
             elseif ddy > thresh && !done
               c = 1
               done = true
             elseif ss.sph[ir].frad <= 0 && !done           ; Lines
               if ddr < -PARAFUZZ && !done
                 c = -1
                 done = true
               elseif ddr > PARAFUZZ && !done
                 c = 1
                 done = true
               endif
            elseif !done                       ; Circles
               if ddr < -thresh && !done
                 c = -1
                 done = true
               elseif ddr > thresh && !done
                 c = 1
                 done = true
               endif
             endif
             if c == 0
               ir = ir + 1
             else
               ss.sph[l] = ss.sph[ir]
               ir = ir +1
               l = l + 1
             endif
           endwhile
           j = j-(ir-l)+1
 
     ; move to final array, skipping dups
         ifinal = 0
         istack = 0
         nold = nfinal
         while ifinal < nold && istack < j && nfinal < flength
           ddx = real(s.sph[ifinal].fcen)-real(ss.sph[istack].fcen)
           ddy = imag(s.sph[ifinal].fcen)-imag(ss.sph[istack].fcen)
           ddr = s.sph[ifinal].frad-ss.sph[istack].frad
           c = 0
           done = false
 
         ; Lines precede circles
           if s.sph[ifinal].frad <= 0 && ss.sph[istack].frad > 0
             c = -1
             done = true
           elseif ss.sph[istack].frad <= 0 && s.sph[ifinal].frad > 0
             c = 1
             done = true
           elseif ddx < -thresh
             c = -1
             done = true
           elseif ddx > thresh && !done
             c = 1
             done = true
           elseif ddy < -thresh && !done
             c = -1
             done = true
           elseif ddy > thresh && !done
             c = 1
             done = true
           elseif s.sph[ifinal].frad <= 0 && !done           ; Lines
             if ddr < -PARAFUZZ && !done
               c = -1
               done = true
             elseif ddr > PARAFUZZ && !done
               c = 1
               done = true
             endif
           elseif !done                       ; Circles
             if ddr < -thresh && !done
               c = -1
               done = true
             elseif ddr > thresh && !done
               c = 1
               done = true
             endif
           endif
           if c < 0
             ifinal = ifinal + 1
           endif
           if c == 0
             istack = istack + 1
           endif
           if c > 0
             s.sph[nfinal] = ss.sph[istack]
             nfinal = nfinal + 1
             istack = istack + 1
           endif
         endwhile
         while istack < j  && nfinal < flength
           s.sph[nfinal] = ss.sph[istack]
           nfinal = nfinal + 1
           istack = istack + 1
         endwhile
         ir = nold
         while ir <nfinal
           ss.sph[ir-nold] = s.sph[ir]
           ir = ir+1
         endwhile
         j = nfinal-nold
 
  ; sort final array by circle identity
         l = round((nfinal)/2)+1
         ir = nfinal
         continue = true
         repeat
           if l > 1
             l = l-1
             rra = s.sph[l-1]
           else
             rra = s.sph[ir-1]
             s.sph[ir-1] = s.sph[0]
             ir = ir-1
             if ir == 0
               s.sph[0] = rra
               continue = false
             endif
           endif
 
           if continue == true
             si = l
             sj = 2*l
           endif
 
           while (sj <= ir) && (continue == true)
             if sj < ir
               ddx = real(s.sph[sj-1].fcen)-real(s.sph[sj].fcen)
               ddy = imag(s.sph[sj-1].fcen)-imag(s.sph[sj].fcen)
               ddr = s.sph[sj-1].frad-s.sph[sj].frad
               c = 0
               done = false
 
             ; Lines precede circles
               if s.sph[sj-1].frad <= 0 && s.sph[sj].frad > 0
                 c = -1
                 done = true
               elseif s.sph[sj].frad <= 0 && s.sph[sj-1].frad > 0
                 c = 1
                 done = true
               elseif ddx < -thresh
                 c = -1
                 done = true
               elseif ddx > thresh && !done
                 c = 1
                 done = true
               elseif ddy < -thresh && !done
                 c = -1
                 done = true
               elseif ddy > thresh && !done
                 c = 1
                 done = true
               elseif s.sph[sj-1].frad <= 0 && !done           ; Lines
                 if ddr < -PARAFUZZ && !done
                   c = -1
                   done = true
                 elseif ddr > PARAFUZZ && !done
                   c = 1
                   done = true
                 endif
               elseif !done                       ; Circles
                 if ddr < -thresh && !done
                   c = -1
                   done = true
                 elseif ddr > thresh && !done
                   c = 1
                   done = true
                 endif
               endif
               if c < 0
                 sj = sj + 1
               endif
             endif
             ddx = real(rra.fcen)-real(s.sph[sj-1].fcen)
             ddy = imag(rra.fcen)-imag(s.sph[sj-1].fcen)
             ddr = rra.frad-s.sph[sj-1].frad
             c = 0
             done = false
 
           ; Lines precede circles
             if rra.frad <= 0 && s.sph[sj-1].frad > 0
               c = -1
               done = true
             elseif s.sph[sj-1].frad <= 0 && rra.frad > 0
               c = 1
               done = true
             elseif ddx < -thresh
               c = -1
               done = true
             elseif ddx > thresh && !done
               c = 1
               done = true
             elseif ddy < -thresh && !done
               c = -1
               done = true
             elseif ddy > thresh && !done
               c = 1
               done = true
             elseif rra.frad <= 0 && !done           ; Lines
               if ddr < -PARAFUZZ && !done
                 c = -1
                 done = true
               elseif ddr > PARAFUZZ && !done
                 c = 1
               done = true
               endif
             elseif !done                       ; Circles
               if ddr < -thresh && !done
                 c = -1
                 done = true
               elseif ddr > thresh && !done
                 c = 1
                 done = true
               endif
             endif
             if c < 0
               s.sph[si-1] = s.sph[sj-1]
               si = sj
               sj = sj + sj
             else
               sj = ir + 1
             endif
           endwhile
           if (continue == true)
             s.sph[si-1] = rra
           endif
         until continue == false
     iii = j-1
     if j == 0
       lev = mlevel+1
     else
       lev = lev+1
     endif
   endwhile
 
   fk = nfinal-1
   ss = 0
   heapsort(s,fk)
   vanquish(s,fk)
   SizeAndRot()
   if @zsort
     heapsortz(s, fk)
   endif
  endfunc
 
 ; Size and rotate spheres
  func SizeAndRot()
    if @cusptype == "Nearby group" && @trans
        sscale = 0.05
    elseif @cusptype == "Double Cusp" && @trans
      if @cusp == "1/15"
        sscale = 0.015
      elseif @cusp == "0/1"
        sscale = 0.5
      elseif @cusp == "2/19"
        sscale = 0.025
      elseif @cusp == "1/9"
        sscale = 0.04
      elseif @cusp == "7/43"
        sscale = 0.035
      elseif @cusp == "2/5"
        sscale = 0.2
      endif
    endif
    int i = 0
    while i < fk
      if s.sph[i].frad > 0
        ; size routines
        if s.sph[i].frad > @maxsphere
          s.sph[i].frad = 0
        endif
        if @cusptype == "Maskit"
          s.sph[i].fcen = (s.sph[i].fcen+(0,-1))*#magn*scale*sscale
        else
          s.sph[i].fcen = s.sph[i].fcen*#magn*scale*sscale
        endif
        s.sph[i].fz = s.sph[i].fz*#magn*scale*sscale
        s.sph[i].frad = s.sph[i].frad*#magn*scale*sscale
 
   ;         rotation around the z axis
   ;
        float fx = real(s.sph[i].fcen)
        float fy = imag(s.sph[i].fcen)
        float fz = s.sph[i].fz
        float xx = real(s.sph[i].fcen)
        fx = fx*cos(#angle) - fy*sin(#angle)
        fy = fy*cos(#angle) + xx*sin(#angle)
        s.sph[i].fcen = fx + flip(fy)
 
   ;         rotation around the y axis
   ;
        xx = fx
        fx = fz*sin(@yang*#pi/180) + fx*cos(@yang*#pi/180)
        fz = fz*cos(@yang*#pi/180) - xx*sin(@yang*#pi/180)
        s.sph[i].fcen = fx + flip(fy)
        s.sph[i].fz = fz
 
   ;         rotation around the x axis
   ;
        float yy = fy
        fy = fy*cos(@xang*#pi/180) - fz*sin(@xang*#pi/180)
        fz = yy*sin(@xang*#pi/180) + fz*cos(@xang*#pi/180)
        s.sph[i].fcen = fx + flip(fy)
        s.sph[i].fz = fz
      endif
 
      i = i + 1
    endwhile
  endfunc
 
   float PARAFUZZ
   float LINEFUZZ
   complex im
   Mobius mob[6]
   int gi
   int iend
   float sscale
   spherearray ss
 
 default:
   title = "Mobius Inversions"
   int param v_mobinvert
     caption = "Version (Mobius Inversions)"
     default = 101
     hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used."
     visible = @v_mobinvert < 101
   endparam
 
   float param minsphere
     caption = "Min sphere size"
     default = 0.001
     min = 0.000001
   endparam
   float param maxsphere
     caption = "Max sphere size"
     default = 0.7
   endparam
   int param maxLevel
     caption = "Recursion level"
     default = 100
     min = 0
     max = 200
   endparam
   float param radadj
     caption = "Gen Rad adj"
     default = 1.0
     visible = false
   endparam
   float param bradadj
     caption = "Base Rad adj"
     default = 1.0
     visible = false
   endparam
   bool param reflect
     caption = "Reflect Object"
     default = false
   endparam
   float param xang
     caption = "X Axis Rotation"
     default = 0.0
   endparam
   float param yang
     caption = "Y Axis Rotation"
     default = 0.0
   endparam
   heading
     text = "For rotation around the Z Axis use the rotation angle on the Location Tab."
   endheading
   heading
     text = "Use these methods with care, as only certain combinations of \
             parameters will give satisfactory results. Use of the help file \
             is essential. Use of the 'Explore' feature is not recommended."
     visible = @cusptype == "Grandma's Special" || @cusptype == "Grandma's Special #2" \
               || @cusptype == "Maskit"
   endheading
   heading
     text = "The Grandma algorithms come from 'Indra's Pearls' by Mumford, \
             Series and Wright. Consult the book for use of the trace parameters."
     visible = @cusptype == "Grandma's Special" || @cusptype == "Grandma's Special #2"
   endheading
   heading
     text = "The Maskit algorithm comes from 'Indra's Pearls' by Mumford, \
             Series and Wright. Consult the book for use of the Maskit parameters."
     visible = @cusptype == "Maskit"
   endheading
   param cusptype
     caption = "Group"
     default = 7
     enum = "Accident" "Indra 11.1" "Degenerate Spiral" "Double Cusp" \
            "Grandma's Special" "Grandma's Special #2" "Maskit" "Nearby Group" \
            "Troel's Point"
   endparam
   param cusp
     caption = "Cusp"
     default = 5
     enum = "0/1" "1/100" "2/99" "3/100" "4/99" "1/15" "1/10" "2/19" "1/9" \
            "7/43" "3/10" "2/5" "1/2" "21/34"
     visible = @cusptype == "Double Cusp"
   endparam
   param showcusp
     caption = "Cusp view"
     default = 1
     enum = "Both" "Standard" "Alternate"
     visible = !(@cusptype == "Maskit"  || @cusptype == "Accident")
   endparam
   param showcuspm
     caption = "Cusp view"
     default = 1
     enum = "Both" "Standard" "Alternate" "Predefined"
     visible = @cusptype == "Maskit"
   endparam
   param showcuspm2
     caption = "Cusp view"
     default = 1
     enum = "Both" "Standard" "Alternate"
     visible = @cusptype == "Maskit" && @showcuspm == "Predefined"
   endparam
 
   heading
     text = "Slowly increase 'gen circle' until good circle packing is observed."
     visible = @cusptype == "Grandma's special"&& @showcusp != "Standard"
   endheading
   param @cir
     caption ="Gen circle"
     default = 0.01
     hint = "Slowly increase until good circle packing is observed."
     visible = @cusptype == "Grandma's special"&& @showcusp != "Standard"
   endparam
   heading
     text = "These parameters set the size and position of the 'Standard' \
             generator circles."
     visible = @cusptype == "Grandma's special #2"&& @showcusp != "Alternate"
   endheading
   param @cirs
     caption ="Gen circle std"
     default = (0,0)
     visible = (@cusptype == "Grandma's special #2"&& @showcusp != "Alternate")
   endparam
   param @rads
     caption ="Rad ref std"
     default = 1.0
     visible = (@cusptype == "Grandma's special #2"&& @showcusp != "Alternate")
   endparam
   bool param @lines
     caption ="Generate as line"
     default = false
     visible = (@cusptype == "Grandma's special #2"&& @showcusp != "Alternate")
   endparam
   heading
     text = "These parameters set the size and position of the 'Alternate' \
             generator circles."
     visible = @cusptype == "Grandma's special #2"&& @showcusp != "Standard"
   endheading
   param @cira
     caption ="Gen circle alt"
     default = (0,0)
     visible = (@cusptype == "Grandma's special #2"&& @showcusp != "Standard")
   endparam
   param @rada
     caption ="Rad ref alt"
     default = 1.0
     visible = (@cusptype == "Grandma's special #2"&& @showcusp != "Standard")
   endparam
   heading
     text = "This parameter set the size and position of the 'Alternate' \
             generator circle."
     visible = @cusptype == "Maskit"&& @showcusp != "Standard"
   endheading
   param @radam
     caption ="Rad ref"
     default = 0.5
     visible = @cusptype == "Maskit"&& @showcuspm != "Standard" && \
               @showcuspm != "Predefined"
   endparam
   bool param @linea
     caption ="Generate as line"
     default = false
     visible = (@cusptype == "Grandma's special #2"&& @showcusp != "Standard")
   endparam
   complex param m
     caption = "Maskit param"
     default = (0.0112785606117653,1.95859103011179)
     visible = @cusptype == "Maskit" && @showcuspm != "Predefined"
   endparam
   param cuspm
     caption = "Cusp"
     default = 1
     enum = "0/1" "1/15" "1/10" "2/19" "1/9" "7/43" "3/10" "2/5" "1/2" "21/34"
     visible = @cusptype == "Maskit" && @showcuspm == "Predefined"
   endparam
   complex param ta
     caption = "Trace a"
     default = (1.95,0.02)
     visible = @cusptype == "Grandma's Special" || @cusptype == "Grandma's Special #2"
   endparam
   complex param tb
     caption = "Trace b"
     default = 3
     visible = @cusptype == "Grandma's Special" || @cusptype == "Grandma's Special #2"
   endparam
   bool param alttab
     caption ="Alternate trace ab"
     default = false
     visible = @cusptype == "Grandma's Special" || @cusptype == "Grandma's Special #2"
   endparam
   complex param tabAB
     caption = "Trace abAB"
     default = -2
     visible = @cusptype == "Grandma's Special #2"
   endparam
   bool param trans
     caption = "Transpose transform"
     default = false
     visible = @cusptype != "Maskit"
   endparam
 }
 


Constructor Summary
InvertMobius()
           
InvertMobius(Generic pparent)
          constructor for Mobius (Kleinian) sphere inversions
 
Method Summary
 void InitSpheresK()
          Initialize spheres, inversion matricies and trace values
 Sphere Inverse(Sphere target, Sphere inverter, int level, int gen)
          Inverts a sphere - Mobius inversion The return value may be a sphere or a plane.
 void Recurse()
          Top level recursion function
 void SizeAndRot()
          Size and rotate spheres
 
Methods inherited from class reb:InvertSphere
findbound, heapify, heapifyz, heapsort, heapsortz, InitSpheres, recurse2, siftdown, siftdownz, swap, vanquish
 
Methods inherited from class common:Generic
GetParent
 
Methods inherited from class Object
 

Constructor Detail

InvertMobius

public InvertMobius(Generic pparent)
constructor for Mobius (Kleinian) sphere inversions


InvertMobius

public InvertMobius()
Method Detail

InitSpheresK

public void InitSpheresK()
Initialize spheres, inversion matricies and trace values


Inverse

public Sphere Inverse(Sphere target,
                      Sphere inverter,
                      int level,
                      int gen)
Inverts a sphere - Mobius inversion The return value may be a sphere or a plane.

Overrides:
Inverse in class InvertSphere
Parameters:
target - = sphere to be inverted
inverter - = not used.
level - = recursion level
gen - = index of generator matrix

Recurse

public void Recurse()
Top level recursion function

Overrides:
Recurse in class InvertSphere

SizeAndRot

public void SizeAndRot()
Size and rotate spheres