/*
 * MandelVroom 2.0
 *
 * (c) Copyright 1987,1989  Kevin L. Clague, San Jose, CA
 *
 * All rights reserved.
 *
 * Permission is hereby granted to distribute this program's source
 * executable, and documentation for non-comercial purposes, so long as the
 * copyright notices are not removed from the sources, executable or
 * documentation.  This program may not be distributed for a profit without
 * the express written consent of the author Kevin L. Clague.
 *
 * This program is not in the public domain.
 *
 * Fred Fish is expressly granted permission to distribute this program's
 * source and executable as part of the "Fred Fish freely redistributable
 * Amiga software library."
 *
 * Permission is expressly granted for this program and it's source to be
 * distributed as part of the Amicus Amiga software disks, and the
 * First Amiga User Group's Hot Mix disks.
 *
 * contents:  this file contains the routines to calculate Mandelbrot and
 * Julia pictures in Amiga Fast Floating Point.
 */

#include "mandp.h"

extern SHORT MaxOrbit;

static LONG MaxIter;

/*
 * Fast Floating Point Mandelbrot Generator
 */
MandelbrotFFP( Pict, StartX, StartY, GapX, GapY )
  /*
   * These parameters are IEEE 32 bit floating point numbers
   */

  struct Picture *Pict;

  LONG  StartX, StartY, GapX, GapY;
{
  register int i, j, k;
  register SHORT *CountPtr;

  register float curx, cury;

  float startx, starty, gapx, gapy;
  float SPFieee();

  UBYTE MandFlag;

  /* Here we convert the IEEE parameters into FFP format */

  startx = SPFieee( StartX );
  starty = SPFieee( StartY );
  gapx   = SPFieee( GapX );
  gapy   = SPFieee( GapY );

  MaxIter = Pict->MaxIteration;

  if (Pict->Flags & NO_RAM_GENERATE)
    CountPtr = Pict->Counts;
  else
    CountPtr = Pict->Counts + Pict->CurLine*Pict->CountX;

  /* start in the upper left hand corner */

  cury = starty;

  /*
   * for each row
   */
  for (i = Pict->CurLine; i < Pict->CountY; i++) {

    curx = startx;

    MandFlag = 0;

    if ( Pict->Flags & NO_RAM_GENERATE )
      CountPtr = Pict->Counts;

    /*
     * for each collumn
     */

    for (j = 0; j < Pict->CountX; j++) {

     /* if we've hit a Mandelbrot point, then use the convergence detector,
      *  to reduce compute time
      */

      if (*CountPtr == 0) {

        if ( MandFlag ) {

          k = FFP_Trace_Height( 0.0, 0.0, curx, cury );

        } else {

          k = FFP_Height( 0.0, 0.0, curx, cury ); /* Normal calculator */
        }

        MandFlag = k == Pict->MaxIteration;

        *CountPtr = k;        /* save it in ram for recoloring */
      }
      CountPtr++;

      curx += gapx;

      ChildPause( Pict );
    }
    cury += gapy;

    CheckEOL( Pict );
  }
} /* Mandelbrot */

/***************************************************************************
 *
 *           Julia Curve Generator that uses Amiga Fast Floating Point
 *
 **************************************************************************/

/*
 * Fast Floating Point Julia Generator
 */
JuliaFFP( Pict, JuliaX, JuliaY, StartX, StartY, GapX, GapY )

  /*
   * These parameters are IEEE 32 bit floating point numbers
   */

  register struct Picture *Pict;
  LONG  JuliaX, JuliaY;
  LONG  StartX, StartY;
  LONG  GapX,   GapY;
{
  register int i, j, k;
  register float curx, cury;

  register SHORT *CountPtr;

  float juliax, juliay, gapx, gapy,startx;
  float SPFieee();

  UBYTE ConvFlag;

  MaxIter = Pict->MaxIteration;

  if (Pict->Flags & NO_RAM_GENERATE)
    CountPtr = Pict->Counts;
  else
    CountPtr = Pict->Counts + Pict->CurLine*Pict->CountX;

  /* Here we convert the IEEE parameters into FFP format */

  juliax = SPFieee( JuliaX );
  juliay = SPFieee( JuliaY );

  gapx = SPFieee( GapX );
  gapy = SPFieee( GapY );

  /* start in the upper left hand corner */

  startx = SPFieee(StartX);
  cury   = SPFieee(StartY);

  /* move down to next not generated line */

  cury += Pict->CurLine*gapy;

  /*
   * for each row
   */
  for (i = Pict->CurLine; i < Pict->CountY; i++) {

    curx = startx;
    ConvFlag = 0;

    if ( Pict->Flags & NO_RAM_GENERATE )
      CountPtr = Pict->Counts;

    /*
     * for each collumn
     */

    for (j = 0; j < Pict->CountX; j++) {

     /* if we've hit a Mandelbrot point, then use the convergence detector,
      *  to reduce compute time
      */

      if (*CountPtr == 0) {

        if (ConvFlag) {
                                           /* try the convergence method */
          k = FFP_Trace_Height( curx, cury, juliax, juliay);
        } else {
          k = FFP_Height( curx, cury, juliax, juliay);      /* Normal calc */
        }

        ConvFlag = k == Pict->MaxIteration;

        *CountPtr = k;        /* save it in ram for recoloring */
      }
      CountPtr++;

      curx += gapx;

      ChildPause( Pict );
    }
    cury += gapy;

    CheckEOL( Pict );
  }
} /* JuliaFFP */

/*
 * This function calculated the 'potential' of a given location
 * in the complex plane.
 */
int
FFP_Height( posx, posy, juliax, juliay )
  float posx;    /* Location in screen coordinate space         */
  float posy;    /* Location in screen coordinate space         */
  float juliax;  /* Real part of location on complex plane      */
  float juliay;  /* Imaginary part of location on complex plane */
{
  register float cura, curb;
  register float cura2, curb2;
  register LONG k;

#ifdef CHECK_TASK_STACK

  CheckStack();

#endif

  cura = cura2 = posx;
  curb = curb2 = posy;

  cura2 *= cura2;
  curb2 *= curb2;

  for (k = 0; k < MaxIter; k ++ ) {

    curb *= cura;                     /* b = 2 * a * b  + juliay    */
    curb += curb + juliay;

    cura = cura2 - curb2 + juliax;    /* a = a2 - b2 + juliax       */

    cura2 = cura * cura;              /* a2 = a * a                 */
    curb2 = curb * curb;              /* b2 = b * b                 */

    if (cura2+curb2 >= 16.0)          /* quit if a2+b2 is too big   */
      return( k );
  }
  return( k );
}

/*
 * this is also a Mandelbrot potential calculator, that is tailored to
 * identify points that are 'in' the Mandelbrot set.  Points outside the
 * Mandelbrot set diverge.  Points inside converge.
 * This code uses a trace table mechanism to detect convergence.
 */
int
FFP_Trace_Height( posx, posy, curx, cury )

  float posx;    /* Real part of location on complex plane      */
  float posy;    /* Imaginary part of location on complex plane */

  float curx;    /* Real part of location on complex plane      */
  float cury;    /* Imaginary part of location on complex plane */
{
  LONG  k, l;

  LONG TraceSize = 32;

  float olda[32], oldb[32]; /* Convergence trace table */

  register float cura, curb, cura2, curb2;
  register float *RealTrace, *ImagTrace;

#ifdef CHECK_TASK_STACK

  CheckStack();

#endif

  cura = cura2 = posx;
  curb = curb2 = posy;

  cura2 *= cura2;
  curb2 *= curb2;

  for (k = 0; k < MaxIter; k += TraceSize) {

    RealTrace = olda;
    ImagTrace = oldb;

    for (l = 0; l < TraceSize; l++) {

      *(RealTrace++) = cura;
      *(ImagTrace++) = curb;

      curb *= cura;
      curb += curb + cury;

      cura = cura2 - curb2 + curx;

      cura2 = cura * cura;
      curb2 = curb * curb;

      if (cura2+curb2 >= 16.0)
        return( k + l );
    }

    /* Scope out trace table for convergence */

    RealTrace = olda;
    ImagTrace = oldb;

    for (l = 0; l < TraceSize; l++)

      if (cura == *(RealTrace++) && curb == *(ImagTrace++)) {
        k += MaxIter;
        return( MaxIter );
      }
  }
  return( MaxIter );
}

DrawOrbitFFP( Pict, creal, cimag, zreal, zimag )
  register struct Picture *Pict;
  LONG zreal, zimag, creal, cimag;
{
  register struct RastPort *Rp;
  register float cura,  curb;
           float cura2, curb2;
           float Creal, Cimag;

  float x_scale, y_scale;
  int x_center, y_center;
  int width, height;
  int x, y;
  register int k;

  struct Window *Window;

  Window = OrbitWind;

  Rp = Window->RPort;

  width  = (Window->Width-Pict->LeftMarg-Pict->RightMarg);
  height = (Window->Height-Pict->TopMarg-Pict->BotMarg);

  x_center = width/2 + Pict->LeftMarg;
  y_center = height/2 + Pict->TopMarg;

  y_scale = x_scale = (float) height / 2.0;
/*x_scale *= AspectRatio( Pict );*/

  if ( Pict->pNode.ln_Type == MANDPICT ) {

    Creal = SPFieee(creal);
    Cimag = SPFieee(cimag);

    cura = cura2 = SPFieee(zreal);
    curb = curb2 = SPFieee(zimag);

  } else {

    Creal = SPFieee(zreal);
    Cimag = SPFieee(zimag);

    cura = cura2 = SPFieee(creal);
    curb = curb2 = SPFieee(cimag);

  }

  cura2 *= cura2;
  curb2 *= curb2;

  SetAPen( Rp, 0 );
  RectFill( Rp, Pict->LeftMarg, Pict->TopMarg, width+Pict->LeftMarg-1, height+Pict->TopMarg-1);

  SetAPen( Rp, HIGHLIGHTPEN);

  for (k = 0; k < MaxOrbit; k++ ) {

    curb *= cura;
    curb += curb + Cimag;

    cura  = cura2 - curb2 + Creal;

    cura2 = cura * cura;
    curb2 = curb * curb;

    if (cura2+curb2 >= 16.0)
      return( k );

    /* map real and imaginary parts into window coordinates */

    x = x_center + (int)(x_scale*cura);
    y = y_center + (int)(y_scale*curb);

    if ( x >= Pict->LeftMarg && x < Window->Width-Pict->RightMarg &&
         y >= Pict->TopMarg  && y < Window->Height-Pict->BotMarg ) {

      /* plot pixel location */
      WritePixel( Rp, x, y);
    }

  }
  return( k );
}