Struggling to Understand Algorithm for Displaying Polynomial Matings with Julia Sets I'm working on creating a program that visualizes projected Julia Sets on a Riemann Sphere (such as my video here) when I came across this website visualizing matings between Julia Sets, and I want to recreate them for my own program (such as this video). However, with any resource that I've read that explains the process, I can't seem to wrap my mind around what's going on... I'm not sure if I simply don't yet have the formal education required (my knowledge of complex analysis is only limited to visualizing iterated fractals), or if these sources are just hard to understand.
What I want to learn specifically about is what is described here (from the previous website - what's in bold is what I want to learn, and what's italicized is what I have a hard time conceptually understanding):
"A progressive interpolation was introduced, between the two polynomial Julia sets and their mating. It consists in gluing equipotentials together and gives a holomorphic dynamical system between different spheres (this was observed by Milnor). This dynamical systems gives an easy method for drawing a conformally correct picture of the deformation of the polynomial Julia sets under the equipotential gluing: this method was explained to me by Buff. The result is an image which depends on the potential. This is what the movies show: the potential starts high and slowly approaches 0."
Essentially, what I'm looking for is given:

*

*some point z on the complex plane (I already know how to project this
onto the Riemann Sphere)

*two filled Julia Set coordinates $c_1$ and $c_2$ (for example, the Basilica and Rabbit - eventually I hope to move beyond two)

*some value t that represents the value of the potential that decreases to 0 (for the mating animation)

*some value n that represents the maximum escape-time iterations

*some value b that represents the bailout value

... do some math that calculates the color for that point (just like the escape-time algorithm - though this is the limit of my understanding, so I'm hoping that I can visualize the matings in the same way) when it's projected on the Riemann Sphere. Is this possible? I would be grateful for anything to help my understanding with this! If I'm in too far in over my head with this kind of math, then I'd also be satisfied with a copy-and-paste solution for my particular goal here.
I already tried reading these papers:

*

*Pasting Together Julia Sets: A Worked Out Example of Mating

*The Medusa Algorithm for Polynomial Matings

*The Thurston Algorithm for Quadratic Matings

*Slow mating and equipotential gluing

*Slow mating of quadratic Julia sets
I did consider putting this on StackOverflow instead, but I think this is more of a math question than a programming one.
EDIT:
After a week of going through Claude's code, I finally figured out an algorithm to which I can display the slow mating in real time! Its implementation in my project is not without a couple of bugs, but I was able to get the basic animation working (I've made some videos to show the mating of Basilica vs. Rabbit, its inverse, and its projection on the Riemann Sphere). The algorithm is as follows:
INITIALIZATION
Constants
R1 >= 5
R2 = R1 * R1
R4 = R2 * R2

Variables
# the two Julia Sets to slow mate
Complex p
Complex q

# mating presets
int mating_iterations
int intermediate_steps

# Julia Set presets
int julia_iterations
float bailout

# image presets
int width
int height
    
# intermediate path segments
Complex x [mating_iterations * intermediate_steps]
Complex y [mating_iterations * intermediate_steps]

# store inverse of pullback function (https://mathr.co.uk/blog/2020-01-16_slow_mating_of_quadratic_julia_sets.html)
Complex ma [mating_iterations * intermediate_steps]
Complex mb [mating_iterations * intermediate_steps]
Complex mc [mating_iterations * intermediate_steps]
Complex md [mating_iterations * intermediate_steps]

# what's sent to the GPU
Complex ma_frame [mating_iterations];
Complex mb_frame [mating_iterations];
Complex mc_frame [mating_iterations];
Complex md_frame [mating_iterations];

# Compute potentials and potential radii
float t[intermediate_steps]
float R[intermediate_steps]

for s: the count of intermediate segments
{
    t[s] = (s + .5) / intermediate_steps
    
    R[s] = exp(pow(2, 1 - t[s]) * log(R1))
}


p_i = 0     # nth iteration of the p Julia Set
q_i = 0     # nth iteration of the q Julia Set

# Calculate path arrays (Wolf Jung's equations 20 and 21)
for i: each frame in mating_iterations*intermediate_steps
{
    # i = intermediate_steps * n + s
    # for each n:
    #     for each s
    int s = i % intermediate_steps;
    int n = (i - s) / intermediate_steps;    # this is not needed here
    

    # Equation 20
           1 + ((1 - t[s]) * q / R2)                p_i / R[s]
    x[i] = ------------------------- * -------------------------------------
           1 + ((1 - t[s]) * p / R2)   1 + ((1 - t[s]) * q / R4 * (p_i - p))

    # Alternatively, if R1 = 1e10
    x[i] = p_i / R[s]



    # Equation 21
           1 + (1 - t[s]) * q / R2   R[s]
    y[i] = ----------------------- * ---- * (1 + ((1 - t[s]) * p / R4 * (q_i - q)))
           1 + (1 - t[s]) * p / R2   q_i

    # Alternatively, if R1 = 1e10
    y[i] = R[s] / q_i
          
    

    if (s == intermediate_steps - 1)    # last 's' before new 'n'
    {
        p_i = p_i^2 + p
        q_i = q_i^2 + q
    }
}

Prior to point calculation (CPU Render Loop)
# This could've be done using a nested for loop, but I needed to be consistent with my notation so I could understand the algorithm easier

for i: each frame in mating_iterations*intermediate_steps
{
    # i = intermediate_steps * n + s
    # for each n:
    #     for each s
    int s = i % intermediate_steps;
    int n = (i- s) / intermediate_steps;
        
    int first = intermediate_steps + s
    int s_prev = (s + intermediate_steps - 1) % intermediate_steps
        
    if (n > 0)
    {
        // Pull back x and y (Wolf Jung's Equation 22)
        for k: count of total mating iterations - current mating iteration (n)
        {
            int k_next = k + 1
            int next = intermediate_steps * k_next + s
            int prev = intermediate_steps * k + s_prev


                         (  1 - y[first]     x[next] - x[first]  )
            z_x[k] = sqrt(  ------------  *  ------------------  )
                         (  1 - x[first]     x[next] - y[first]  )
        
        
                                                                   
                                                     x[first]
                                                 1 - --------
                         (  (1 - y[first])           y[next]   )
            z_y[k] = sqrt(  --------------  *  --------------  )
                         (  (1 - x[first])           y[first]  )
                                                 1 - --------                                                                 
                                                     y[next]
        
            // choose sign by continuity
            if (length(-z_x[k] - x[prev]) < length(z_x[k] - x[prev]))
            {
                z_x[k] = -z_x[k]
            }
            if (length(-z_y[k] - y[prev]) < length(z_y[k] - y[prev]))
            {
                z_y[k] = -z_y[k]
            }
        }
        
        // copy results into path arrays
        for k: count of total mating iterations - current iteration (n)
        {
            x[intermediate_steps * k + s] = z_x[k]
            y[intermediate_steps * k + s] = z_y[k]
        }
    }
    
    a = x[intermediate_steps + s]
    b = y[intermediate_steps + s]
    ma[i] = b * (1 - a)
    mb[i] = a * (b - 1)
    mc[i] = 1 - a
    md[i] = b - 1
    
    for k: 0 to current mating iteration (n)
    {
        ma_frame[k] = ma[intermediate_steps * k + s]
        mb_frame[k] = mb[intermediate_steps * k + s]
        mc_frame[k] = mc[intermediate_steps * k + s]
        md_frame[k] = md[intermediate_steps * k + s]
    }

    # SEND VARIABLES TO GPU
        julia_iterations
        bailout
        p
        q
        R (taken from 'R[s]')
        current_mating_iteration (taken from 'n')
        ma_frame
        mb_frame
        mc_frame
        md_frame
}

Apply for each point on the complex plane (GPU Fragment Shader: for each pixel on the screen)
z = point on complex plane

for k: starting from current_mating_iteration and decreasing to zero
{
        ma_frame[k] * z + mb_frame[k]
    z = -----------------------------
        mc_frame[k] * z + md_frame[k]
}
    
if (length(z) < 1)
{
    c = p
    w = R * z
}
else
{
    c = q
    w = R / z    # note: this is complex division
}


for i: the rest of the regular Julia Set iterations (julia_iterations - n)
{
    break if (length(z) > bailout)
    
    w = w^2 + c
}

pixel_color = based on w

 A: From what you are saying I am unfortunately almost certain that you do not have enough mathematical background to really understand the theory of matings (which is something quite advanced, not realistically understandable before grad school). You typically want to have taken an advanced course in complex analysis, topology, and some background in complex dynamics.
It is difficult to explain shortly the notion of matings, but I'll try to give a quick summary of my own understanding of the topic. An equipotential is a line where the Green function is constant and strictly positive (if you don't know what the Green function is, you need background in complex dynamics). If the Julia set is connected then this equipotential is homeomorphic to a circle. In particular, it splits the sphere into two components, one containing the Julia set and the other containing $\infty$, and each of those components are homeomorphic to a disk. If you take two different connected Julia sets and two such equipotentials, you can glue the components containing the Julia sets by identifying the two equipotentials (if you don't understand this sentence you need background in topology). What you get then is something homeomorphic to a sphere (let's call it $S_1$, and you get a continuous map $f$ defined on $S_1$ that coincides with the restrictions of both polynomials outside of the glueing line. However the range of $f$ is not $S_1$ but rather a similar object obtained by glueing two different equipotentials (the ones that are images of the previous ones by the two polynomials). So you get a continuous map $f: S_1 \to S_2$, where $S_1$ and $S_2$ are topological spaces homeomorphic to spheres.
With some work and a really deep theorem that I won't even try to explain here, called the measurable Riemann mapping theorem, you can somehow get a holomorphic map $g: \hat{\mathbb C} \to \hat{ \mathbb C}$ from this whole business. The map $g$ is conjugated to $f$ by homeomorphisms that map the $S_i$ to the Riemann sphere. However, you shouldn't consider that $g$ is a dynamical system since it is conjugated to $f$, which has different domain and range. But if you used equipotentials $G=t$ for $S_1$, then you used equipotentials $G=dt$ for $S_2$, and as $t \to 0$ the difference between these two equipotentials shrinks to zero. So you want to prove that the map $g_t$ that you get with this whole procedure has a limit when $t \to 0$. This is not true in general, but when it is the case that limit is exactly what's called the mating (in one sense) between the two polynomials.
A: The theory is also beyond my level of education, but a practical implementation to make pictures is relatively straight forward following Chapter 5 of Wolf Jung's paper "The Thurston Algorithm for Quadratic Matings" (preprint linked in question).  However, an important thing missing in my code is detecting homotopy violations, so there is no proof that the images are correct.
In my code, the slow mating is computed as per Wolf Jung's chapter 5, pulling back the critical orbits using continuity to choose the sign of the square root.  Pulling back an orbit means the next orbit's $z_n$ depends in some way on the previous orbit's $z_{n+1}$.  The process has a sequence of orbits, where the orbit at time $t+s+1$ depends on the orbits at time $t + s$ (for choosing roots by continuity) and time $t + 1$ (for the input to the square root function).  Increasing the granularity $s$ presumably makes the continuity test more reliable.
To render images, the pull-back in Wolf Jung's paper is inverted to a sequence of functions of the form $z\to\frac{az^2+b}{cz^2+d}$, which are composed in reverse order starting from the desired pixel coordinates. Then choose the hemisphere based on $|z|<1$ or $|z|>1$, find $w=Rz$ or $w=R/z$ and $c=c_1$ or $c=c_2$ depending on hemisphere, and continue to iterate $w→w^2+c$ until escape (or maximum iteration count is reached).
Here is a scrappy diagram of the process that I made, which is how I initially understood how it works.  The top left triangle (green) is computed for the critical orbits, with the aim to calculate the coefficients of the inverse of the bottom diagonal.  Then the red path is computed per pixel.  The subdiagram at the right shows the continuity checking process.

For colouring filaments with distance estimation I use dual complex numbers for automatic differentiation, for colouring interior I use a function of the final $w$ to adjust hue.  To keep images stable for animations, the same number of total iterations is needed for the interior pixels each frame.
