/*
 * This software is copyrighted as noted below.  It may be freely copied,
 * modified, and redistributed, provided that the copyright notice is
 * preserved on all copies.
 *
 * There is no warranty or other guarantee of fitness for this software,
 * it is provided solely "as is".  Bug reports or fixes may be sent
 * to the author, who may or may not act on them as he desires.
 *
 * You may not include this software in a program or other software product
 * without supplying the source, or without informing the end-user that the
 * source is available for no extra charge.
 *
 * If you modify this software, you should include a notice giving the
 * name of the person performing the modification, the date of modification,
 * and the reason for such modification.
 */
/*
 * raindrops.c - calculates raindrops hitting [more] water.  
 *
 * Author:      Raul Rivero
 *              Mihai P. Datcu
 *              Mathematics Dept.
 *              University of Oviedo
 * Date:        Fri 20 Aug 1993
 * Copyright (c) 1993, Mihai P. Datcu && Raul Rivero
 *
 */
#include <lug.h>
#include <lugfnts.h>

#ifdef WRITE_TIFF
#	define SCALE           .75
#else
#	define SCALE           .1
#endif
#define XSIZE			256
#define YSIZE			XSIZE
#define NOFRAMES                32
#define X0			512
#define Y0			512
#define SIGMA			(64.*4.)

float gasdev(), rand1();
float *smooth(), *apply_convolve();

main(argc, argv)
int argc;
char **argv;
{
  register i,j;
  float *buffer, *aux;
  int totalsize = XSIZE * YSIZE;
#ifdef WRITE_TIFF
  bitmap_hdr in;
#endif
  char name[80];

  if ( argc < 2 ) {
    fprintf( stderr, "Use: %s [frame_number]\n", argv[0]);
    exit( 1 );
  }

#ifdef WRITE_TIFF
  /*
   * We'll generate a TIFF image, so fill our
   * header.
   */
  in.magic = LUGUSED;
  in.xsize = XSIZE;
  in.ysize = YSIZE;
  in.depth = 8;
  in.colors = ( 1 << in.depth );
  in.cmap = create_bw_pallete();
  /* Get memory for raster information */
  in.r = (byte *) Malloc( totalsize );
#endif  /* WRITE_TIFF */

  /*
   * Get memory for main buffer. Although all operators will
   * be doubles, the buffer will be float 'cos Rayshade uses
   * that format.
   */
  buffer = (float *) Malloc( sizeof(float) * totalsize );

  /*
   * This could be a loop, but ...
   */
  i = atoi( argv[1] );

#ifdef ADD_NOISE
  /*
   * Reset the buffer.
   */
  fprintf( stderr, "Adding noise..." ), fflush( stderr );
  (void)rand1(-i);
  for ( j = 0, aux = buffer; j < totalsize; j++ )
      *aux++ = 0.005 * SCALE * gasdev(i);

  /*
   * Now, smooth it.
   */
  fprintf( stderr, "\nSmoothing noise..." ), fflush( stderr );
  aux = smooth( buffer, XSIZE, YSIZE );
  Free( buffer );
  buffer = aux;
#endif  /* ADD_NOISE */ 

  fprintf( stderr, "\nTime = %d", i);
  /*
   * Calculate all the hits.
   */
  /* 1st loop */
  calculate_wave( buffer, 100, 150, i );
  calculate_wave( buffer, 125,  75, i - 10 );
  calculate_wave( buffer, 225, 125, i - 20 );
  calculate_wave( buffer, 200,  50, i - 30 );
  calculate_wave( buffer,  25, 225, i - 40 );
  calculate_wave( buffer, 150, 175, i - 50 );
  calculate_wave( buffer,  50,  25, i - 60 );
  calculate_wave( buffer,  25, 125, i - 70 );
  /* 2nd loop */
  calculate_wave( buffer, 100, 150, i - 80 );
  calculate_wave( buffer, 125,  75, i - 90 );
  calculate_wave( buffer, 225, 125, i - 100 );
  calculate_wave( buffer, 200,  50, i - 110 );
  calculate_wave( buffer,  25, 225, i - 120 );
  calculate_wave( buffer, 150, 175, i - 130 );
  calculate_wave( buffer,  50,  25, i - 140 );
  calculate_wave( buffer,  25, 125, i - 150 );

  /* Smooth the resulting surface */
  fprintf( stderr, "\nSmoothing surface..." ), fflush( stderr );
  aux = smooth( buffer, XSIZE, YSIZE );
  Free( buffer );
  buffer = aux;

#ifdef WRITE_TIFF
  /*
   * Convert to an image ( from [0..1] to [0..255] ).
   */
  (void)tobyte( buffer, in.r, totalsize );
  /*
   * Write the TIFF file.
   */
  sprintf( name, "%03d.tif", i );
  fprintf( stderr, "\nWriting file %s \n", name );
  write_tiff_file( name, &in );
#else
  /*
   * Write the Rayshade's heightfield file.
   */
  sprintf( name, "%03d.hf", i );
  fprintf( stderr, "\nWriting file %s \n", name );
  write_heightfield_file( name, buffer, XSIZE );
#endif  /* WRITE_TIFF */
}

calculate_wave( buffer, xpos, ypos, time )
register float *buffer;
int xpos, ypos, time;
{
  if ( time >= 0 && time <= 15) {
    register int x, y;
    register double I = ((double) time) / ((double)NOFRAMES);
    register double X, Y;
    register double distance;
    register double PI_2 = 2. * M_PI;
    register double shift = sqrt( 2 * pow(time*8., 2.) );
    register double time_scale = SCALE * exp( -2.5 * ((double)time) / ((double)NOFRAMES) );
    register double scale;
    register double factor, yfactor;

    fprintf( stderr, "\n\tcalculating impact on %dx%d...", xpos, ypos ), fflush( stderr );
    for ( y = 0; y < YSIZE; y++ ) {
      Y = (y - ypos) * (y - ypos);
      yfactor = pow( y - ypos, 2. );
      for ( x = 0; x < XSIZE; x++ ) {
        X = (x - xpos) * (x - xpos);
        distance = sqrt( X + Y );
        factor = exp( -1 * pow( (sqrt(pow(x-xpos, 2.)+yfactor)-shift), 2. ) / SIGMA );
        scale  = factor * time_scale * exp( -7.5 * distance / ((double)XSIZE) );
        *buffer += scale * cos( PI_2 * ( I - distance / 16. ) );
        buffer++;
      }
    }
  }
}

float *
smooth(buffer, xsize, ysize)
float *buffer;
int xsize, ysize;
{
  double blur_matrix[9] = {   1./9.,   1./9.,    1./9.,
			      1./9.,   1./9.,    1./9.,
			      1./9.,   1./9.,    1./9.   };


  return apply_convolve( buffer, xsize, ysize, blur_matrix );
}

float *
apply_convolve(buffer, xsize, ysize, matrix)
float *buffer;
int xsize, ysize;
double matrix[9];
{
  register int i;
  int totalsize = xsize * ysize * sizeof(float);
  float *nbuffer;
  register float *fin;
  register float *base, *line_before;    /* pointers to new [blured] image */
  register float *last, *next_line;      /* pointers to old [zoomed] image */

  /* Allocate memory for new blured image */
  nbuffer= (float *) Malloc( totalsize );

  /*
   * First and last rows, and first and last pixels in
   * each row are not modified, so we copy image into new
   * buffer and keep all information.
   */
  bcopy( buffer, nbuffer, totalsize );

  /* Skip first and last row */
  ysize--;
  for ( i = 1; i < ysize; i++ ) {
    /* Skip first pixel of current line */
    base = nbuffer + i*xsize + 1;
    /* Point to current line */
    last = buffer + i*xsize + 1;
    /* Point to line before */
    line_before = last - xsize;
    /* Point to next line */
    next_line = last + xsize;
    /*
     * Pointer to last pixel for being modified of the current
     * line ( skip last pixel ).
     */
    fin = base + xsize - 1;
    while ( base < fin ) {
      /* Blur the current pixel */
      *base++ = matrix[4] * *last +
		matrix[3] * *(last-1)        + matrix[5] * *(last+1)        +
		matrix[1] * *line_before     + matrix[7] * *next_line       +
		matrix[0] * *(line_before-1) + matrix[2] * *(line_before+1) +
		matrix[6] * *(next_line-1)   + matrix[8] * *(next_line+1);
      /* Update pointers */
      next_line++;
      line_before++;
      last++;
    }
  }

  return nbuffer;
}

#ifdef ADD_NOISE 

/*
 * We need a Gasussian Random Numbers Generator.
 *
 * From Numerical Recipes in C.
 *
 */
float gasdev( idum )
int idum;
{
  static iset = 0;
  static float gset;
  float fac, r, v1, v2;

  if ( iset == 0 ) {
    do {
      v1 = 2. * rand1(idum) - 1.;
      v2 = 2. * rand1(idum) - 1.;
      r = v1*v1 + v2*v2;
    }while ( r >= 1. );
    fac = sqrt( -2. * log(r) / r );
    gset = v1 * fac;
    iset = 1;
    return v2*fac;
  }else {
    iset = 0;
    return gset;
  }
}

#define M1		259200
#define IA1		7141
#define IC1		54773
#define RM1		(1./M1)
#define M2		134456
#define IA2		8121
#define IC2		28411
#define RM2		(1./M2)
#define M3		243000
#define IA3		4561
#define	IC3		51349

float rand1( idum )
int idum;
{
  static long ix1, ix2, ix3;
  static float r[98];
  float temp;
  static iff = 0;
  int j;

  if ( idum < 0 || iff == 0 ) {
    iff = 1;
    ix1 = (IC1 - (idum))  % M1;
    ix1 = (IA1 * ix1 + IC1) % M1;
    ix2 = ix1 % M2;
    ix1 = (IA1 * ix1 + IC1) % M1;
    ix3 = ix1 % M3;
    for ( j = 1; j <= 97; j++ ) {
      ix1 = (IA1 * ix1 + IC1) % M1;
      ix2 = (IA2 * ix2 + IC2) % M2;
      r[j] = (ix1 + ix2 * RM2) * RM1;
    }
  }
  ix1 = (IA1 * ix1 + IC1) % M1;
  ix2 = (IA2 * ix2 + IC2) % M2;
  ix3 = (IA3 * ix3 + IC3) % M3;
  j = 1 + ( (97*ix3) / M3);
  if ( j > 97 || j < 1)
    perror( "RAN1: This cannot happen." );
  temp = r[j];
  r[j] = (ix1+ix2*RM2) * RM1;
  return temp;
}

#endif  /* ADD_NOISE */

/*
 * The output routines, if you choose the TIFF file format
 * you need the LUG library.
 *               |
 *               |
 *               --> ftp://ftp.uniovi.es/uniovi/mathdepth/src/liblug-1.0.tar.gz
 *
 * Else, the program'll generate Rayshade's heightfields.
 *
 */

#ifdef WRITE_TIFF

tobyte( in, out, count )
float *in;
byte *out;
int count;
{
  float aux;

  while ( count-- ) {
    aux = 255. * *in++;
    *out++ = CORRECT( aux );
  }
}

#else

write_heightfield_file( name, buffer, size )
char *name;
float *buffer;
int size;
{
  FILE *handle;

  /* Open the file descriptor */
  if ( name != NULL )
    handle = Fopen( name, "wb" );
  else handle = stdout;

  /* Write the bitmap */
  Fwrite( &size, sizeof(int), 1, handle );
  Fwrite( buffer, sizeof(float), size*size, handle );

  /* Close the file */
  Fclose( handle );

}

#endif  /* WRITE_TIFF */