/* entheh */

#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "main.h"
#include "ripples.h"



/* Macro to calculate two-dimensional offsets in ripple tiles. */
#define P(x, y) ((y)*RIPPLE_D + (x))



RIPPLE_INFO *create_ripple_info() {

 RIPPLE_INFO *ri = malloc(sizeof(RIPPLE_INFO));

 if (!ri) return NULL;

 ri->frame = 0;
 ri->swapped = 0;
 ri->dist = 0;

 return ri;
}



void destroy_ripple_info(RIPPLE_INFO *ripple_info) {

 R_DISTURBANCE *rd = ripple_info->dist, *next;

 while (rd) {
  next = rd->next;
  destroy_disturbance(rd);
  rd = next;
 }

 free(ripple_info);
}



RIPPLES *create_ripples() {

 static int frame = 0;

 RIPPLES *r = malloc(sizeof(RIPPLES));

 if (!r) return NULL;

 memset(r->prev, 0, sizeof(r->prev));
 memset(r->curr, 0, sizeof(r->curr));
 r->frame = frame;

 frame = (frame + 1) & (RIPPLE_FRAMES - 1);

 return r;
}

void destroy_ripples(RIPPLES *r) {

	if (r)
		free(r);

	return;
}



void propagate_ripples(RIPPLE_INFO *ripple_info,
                       RIPPLES *ripples,
                       RIPPLES *rip_l,
                       RIPPLES *rip_r,
                       RIPPLES *rip_t,
                       RIPPLES *rip_b) {

 int x, y;

 ripple_p *p;
 ripple_p *c, *cl, *cr, *ct, *cb;

 if (ripples == 0 || ripples->frame != ripple_info->frame) return;

 if (ripple_info->swapped) {
  p = ripples->curr;
  c = ripples->prev;
  cl = rip_l ? rip_l->prev : 0;
  cr = rip_r ? rip_r->prev : 0;
  ct = rip_t ? rip_t->prev : 0;
  cb = rip_b ? rip_b->prev : 0;
 } else {
  p = ripples->prev;
  c = ripples->curr;
  cl = rip_l ? rip_l->curr : 0;
  cr = rip_r ? rip_r->curr : 0;
  ct = rip_t ? rip_t->curr : 0;
  cb = rip_b ? rip_b->curr : 0;
 }

 if (cl) cl += P(RIPPLE_D-1,0);
 if (ct) ct += P(0,RIPPLE_D-1);

 /*
 Here are the individual steps, excluding damping:

 accel = ((c[P(-1,0)] + c[P(1,0)] +
           c[P(0,-1)] + c[P(0,1)]) >> 2) - *c;

 velocity = *c - *p;

 *p = *c + velocity + accel;

 These have been optimised into a single expression.
 The extra line is for damping.
 */
 #define PROPAGATE_RIPPLE_P(l, r, t, b) {                         \
  register ripple_p n = *c - *p + (((l) + (r) + (t) + (b)) >> 2); \
  *p = n - (n >> 5);                                              \
 }

 /* Top Left */
 PROPAGATE_RIPPLE_P(cl ? *cl : *c, c[P(1,0)], ct ? *ct++ : *c, c[P(0,1)]);
 p++;
 c++;
 if (cl) cl += P(0,1);

 /* Top Row */
 for (x = 1; x < RIPPLE_D - 1; x++) {
  PROPAGATE_RIPPLE_P(c[P(-1,0)], c[P(1,0)], ct ? *ct++ : *c, c[P(0,1)]);
  p++;
  c++;
 }

 /* Top Right */
 PROPAGATE_RIPPLE_P(c[P(-1,0)], cr ? *cr : *c, ct ? *ct++ : *c, c[P(0,1)]);
 p++;
 c++;
 if (cr) cr += P(0,1);


 for (y = 1; y < RIPPLE_D - 1; y++) {

  /* Left */
  PROPAGATE_RIPPLE_P(cl ? *cl : *c, c[P(1,0)], c[P(0,-1)], c[P(0,1)]);
  p++;
  c++;
  if (cl) cl += P(0,1);

  /* Middle */
  for (x = 1; x < RIPPLE_D - 1; x++) {
   PROPAGATE_RIPPLE_P(c[P(-1,0)], c[P(1,0)], c[P(0,-1)], c[P(0,1)]);
   p++;
   c++;
  }

  /* Right */
  PROPAGATE_RIPPLE_P(c[P(-1,0)], cr ? *cr : *c, c[P(0,-1)], c[P(0,1)]);
  p++;
  c++;
  if (cr) cr += P(0,1);

 }


 /* Bottom Left */
 PROPAGATE_RIPPLE_P(cl ? *cl : *c, c[P(1,0)], c[P(0,-1)], cb ? *cb++ : *c);
 p++;
 c++;

 /* Bottom Row */
 for (x = 1; x < RIPPLE_D - 1; x++) {
  PROPAGATE_RIPPLE_P(c[P(-1,0)], c[P(1,0)], c[P(0,-1)], cb ? *cb++ : *c);
  p++;
  c++;
 }

 /* Bottom Right */
 PROPAGATE_RIPPLE_P(c[P(-1,0)], cr ? *cr : *c, c[P(0,-1)], cb ? *cb : *c);

}



void finalise_ripples(RIPPLE_INFO *ripple_info) {

 R_DISTURBANCE *rd, *next;

 ripple_info->frame++;

 if (ripple_info->frame >= RIPPLE_FRAMES) {
  ripple_info->frame = 0;
  ripple_info->swapped ^= 1;

  rd = ripple_info->dist;
  ripple_info->dist = 0;

  while (rd) {
   apply_disturbance(ripple_info->swapped ? rd->ripples->prev :
                                            rd->ripples->curr, rd);
   next = rd->next;
   destroy_disturbance(rd);
   rd = next;
  }
 }
}



void create_disturbance(RIPPLE_INFO *ripple_info,
                        RIPPLES *ripples,
                        float x, float y, float r, float h) {

 R_DISTURBANCE *rd;
 float t;

 if (ripples == 0) return;

 rd = ripple_info->dist;

 while (rd) {
  if (rd->ripples == ripples &&
      x - rd->x > -(1.0 / RIPPLE_D) &&
      x - rd->x <  (1.0 / RIPPLE_D) &&
      y - rd->y > -(1.0 / RIPPLE_D) &&
      y - rd->y <  (1.0 / RIPPLE_D) &&
      r - rd->r > -(1.0 / RIPPLE_D) &&
      r - rd->r <  (1.0 / RIPPLE_D)) {

   rd->h += h;

   t = h / rd->h;

   rd->x += t * (x - rd->x);
   rd->y += t * (y - rd->y);

   return;
  }
  rd = rd->next;
 }

 rd = malloc(sizeof(R_DISTURBANCE));

 if (!rd) return;

 rd->ripples = ripples;
 rd->x = x;
 rd->y = y;
 rd->r = r;
 rd->h = h;

 rd->next = ripple_info->dist;
 ripple_info->dist = rd;
}



void apply_disturbance(ripple_p *curr, R_DISTURBANCE *rd) {

 float k;
 float xl, yl;
 int xli, xhi, yli, yhi;
 int x, y;

 float zdy, dzdy;
 float z, dz, ddz;
 float zdx, dzdx;

 rd->x *= RIPPLE_D;
 rd->y *= RIPPLE_D;
 rd->r *= RIPPLE_D;

 /* The parabola is defined by:

    z = h - k*x*x - k*y*y

    where x and y are relative to the centre.

    It must be true that, at x*x + y*y = r*r, z = 0.
    0 = h - k*r*r
    h = k*r*r
    k = h/(r*r)
 */

 k = rd->h / (rd->r * rd->r);

 /* Set up us the x range. */

 xhi = (int)ceil(rd->x + rd->r);
 xhi = MIN(xhi, RIPPLE_D);

 xl = ceil(rd->x - rd->r);
 xl = MAX(xl, 0);
 xli = (int)xl;

 xl -= rd->x;

 /* Set up us the y range. */

 yhi = (int)ceil(rd->y + rd->r);
 yhi = MIN(yhi, RIPPLE_D);

 yl = ceil(rd->y - rd->r);
 yl = MAX(yl, 0);
 yli = (int)yl;

 yl -= rd->y;

 /* Set up us the pointer to the ripples. */
 curr += P(0,yli);

 /* Set up us zdy, and its rates of change with respect to y. */
 zdy = -k*yl*yl;
 dzdy = -k*(yl+1.0)*(yl+1.0) - zdy;

 /* Set up us zdx, and its rates of change with respect to x. */
 zdx = -k*xl*xl;
 dzdx = -k*(xl+1.0)*(xl+1.0) - zdx;

 /* ddz is the same with respect to each coordinate. */
 ddz = -2.0*k;


 for (y = yli; y < yhi; y++) {

  z = rd->h + zdy + zdx;
  dz = dzdx;

  for (x = xli; x < xhi; x++) {

   curr[x] += z; /* YAY! I FINALLY GET TO WRITE THIS LINE!!! :-) */

   z += dz;
   dz += ddz;
  }

  curr += P(0,1);

  zdy += dzdy;
  dzdy += ddz;
 }

}


/*
#define SAVE_RIPPLE_BMP
*/

#ifdef SAVE_RIPPLE_BMP
#include <stdio.h>
#endif

void create_ripple_colours(RIPPLE_COLOURS *rc,
                           int hr, int hg, int hb,
                           int sr, int sg, int sb,
                           int lr, int lg, int lb) {
 int x, y;
 unsigned char *c = &rc->c[0][0];

#ifdef SAVE_RIPPLE_BMP

 static int bmp_num = 0;
 BITMAP *bmp = create_bitmap_ex(24, 256, 256);
 if (!bmp) exit(37);

#endif

 /* Subtract the highlight from the sky colour. It may be used as an addition
    later.
 */
 hr -= sr;
 hg -= sg;
 hb -= sb;

 for (y = -128; y < 128; y++) {
  for (x = -128; x < 128; x++) {

   #define RIPPLE_Z 64

   float d = sqrt(x*x + y*y + RIPPLE_Z*RIPPLE_Z);
   float vx = x / d;
   float vy = y / d;
   float vz = RIPPLE_Z / d;

   float dot = 0.8164965809 * vz - 0.4082482905 * (vx + vy);

   #define H_THRESHOLD 0.99

   float t = 0.5 * (dot + 1.0);

   float r, g, b;

   if (vz >= 0.7071067812) {
    r = t * sr;
    g = t * sg;
    b = t * sb;
    if (dot >= H_THRESHOLD) {
     t = (dot - H_THRESHOLD) * (1.0 / (1.0 - H_THRESHOLD));
     r += t * hr;
     g += t * hg;
     b += t * hb;
    }
   } else {
    r = t * lr;
    g = t * lg;
    b = t * lb;

    ASSERT(dot < H_THRESHOLD);
   }

   ASSERT((int)r >= 0 && (int)r < 256);
   ASSERT((int)g >= 0 && (int)g < 256);
   ASSERT((int)b >= 0 && (int)b < 256);

   *c++ = makecol((int)r, (int)g, (int)b);

#ifdef SAVE_RIPPLE_BMP
   putpixel(bmp, 128 + x, 128 + y, makecol24((int)r, (int)g, (int)b));
#endif
  }
 }

#ifdef SAVE_RIPPLE_BMP
 {
  unsigned char ripple_bmp[256];
  sprintf(ripple_bmp, "data/files/truecol/ripple%c.bmp", '0' + bmp_num++);
  save_bitmap(ripple_bmp, bmp, 0);
 }
#endif
}



static int dr_scale;
static int dr_steps;
static BITMAP *dr_bmp;
static RIPPLE_COLOURS *dr_rc;
static int dr_y; /* Current row */

/* row*[] hold normals and deltas of current row. */
static int rowx[RIPPLE_D + 1];
static int rowy[RIPPLE_D + 1];
static int rowdx[RIPPLE_D + 1];
static int rowdy[RIPPLE_D + 1]; /* LOL */



static ripple_p dr_rl_rr[RIPPLE_D][RIPPLE_D];


static ripple_p *create_rl(ripple_p *r) {
 int yr;

 for (yr = 0; yr < RIPPLE_D; yr++) {
  dr_rl_rr[yr][0] = (r[0] << 1) - r[1];
  r += RIPPLE_D;
 }

 return &dr_rl_rr[0][0];
}


static ripple_p *create_rr(ripple_p *r) {
 int yr;

 for (yr = 0; yr < RIPPLE_D; yr++) {
  dr_rl_rr[yr][1] = (r[RIPPLE_D-1] << 1) - r[RIPPLE_D-2];
  r += RIPPLE_D;
 }

 return &dr_rl_rr[0][1];
}



static ripple_p *create_rt(ripple_p *r) {
 static ripple_p rt[RIPPLE_D];
 int xr;

 for (xr = 0; xr < RIPPLE_D; xr++) {
  rt[xr] = (r[P(0,0)] << 1) - r[P(0,1)];
  r++;
 }

 return rt;
}



static ripple_p *create_rb(ripple_p *r) {
 static ripple_p rb[RIPPLE_D];
 int xr;

 for (xr = 0; xr < RIPPLE_D; xr++) {
  rb[xr] = (r[P(0,RIPPLE_D-1)] << 1) - r[P(0,RIPPLE_D-2)];
  r++;
 }

 return rb;
}



static ripple_p *create_rtl(ripple_p *r) {
 static ripple_p rtl;
 rtl = (r[P(0,0)] << 1) - r[P(1,1)];
 return &rtl;
}



static ripple_p *create_rtr(ripple_p *r) {
 static ripple_p rtr;
 rtr = (r[P(RIPPLE_D-1,0)] << 1) - r[P(RIPPLE_D-2,1)];
 return &rtr;
}



static ripple_p *create_rbl(ripple_p *r) {
 static ripple_p rbl;
 rbl = (r[P(0,RIPPLE_D-1)] << 1) - r[P(1,RIPPLE_D-2)];
 return &rbl;
}



static ripple_p *create_rbr(ripple_p *r) {
 static ripple_p rbr;
 rbr = (r[P(RIPPLE_D-1,RIPPLE_D-1)] << 1) - r[P(RIPPLE_D-2,RIPPLE_D-2)];
 return &rbr;
}



static void draw_ripple_row() {

 int xr;

 int nx, ny, ndx, ndy; /* Used for interpolation within a row */
 int xi, yi; /* Interpolation counters */
 unsigned char *bp; /* Bitmap pixel pointer */

 for (xr = 0; xr < RIPPLE_D + 1; xr++) {
  rowdx[xr] = (rowdx[xr] - rowx[xr]) >> dr_scale;
  rowdy[xr] = (rowdy[xr] - rowy[xr]) >> dr_scale;
 }

 /* For each row of pixels... */
 for (yi = dr_steps; yi; yi--) {

  /* Render the row. */
  bp = dr_bmp->line[dr_y++];

  nx = rowx[0];
  ny = rowy[0];

  for (xr = 0; xr < RIPPLE_D; xr++) {
   ndx = (rowx[xr+1] - nx) >> dr_scale;
   ndy = (rowy[xr+1] - ny) >> dr_scale;

   for (xi = dr_steps; xi; xi--) {

    *bp++ = dr_rc->c[(ny >> 8) & 0xFF][(nx >> 8) & 0xFF];

    nx += ndx;
    ny += ndy;
   }
  }

  /* Interpolate the row*[] arrays. */
  for (xr = 0; xr < RIPPLE_D + 1; xr++) {
   rowx[xr] += rowdx[xr];
   rowy[xr] += rowdy[xr];
  }
 }
}



void draw_ripples(RIPPLE_COLOURS *rc, BITMAP *bmp, int scale,
                  RIPPLE_INFO *ripple_info,
                  RIPPLES *ripples,
                  RIPPLES *rip_l,
                  RIPPLES *rip_r,
                  RIPPLES *rip_t,
                  RIPPLES *rip_b,
                  RIPPLES *rip_tl,
                  RIPPLES *rip_tr,
                  RIPPLES *rip_bl,
                  RIPPLES *rip_br) {

 ripple_p *r, *rl, *rr, *rt, *rb, *rtl, *rtr, *rbl, *rbr;

 int xr, yr; /* Ripple point coordinates */

 if (ripples == 0) return;

 ASSERT(ripple_info);

 dr_scale = scale;
 dr_steps = 1 << scale;
 dr_bmp = bmp;
 dr_rc = rc;

 if (rip_l == 0 && rip_t == 0) rip_tl = 0;
 if (rip_r == 0 && rip_t == 0) rip_tr = 0;
 if (rip_l == 0 && rip_b == 0) rip_bl = 0;
 if (rip_r == 0 && rip_b == 0) rip_br = 0;

 if (ripple_info->swapped) {
  r = ripples->prev;
  rl = rip_l ? rip_l->prev : 0;
  rr = rip_r ? rip_r->prev : 0;
  rt = rip_t ? rip_t->prev : 0;
  rb = rip_b ? rip_b->prev : 0;
  rtl = rip_tl ? rip_tl->prev : 0;
  rtr = rip_tr ? rip_tr->prev : 0;
  rbl = rip_bl ? rip_bl->prev : 0;
  rbr = rip_br ? rip_br->prev : 0;
 } else {
  r = ripples->curr;
  rl = rip_l ? rip_l->curr : 0;
  rr = rip_r ? rip_r->curr : 0;
  rt = rip_t ? rip_t->curr : 0;
  rb = rip_b ? rip_b->curr : 0;
  rtl = rip_tl ? rip_tl->curr : 0;
  rtr = rip_tr ? rip_tr->curr : 0;
  rbl = rip_bl ? rip_bl->curr : 0;
  rbr = rip_br ? rip_br->curr : 0;
 }

 /* Points are arranged: ab */
 /*                      cd */
 #define COMPUTE_NORMAL(x, y, a, b, c, d) { \
  (x) = (128 << 8) + (a) - (d) - (b) + (c); \
  (y) = (128 << 8) + (a) - (d) + (b) - (c); \
 }


 if (rl)  rl  += P(RIPPLE_D-1, 0         ); else rl  = create_rl(r);
 if (rr)  rr  += P(0         , 0         ); else rr  = create_rr(r);
 if (rt)  rt  += P(0         , RIPPLE_D-1); else rt  = create_rt(r);
 if (rb)  rb  += P(0         , 0         ); else rb  = create_rb(r);
 if (rtl) rtl += P(RIPPLE_D-1, RIPPLE_D-1); else rtl = create_rtl(r);
 if (rtr) rtr += P(0         , RIPPLE_D-1); else rtr = create_rtr(r);
 if (rbl) rbl += P(RIPPLE_D-1, 0         ); else rbl = create_rbl(r);
 if (rbr) rbr += P(0         , 0         ); else rbr = create_rbr(r);


 /* Initialise top row of normals. */

 COMPUTE_NORMAL(rowx[0], rowy[0], rtl[0], rt[0],
                                  rl [0], r [0]);
 rt++;
 r++;

 for (xr = 1; xr < RIPPLE_D; xr++) {
  COMPUTE_NORMAL(rowx[xr], rowy[xr], rt[-1], rt[0],
                                     r [-1], r [0]);
  rt++;
  r++;
 }

 COMPUTE_NORMAL(rowx[RIPPLE_D], rowy[RIPPLE_D], rt[-1], rtr[0],
                                                r [-1], rr [0]);

 rl += RIPPLE_D;
 rr += RIPPLE_D;

 dr_y = 0;

 /* For each row in the ripple tile... */
 for (yr = 1; yr < RIPPLE_D; yr++) {

  /* Initialise next row of normals. */

  COMPUTE_NORMAL(rowdx[0], rowdy[0], rl[P(0,-1)], r[P(0,-1)],
                                     rl[P(0, 0)], r[P(0, 0)]);
  r++;

  for (xr = 1; xr < RIPPLE_D; xr++) {
   COMPUTE_NORMAL(rowdx[xr], rowdy[xr], r[P(-1,-1)], r[P(0,-1)],
                                        r[P(-1, 0)], r[P(0, 0)]);
   r++;
  }

  COMPUTE_NORMAL(rowdx[RIPPLE_D], rowdy[RIPPLE_D], r[P(-1,-1)], rr[P(0,-1)],
                                                   r[P(-1, 0)], rr[P(0, 0)]);

  rl += RIPPLE_D;
  rr += RIPPLE_D;

  draw_ripple_row();
 }

 /* Initialise final row of normals. */

 COMPUTE_NORMAL(rowdx[0], rowdy[0], rl [P(0,-1)], r [P(0,-1)],
                                    rbl[P(0, 0)], rb[P(0, 0)]);
 r++;
 rb++;

 for (xr = 1; xr < RIPPLE_D; xr++) {
  COMPUTE_NORMAL(rowdx[xr], rowdy[xr], r [P(-1,-1)], r [P(0,-1)],
                                       rb[P(-1, 0)], rb[P(0, 0)]);
  r++;
  rb++;
 }

 COMPUTE_NORMAL(rowdx[RIPPLE_D], rowdy[RIPPLE_D], r [P(-1,-1)], rr [P(0,-1)],
                                                  rb[P(-1, 0)], rbr[P(0, 0)]);

 draw_ripple_row();

}
