/*! \file bmark-teem.c
 *
 * \author John Reppy
 */

/*
 * COPYRIGHT (c) 2012 The Diderot Project (http://diderot-language.cs.uchicago.edu)
 * All rights reserved.
 */

// unu quantize -b 8 -i bmark-teem.nrrd -o bmark-teem.png 

#include <stdio.h>
#include <math.h>
#include <stdbool.h>
#include "teem/nrrd.h"
#include "teem/ell.h"
#include "teem-defs.h"

#if 1
#  define REAL  	float
#  define POW		powf
#  define SIN		sinf
#  define nrrdTypeReal	nrrdTypeFloat
#else
#  define REAL  	double
#  define POW		pow
#  define SIN		sin
#  define nrrdTypeReal	nrrdTypeDouble
#endif

#define RESOLUTION	2000
#define MAX_ITER	2000
#define FOV		0.15

STATIC_INLINE REAL lerp (REAL out_min, REAL out_max, REAL in_min, REAL in, REAL in_max)
{
    return (out_min * (in_max - in) + out_max * (in - in_min)) / (in_max - in_min);
}

STATIC_INLINE REAL sqr (REAL x) { return x*x; }

STATIC_INLINE REAL dot2 (REAL vec1[2], REAL vec2[2])
{
  return (vec1[0] * vec2[0] + vec1[1] * vec2[1]);
}

int main (int argc, const char *argv[])
{
    const char *me = argv[0];
    char *outFile = "bmark-teem.nrrd";
    REAL center[2] = {-1.2, -0.3};

    airArray *mop;
    mop = airMopNew();

    Nrrd *nOut = nrrdNew();
    airMopAdd(mop, nOut, (airMopper)nrrdNuke, airMopAlways);  
    if (nrrdAlloc_va(nOut, nrrdTypeReal, 3, 
		     AIR_CAST(size_t, 3),
		     AIR_CAST(size_t, RESOLUTION),
		     AIR_CAST(size_t, RESOLUTION)))
    {
	char *err = biffGetDone(NRRD);
	airMopAdd (mop, err, airFree, airMopAlways);
	fprintf (stderr, "%s: Trouble allocating:\n%s", me, err);
	airMopError (mop);
	return -1;
    }
    REAL *outData = AIR_CAST(REAL *, nOut->data);

    double t0 = airTime(); // start timing

    for (int ciIdx = 1; ciIdx <= RESOLUTION;  ciIdx++) {
	for (int crIdx = 1;  crIdx <= RESOLUTION;  crIdx++) {
	    REAL c[2] = {
		lerp(center[0]-FOV, center[0]+FOV, 1.0, (REAL)crIdx, (REAL)RESOLUTION),
		lerp(center[1]-FOV, center[1]+FOV, 1.0, (REAL)ciIdx, (REAL)RESOLUTION)
	      };
	    REAL z[2] = {0.0, 0.0};
	    int iter = 0;
	    while (true) {
	      // z = z^2 + c
		REAL re = sqr(z[0]) - sqr(z[1]);
		REAL im = 2.0 * z[0] * z[1];
		ELL_2V_SET(z, re + c[0], im + c[1]);
		iter += 1;
		if (dot2(z, z) > 4.0) {
		  // point has escaped; set color based on iteration
		    REAL t = 11.0* POW((REAL)iter, 0.2);
		    REAL red = (SIN(t) + 1.7) / 2.7;
		    REAL green = (SIN(t + 2.0*M_PI/3.0) + 1.7) / 2.7;
		    REAL blue = (SIN(t - 2.0*M_PI/3.0) + 1.7) / 2.7;
		    ELL_3V_SET(outData + 3*(crIdx-1 + RESOLUTION*(ciIdx-1)), red, green, blue);
		    break;
		}
		else if (iter > MAX_ITER) {
		    ELL_3V_SET(outData + 3*(crIdx-1 + RESOLUTION*(ciIdx-1)), 0.0, 0.0, 0.0);
		    break;
		}
	    }
	}
    }

    double totalTime = airTime() - t0;
    printf("usr=%f\n", totalTime);

    if (nrrdSave(outFile, nOut, NULL)) {
	char *err = biffGetDone(NRRD);
	airMopAdd(mop, err, airFree, airMopAlways);
	fprintf(stderr, "%s: Trouble saving:\n%s", me, err);
	airMopError(mop);
	return -1;
    }

    airMopOkay(mop);
    return 0;

}