Mandelbrot set, transformed into a teardrop shape, in R and Rcpp

mandeldrop.cpp

#include <Rcpp.h>
using namespace Rcpp;

// Mandelbrot algorithm based on: https://stackoverflow.com/q/48069990
// Smooth colouring and AA algorithm from iq: https://iquilezles.org/articles/msetsmooth/
// [[Rcpp::export]]
Rcpp::NumericMatrix mandeldrop(
    const double x_min,
    const double x_max,
    const double y_min,
    const double y_max,
    const int res_x,
    const int res_y,
    const int nb_iter
) {
  Rcpp::NumericMatrix ret(res_x, res_y);
  double x_step = (x_max - x_min) / res_x;
  double y_step = (y_max - y_min) / res_y;
  int r,c;
  int AA = 2;
  
  for( int m=0; m < AA; m++ ) {
    for( int n=0; n < AA; n++ ) {
      for (r = 0; r < res_y; r++) {
        for (c = 0; c < res_x; c++) {
          double zx = 0.0, zy = 0.0, new_zx;
          double rcx = x_min + c*x_step + m/AA, rcy = y_min + r*y_step + n/AA;
          double cx = 0.2*sinh(rcy)/(cosh(rcy)-cos(rcx));
          double cy = 0.2*sin(rcx)/(cosh(rcy)-cos(rcx));
          int n = 0;
          for (n=0; (zx*zx + zy*zy < 4.0 ) && ( n < nb_iter ); n++ ) {
            new_zx = zx*zx - zy*zy + cx;
            zy = 2.0*zx*zy + cy;
            zx = new_zx;
          }
          ret(c,r) += n - log2(log2((zx*zx+zy*zy))) + 4.0;
        }
      }
    }
  }
  
  return ret/(AA*AA);
}

mandeldrop.R

Rcpp::sourceCpp("mandeldrop.cpp")

xlims=c(-0.65,0.65)
ylims=c(-0.6,2.2)
m <- mandeldrop(xlims[[1]], xlims[[2]], ylims[[1]], ylims[[2]], 1560, 2160, 512)

# Colour palette from https://stackoverflow.com/q/48069990
rainbow=c(rgb(0.47,0.11,0.53),rgb(0.27,0.18,0.73),rgb(0.25,0.39,0.81),rgb(0.30,0.57,0.75),rgb(0.39,0.67,0.60),rgb(0.51,0.73,0.44),rgb(0.67,0.74,0.32),rgb(0.81,0.71,0.26),rgb(0.89,0.60,0.22),rgb(0.89,0.39,0.18),rgb(0.86,0.13,0.13))
cols=c(colorRampPalette(rainbow)(100),rev(colorRampPalette(rainbow)(100)),"black")

par(bg="black")
image(
  log(m),
  col=cols,
  asp=diff(ylims)/diff(xlims),
  axes=F,
  useRaster=T
)