# # External Ray In

Tomoki Kawahira “An algorithm to draw external rays of the Mandelbrot set” (2009) describes an algorithm to trace external rays.

The next point $$r$$ along an external ray with current doubled angle $$\theta$$ measured in turns, current depth $$p$$ and current radius $$R$$ satisfies:

$F^p(0,r)=\lambda R e^{2 \pi i \theta}$

where $$\lambda < 1$$ controls the sharpness of the ray. Applying Newton’s method in one complex variable:

$r_{m+1} = r_m - \frac{F^p(0,r_m) - \lambda R e^{2 \pi i \theta}} {\frac{\partial}{\partial c}F^p(0,r_m)}$

When crossing dwell bands, double $$\theta$$ and increment $$p$$, resetting the radius $$R$$. Stop tracing when close to the target (for example when within the basin of attraction for Newton’s method for nucleus).

## # 1 C99 Code

#include <complex.h>
#include <math.h>
#include <gmp.h>

struct m_exray_in {
mpq_t angle;
mpq_t one;
int sharpness;
double er;
double _Complex c;
int j;
int k;
};

struct m_exray_in *m_exray_in_new
(const mpq_t angle, int sharpness)
{
m_d_exray_in *ray = malloc(sizeof(*ray));
mpq_init(ray->angle);
mpq_set(ray->angle, angle);
mpq_init(ray->one);
mpq_set_ui(ray->one, 1, 1);
ray->sharpness = sharpness;
ray->er = 65536.0;
double a = twopi * mpq_get_d(ray->angle);
ray->c = ray->er * cexp(I * a);
ray->k = 0;
ray->j = 0;
return ray;
}

void m_exray_in_delete(struct m_exray_in *ray)
{
mpq_clear(ray->angle);
mpq_clear(ray->one);
free(ray);
}

double _Complex m_exray_in_step
(struct m_exray_in *ray, int n)
{
if (ray->j >= ray->sharpness)
{
mpq_mul_2exp(ray->angle, ray->angle, 1);
if (mpq_cmp_ui(ray->angle, 1, 1) >= 0)
mpq_sub(ray->angle, ray->angle, ray->one);
ray->k = ray->k + 1;
ray->j = 0;
}
double r = pow(ray->er,
pow(0.5, (ray->j + 0.5) / ray->sharpness));
double a = twopi * mpq_get_d(ray->angle);
double _Complex target = r * cexp(I * a);
double _Complex c = ray->c;
for (int m = 0; m < n; ++m)
{
double _Complex z = 0;
double _Complex dc = 0;
for (int p = 0; p <= ray->k; ++p)
{
dc = 2 * z * dc + 1;
z = z * z + c;
}
c = c - (z - target) / dc;
}
ray->j = ray->j + 1;
ray->c = c;
return c;
}