Hi Adam,

Please clarify what you need. Is it the Mandelbrot code for Forward differencing?

Maybe this code is what you're looking for...

`/*`

FWDDIFF.CPP a module for determining the slope using forward differencing calculations of fractals.

Written in Microsoft Visual 'C++' by Paul de Leeuw.

https://github.com/hrkalona/Fractal-Zoomer/blob/master/src/fractalzoomer/core/ThreadDraw.java#L3640

https://github.com/hrkalona/Fractal-Zoomer/blob/master/src/fractalzoomer/main/app_settings/BumpMapSettings.java

*/

#include "slope.h"

/**************************************************************************

Calculate functions

**************************************************************************/

void CSlope::DoSlopeFwdDiffFn(Complex *z, Complex *q)

{

Complex sqr, temp;

Complex t = { 1.0, 0.0 };

double real_imag;

int k;

int degree = (int)param[5];

int degree1, degree2;

// PaletteStart = (int)param[4];

switch (subtype)

{

case 0: // Mandelbrot

sqr.x = z->x * z->x;

sqr.y = z->y * z->y;

real_imag = z->x * z->y;

z->x = q->x + sqr.x - sqr.y;

z->y = q->y + real_imag + real_imag;

break;

// other cases removed

}

}

/**************************************************************************

Get the gradients in the x and y directions

**************************************************************************/

double CSlope::getGradientX(double *wpixels, int index, int width)

{

int x = index % width;

double it = *(wpixels + index);

if (x == 0) {

return (*(wpixels + index + 1) - it) * 2;

}

else if (x == width - 1) {

return (it - *(wpixels + index - 1)) * 2;

}

else {

double diffL = it - *(wpixels + index - 1);

double diffR = it - *(wpixels + index + 1);

return diffL * diffR >= 0 ? 0 : diffL - diffR;

}

}

double CSlope::getGradientY(double *wpixels, int index, int width, int height)

{

int y = index / width;

double it = *(wpixels + index);

if (y == 0) {

return (it - *(wpixels + index + width)) * 2;

}

else if (y == height - 1) {

return (*(wpixels + index - width) - it) * 2;

}

else {

double diffU = it - *(wpixels + index - width);

double diffD = it - *(wpixels + index + width);

return diffD * diffU >= 0 ? 0 : diffD - diffU;

}

}

/**************************************************************************

Brightness Scaling

**************************************************************************/

int CSlope::changeBrightnessOfColorScaling(int rgb, double delta)

{

int new_color = 0;

double bump_transfer_factor = param[0]; // future add to dialogue box

// double mul = getBumpCoef(delta);

double mul = (1.5 / (fabs(delta * bump_transfer_factor) + 1.5));

if (delta > 0) {

rgb ^= 0xFFFFFF;

int r = rgb & 0xFF0000;

int g = rgb & 0x00FF00;

int b = rgb & 0x0000FF;

int ret = (int)(r * mul + 0.5) & 0xFF0000 | (int)(g * mul + 0.5) & 0x00FF00 | (int)(b * mul + 0.5);

new_color = 0xff000000 | (ret ^ 0xFFFFFF);

}

else {

int r = rgb & 0xFF0000;

int g = rgb & 0x00FF00;

int b = rgb & 0x0000FF;

new_color = 0xff000000 | (int)(r * mul + 0.5) & 0xFF0000 | (int)(g * mul + 0.5) & 0x00FF00 | (int)(b * mul + 0.5);

}

return new_color;

}

/**************************************************************************

Slope Fractal

**************************************************************************/

int CSlope::RunSlopeFwdDiff(HWND hwndIn, void(*plot)(WORD, WORD, DWORD), int user_data(HWND hwnd), char* StatusBarInfo, bool *ThreadComplete, int subtypeIn, int NumThreadsIn, int threadIn, Complex j, double mandel_width, double hor, double vert, double xgap, double ygap,

double rqlim, long threshold, double paramIn[], CTrueCol *TrueCol, CDib *Dib, double *SlopeDataPtr, BYTE juliaflag)

{

Complex z = 0.0; // initial value for iteration Z0

Complex c, q;

BigComplex cBig, zBig, qBig;

BigDouble BigTemp, BigTemp1;

double log_zn, nu;

int lastChecked = -1;

int x, y, i;

DWORD index;

double iterations;

thread = threadIn;

NumThreads = NumThreadsIn;

subtype = subtypeIn;

hwnd = hwndIn;

for (i = 0; i < NUMSLOPEDERIVPARAM; i++)

param[i] = paramIn[i];

if (NumThreads == 0)

StripWidth = xdots;

else

StripWidth = xdots / NumThreads;

StripStart = StripWidth * thread;

*ThreadComplete = false;

for (y = 0; y < ydots; y++)

{

if (BigNumFlag)

cBig.y = BigVert + BigWidth - Big_ygap * y;

else

c.y = vert + mandel_width - y * ygap;

double progress = (double)y / ydots;

if (int(progress * 100) != lastChecked)

{

lastChecked = int(progress * 10);

sprintf(StatusBarInfo, "Progess (%d%%), %d Threads", int(progress * 100), NumThreads);

}

for (x = StripStart; x < StripStart + StripWidth; x++)

{

if (user_data(hwnd) < 0) // trap user input

return -1;

c.x = hor + x * xgap;

{

if (juliaflag)

{

q = j;

z = c;

}

else

{

q = c;

z = 0.0;

}

}

iterations = 0.0;

for (;;)

{

iterations++;

if (iterations >= threshold)

break;

DoSlopeFwdDiffFn(&z, &q);

if ((z.x * z.x + z.y * z.y) >= rqlim)

break;

}

if (iterations < threshold)

{

if (BigNumFlag)

{

// sqrt of inner term removed using log simplification rules.

BigTemp = zBig.x * zBig.x + zBig.y * zBig.y;

BigTemp1 = BigTemp.BigLog();

log_zn = mpfr_get_d(BigTemp1.x, MPFR_RNDN) / 2;

nu = log(log_zn / log(2.0)) / log(2.0);

}

else

{

// sqrt of inner term removed using log simplification rules.

log_zn = log(z.x * z.x + z.y * z.y) / 2;

nu = log(log_zn / log(2.0)) / log(2.0);

}

iterations = iterations + 1 - nu;

}

index = ((DWORD)y * (DWORD)width) + (DWORD)x;

if (x >= 0 && x < xdots - 1 && y >= 0 && y < ydots - 1)

*(SlopeDataPtr + index) = iterations;

plot(x, y, (long)iterations);

}

}

*ThreadComplete = true;

return 0;

}

/**************************************************************************

Slope Fractal

**************************************************************************/

int CSlope::RenderSlope(long threshold, double paramIn[], CTrueCol *TrueCol, CDib *Dib, double *SlopeDataPtr)

{

double dotp, gradAbs, gradCorr, cosAngle, sizeCorr, smoothGrad, lightAngleRadians, lightx, lighty;

for (int i = 0; i < NUMPARAM; i++)

param[i] = paramIn[i];

double bumpMappingDepth = param[3]; // future add to dialogue box

double bumpMappingStrength = param[4]; // future add to dialogue box

double lightDirectionDegrees = param[2]; // future add to dialogue box

int PaletteStart = (int)param[1];

double iterations;

int lastChecked = -1;

DWORD index;

int x, y;

double gradx, grady;

unsigned char r, g, b;

RGBTRIPLE colour;

int modified;

lastChecked = -1;

sizeCorr = 0.0;

lightx = 0.0;

lighty = 0.0;

gradCorr = pow(2, (bumpMappingStrength - DEFAULT_BUMP_MAPPING_STRENGTH) * 0.05);

sizeCorr = ydots / pow(2, (MAX_BUMP_MAPPING_DEPTH - bumpMappingDepth) * 0.16);

lightAngleRadians = lightDirectionDegrees * PI / 180.0;

lightx = cos(lightAngleRadians) * gradCorr;

lighty = sin(lightAngleRadians) * gradCorr;

for (y = 0; y < ydots; y++)

{

for (x = 0; x < xdots; x++)

{

index = ((DWORD)y * (DWORD)width) + (DWORD)x;

if (SlopeDataPtr)

iterations = *(SlopeDataPtr + index);

else

return 0; // oops, we don't actually have forward differencing

if (iterations >= threshold)

{ // interior of Mandelbrot set = inside_color

colour.rgbtRed = (BYTE)TrueCol->InsideRed; // M_waves

colour.rgbtGreen = (BYTE)TrueCol->InsideGreen;

colour.rgbtBlue = (BYTE)TrueCol->InsideBlue;

}

else

{

// modified = rgbs[index];

if (iterations < PaletteStart)

modified = 0x00FFFFFF;

else

modified = 0xFF000000 | ((DWORD)*(TrueCol->PalettePtr + (BYTE)(((long)iterations) % 256) * 3 + 0) << 16) | ((DWORD)*(TrueCol->PalettePtr + (BYTE)(((long)iterations) % 256) * 3 + 1) << 8) | *(TrueCol->PalettePtr + (BYTE)(((long)iterations) % 256) * 3 + 2);

gradx = getGradientX(SlopeDataPtr, index, xdots);

grady = getGradientY(SlopeDataPtr, index, xdots, ydots);

dotp = gradx * lightx + grady * lighty;

// int original_color = modified; // not sure what this is for

if (dotp != 0)

{

gradAbs = sqrt(gradx * gradx + grady * grady);

cosAngle = dotp / gradAbs;

smoothGrad = -2.3562 / (gradAbs * sizeCorr + 1.5) + 1.57;

//smoothGrad = Math.atan(gradAbs * sizeCorr);

modified = changeBrightnessOfColorScaling(modified, cosAngle * smoothGrad);

}

r = (modified >> 16) & 0xFF;

g = (modified >> 8) & 0xFF;

b = modified & 0xFF;

colour.rgbtRed = r;

colour.rgbtGreen = g;

colour.rgbtBlue = b;

}

outRGBpoint(Dib, x, y, colour);

}

}

return 0;

}

Let me know if you need more info.