# 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;
}