/* 
Caves - Tribute to Herbert W. Franke

In this generative audiovisual work in the style of ASCII roguelike games, 
a lone adventurer  descends to explore fractal caves of increasing complexity. 
Sound and animation procedurally  generated by self contained C code. Released on 
the occasion of the  Tribute to Herbert W. Franke, commemorating his love of 
exploration in  diverse fields including speleology(cave exploration), mathematics, 
computer art and science fiction.

by 0xDEAFBEEF
Oct 2022

This self contained C file will compile and run to generate audio and image data.

Code conforms to the ISO c99 standard. Only library dependency is C math library. 
Pass '-lm' to gcc. Assumes architecture is little endian.

1. Create directory to store output images:
# mkdir 'frames'
2. Compile and run:
# gcc main.c -lm && ./a.out

This will produce:
a. Numbered image files in BMP file format stored in 'frames' directory. 24 FPS(frames
per second)
b. Audio file named 'output.wav' in WAV format, stereo, 44.1khz, 16bit depth

This is the 'raw' information representing digital audio signals
and image pixel intensities. 
This raw information can be encoded perceptually using whatever tools exist at the
time of reconstruction. At present, for example, opensource tools such as imagemagick and ffmpeg
can be used to encode images,video and audio in different format for popular distributions.

Examples:

1. Audio can be encoded to MP3.
# ffmpeg -i output.wav  -b 192k output.mp3

2. Images can be assembled into an animation and encoded to MP4, for example:
# ffmpeg -framerate 24 -i frames/frame%05d.bmp -crf 20 video.mp4

3. audio can be encoded with video:
# ffmpeg -i output.wav -i video.mp4 -c:v libx264 -c:a aac -b:a 192k -pix_fmt yuv420p  -shortest audiovisual.mp4

4. Screenshots of particular frames(in this example, frame 60), using imagemagick:
# convert frames/frame0060.bmp -scale 720x720 image.png

*/

char *seed="0x2c68ca657f166e1e7a5db2266ee1edca9daf2da5ecc689c39c1d8cfe0b9aad2e";

unsigned int token_param[8];
void init_token_params() {
    //insert value from smart contract function getTokenParam(tokenID,0)
  token_param[0] = 3338148;

  //params 1 - 7 unused for this series  
  token_param[1] = 0; 
  token_param[2] = 0; 
  token_param[3] = 0; 
  token_param[4] = 0; 
  token_param[5] = 0; 
  token_param[6] = 0;
  token_param[7] = 0;    
}

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>

int myrand();
#define rand() myrand()

#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288
#endif

#define SAMPLE_RATE 44100

long randp();
double drand ( double low, double high );
int ilimit(int low, int high, int x);

typedef struct {
  double freq;
  double phase;
  double amp;
  double pw; 
  int i; 

  // for DSF
  double freq2;
  double phase2;
  double w;
} gen_data;

typedef struct {
  int ftype;
  double fc;  
  double bw;  
  
  double pole;
  double a1;
  double a2;
  double rA; 
  double rB; 
  double rC; 
  double gain;
  double x[3];
  double y[3];

} iirfilter_params;
#define FTYPE_LP 1
#define FTYPE_HP 2
#define FTYPE_RES2 3
#define FTYPE_RES2A 4

double *make_bufn(int n);
double *make_buf(double dur,int*n);
void fill(double a, double *env,int n);
void delay(int b, double *data,int n);
void shift(int b, double *data,int n);
double maxsample(double *data,int n);
void normalize(double *data,int n);
void normalize_stereo(double *lchan, double *rchan,int n);
void mult(double *src, double *dst, int n);
void lerp(double a, double *src, double *dst, int n);
void add(double *src, double *dst, int n);
void scale(double a,double b, double *dst, int n);
void softclip(double gain,double offset,double *data,int n);
void gain(double g,double*data,int n);

double nond_drand ( double low, double high );
void dbgain(double db,double*data,int n);

//misc
double limit(double low,double high,double x);
double dmod(double x, double y);

void stereomixdown(double *lchan,double*rchan,int n);
void msmix(double *src,double dbgain, double pos, double *dstl, double *dstr,int n);

//bsets (breakpoint sets)
typedef struct bset {
  int n; 
  int max_n;
  int sr; 
  int interp;
  double param;
  double *t; 
  double *y; 
  double *iparam;
  int *itype; 
} bset;

#define INTERP_SAMPLEHOLD 1
#define INTERP_LINEAR 2
#define INTERP_CUBIC 3
#define INTERP_GAMMA 4
#define INTERP_EXP 5

void bset_init(bset *b);
void bset_destroy(bset *b);
int bset_search(bset *b, double t);
double bset_interp(bset *b, double t);
double bset_interp_linear(bset *b, double t);
double bset_interp_samplehold(bset *b, double t);
double bset_render(bset *b, double *data,int n);
void bset_dump(bset *b);
void bset_add(bset *b, double t, double y);
void bset_remove(bset *b, int j);
void bset_clear(bset *b);
double bset_length(bset *b);
double bset_last(bset *b);

//generators
void init_gen(gen_data *g);
double sinlu(double w);
double sin_gen(gen_data *g);
void whitenoise(double gain, double *data,int n);

//filters
void lowpass(double fc,double *data,int n);

void res2a_bset(bset *fenv,
		bset *bwenv,
		double fmult,
		double *data,int n);

//effects
void write_wav_header(FILE *f,int numSamples);
void write_wav_data(FILE *f,double *lchannel,double *rchannel,int num_samples);



void gain(double g,double*data,int n) {
  for (int i=0;i<n;i++) {
    data[i] *= g;
  }
}

void dbgain(double db,double*data,int n) {
  double g = exp(db/20 * log(10));
  
  for (int i=0;i<n;i++) {
    data[i] *= g;
  }
}

#define RAND_MAX_P 2147483647

//LCD PRNG with 8 streams
unsigned long long rs[8] = {1,2,3,4,5,6,7,8};
unsigned long long rs_initial[8] = {1,2,3,4,5,6,7,8};
int seed_index = 0;
int randp_seed(int i,unsigned long long entropy) {
  rs[i] = entropy;
  rs_initial[i] = entropy;  
  rs[i] = rs[i] % 2147483647;
  if (rs[i] <= 0) rs[i] += 2147483646;
}

int randp_getstream() { return seed_index; }
void randp_setstream(int i) {
  if (i>=0 && i< 8)  seed_index = i;
  else printf("randp_setstream index out of range");
}
long randps(int i) {  return rs[i] = rs[i] * 16807 % 2147483647;} //from specific stream
long randp() {  return rs[seed_index] = rs[seed_index] * 16807 % 2147483647;}

//resets stream to beginning
void randp_reset(int i) {
  rs[i] = rs_initial[i];
  rs[i] = rs[i] % 2147483647;
  if (rs[i] <= 0) rs[i] += 2147483646;  
}

//advances stream to specific point
void randp_advance(int i, int num) {
  for (int j=0;j<num;j++)
    randps(i);
}

//sample from uniform distribution
double drand ( double low, double high ) {
  return ( (double)randp() * ( high - low ) ) / (double)RAND_MAX_P + low;
}
//random integer
int irand ( int low, int high ) {
  return ( randp() % ( high - low ) )  + low;
}

//sample normal distribution via central limit theorem
double nrand() {
  double x = 0;
  for (int i=0;i<16;i++) {
    x += (double)randp()/RAND_MAX_P;
  }
  x -= 16/2.0;
  x /= sqrt(16/12.0);
  
  return x;
}

//truncated normal distribution
double tnrand(double mean, double sd, double low, double high) {
  int k = 0;
  while (1) {
    double x = nrand() * sd + mean;
    if (x > low && x < high) return x;
    if (k++ > 1000) {
      printf("truncated normal distribution taking more than 1000 rolls\n");
      return mean;
    }
  }
}

double choose2(double c1,double c2) { return (randp() % 2 ==0 ? c1 : c2); }
int ichoose2(int c1,int c2) { return (randp() % 2 ==0 ? c1 : c2); }
double choose3(double c1,double c2,double c3) {
  int c=  randp() % 3;
  if (c==0) return c1;
  if (c==1) return c2;
  if (c==2) return c3;

}

// adds src track with dst track
void add(double *src, double *dst, int n) {
  for (int i=0;i<n;i++) { dst[i] += src[i]; }
}

void invert(double *src,int n) {
  for (int i=0;i<n;i++) { src[i] = 1.0 - src[i]; }
}

//modulo operator for floating point
double dmod(double x, double y) {
  return x - (int)(x/y) * y;
}

// multiplies src track with dst track (ring modulator)
void mult(double *src, double *dst, int n) {
  for (int i=0;i<n;i++) { dst[i] *= src[i]; }
}

// scale and offset
void scale(double a,double b, double *dst, int n) {
  for (int i=0;i<n;i++) { dst[i] = dst[i]*a + b; }
}

void fill(double a, double *env,int n) {
  for (int i=0;i<n;i++) {
    env[i] = a;
  }
}

// shift back in time by 'b', and zero fill the rest
// to account for group delay in filters
void shift(int b, double *data,int n) {
  for (int i=0;i<n;i++) data[i] = data[i+b];
  for (int i=n-b;i<n;i++) data[i] = 0;
}
//delay
void delay(int b, double *data,int n) {
  for (int i=n-1;i>=b;i--) data[i] = data[i-b];
  for (int i=0;i<b;i++) { data[i] = 0; }
}


double *make_buf(double dur,int*n) {
  int num_samples=dur * SAMPLE_RATE;
  if (n!=NULL) { *n = num_samples; }
  return (double *)calloc(num_samples, sizeof(double));
}

double *make_bufn(int n) {
  return (double *)calloc(n, sizeof(double));
}
int *make_ibufn(int n) {
  return (int *)calloc(n, sizeof(int));
}

double *dup_buf(double *src,int n) {
  double *dst= (double *)malloc(n * sizeof(double));
  for (int i=0;i<n;i++)
    dst[i] = src[i];
  return dst;
}
// crossfade src track with dst track
void lerp(double a, double *src, double *dst, int n) {
  for (int i=0;i<n;i++) { dst[i] = a*src[i] + (1.0-a)*dst[i]; }
}
// normalize, optionally pass maximum to use
void normalize(double *data,int n) {
  double m = maxsample(data,n);
  if (m < 1e-009) { printf("Normalize: maxsample < 1e-009\n"); return; }
  gain(1.0/m,data,n);
}
double maxsample(double *data,int n) {
  double m = 0;
  for (int i=0;i<n;i++) {
    double x = fabs(data[i]);
    if (x > m) m = x;
  }
  return m;
}

double mean(double *data, int n) {
  if (n==0) return 0;
  double x = 0;
  for (int i=0;i<n;i++) {
    x += data[i];
  }
  return x/n;
}

void removedc(double *data,int n) {
  double m = mean(data,n);
  for (int i=0;i<n;i++) data[i] -=m;
}

void softclip(double gain,double offset,double *data,int n) {
  for (int i=0;i<n;i++) {
    double x = data[i] * gain + offset;
    if (x < -1.0) {
      data[i] = -(2.0/3.0);
    } else if (x > 1.0) {
      data[i] = 2.0/3.0;
    } else {
      data[i] = x - (x*x*x)/3.0;
    }
    data[i] /= gain;
    data[i] *= 3.0/2.0;
  }
}

void stereomixdown(double *lchan,double*rchan,int n) {
  double m = maxsample(lchan,n);
  double m2 = maxsample(rchan,n);
  if (m < m2) m = m2;
  gain(.98/m,lchan,n);
  gain(.98/m,rchan,n);
  
  FILE *f = fopen("output.wav","wb");
  write_wav_header(f,n);
  write_wav_data(f,lchan,rchan,n);  
  fclose(f);
}

void toBigEndian(int *a) {
  int b=0x0000;
  for (int i=0;i<4;i++) {
    b <<= 8;    
    b |= (*a & 0xff);
    *a >>= 8;
  }
  *a = b;
}

typedef struct typeWavHeader{
  //Riff header
  int ChunkID;
  int ChunkSize;
  int Format;

  // WAVE format
  int Subchunk1ID;
  int Subchunk1Size;
  short AudioFormat;
  short NumChannels;
  int SampleRate;
  int ByteRate;
  short BlockAlign;
  short BitsPerSample;

  //Data
  int Subchunk2ID;
  int Subchunk2Size;
} WavHeader;

void write_wav_header(FILE *f,int numSamples) {
  WavHeader h;

  h.ChunkID = 0x52494646;
  h.ChunkSize = 0;
  h.Format = 0x57415645;
  h.Subchunk1ID = 0x666d7420;
  h.Subchunk1Size = 16;
  h.AudioFormat = 1;
  h.NumChannels = 2; //Stereo=2,Mono=1
  h.SampleRate = 44100;
  h.BitsPerSample = 16; //16 bit
  h.ByteRate = h.SampleRate * h.NumChannels * h.BitsPerSample/8;
  h.BlockAlign = h.NumChannels * h.BitsPerSample/8;

  h.Subchunk2ID = 0x64617461;
  h.Subchunk2Size = numSamples * h.NumChannels * h.BitsPerSample/8;
  
  h.ChunkSize = 36+h.Subchunk2Size; //now calculate this

  //switch big endian
  toBigEndian(&h.ChunkID);
  toBigEndian(&h.Format);
  toBigEndian(&h.Subchunk1ID);
  toBigEndian(&h.Subchunk2ID);
  
  fwrite(&h,sizeof(h),1,f);
}

void write_wav_data(FILE *f,double *lchannel,double *rchannel,int num_samples) {
  short d;
  for (int i =0;i<num_samples;i++) {
    d = lchannel[i] * 32767;
    fwrite(&d,sizeof(short),1,f);    
    d = rchannel[i] * 32767;
    fwrite(&d,sizeof(short),1,f);
  }
}

//digital resonator described by Platt in speech synthesis
void res2a(double fc,double bw, double *data,int n) {
  double C = -1.0 * exp(-2.0*M_PI*bw/SAMPLE_RATE);
  double B = 2.0 * exp(-1.0*M_PI*bw/SAMPLE_RATE) * cos(2*M_PI*fc/SAMPLE_RATE);
  double A = 1.0 - B - C;
  
  double x[3] = {0,0,0}; double y[3] = {0,0,0};
  
  for (int i=0;i<n;i++) {
    x[0] = x[1];  x[1] = x[2]; x[2] = data[i];
    y[0] = y[1]; y[1] = y[2];
    y[2] = A * x[2] + B * y[1] + C *y[0];
    data[i] = y[2];
  }
}

//1st order IIR lowpass at cutoff freq fc  (-3db point)
void lowpass(double fc,double *data,int n) {
  double pole = 1 - 4*fc/SAMPLE_RATE; //location of pole in transfer function H(z)
  double GAIN = 2.0/(1-pole);
  
  double x[2] = {0,0}; double y[2] = {0,0};
  for (int i=0;i<n;i++) {
    x[0] = x[1];
    x[1] = data[i]/GAIN;
    y[0] = y[1];
    y[1] = (x[0] + x[1]) + (pole * y[0]);
    data[i] = y[1];
  }
}

//1st order IIR highpass at cutoff freq fc  (-3db point)
void highpass(double fc,double *data,int n) {
  double pole = 1 - 4*fc/SAMPLE_RATE; //location of pole in transfer function H(z)
  double GAIN = 2.0/(1+pole);
  
  double x[2] = {0,0}; double y[2] = {0,0};
  for (int i=0;i<n;i++) {
    x[0] = x[1];
    x[1] = data[i]/GAIN;
    y[0] = y[1];
    y[1] = (x[0] - x[1]) + (pole * y[0]);
    data[i] = y[1];
  }
}

//exponential envelope
void expenv(double tc, double *env, int env_n) {
  double lambda = 1.0/tc;
  for (int i=0;i<env_n;i++) {
    double t = (double)i/SAMPLE_RATE;
    env[i] = exp(-lambda*t);
  }
}
//render breakpoint set to an array of doubles
double bset_rendern(bset *b, double *data,int n) {
  double t = 0;
  int k=0;
  while (t < b->t[b->n-1] && k<n) {
    data[k++] = bset_interp(b,t);
    t+=1.0/b->sr;
  }
}


void apply_amp(bset *b, double *src, int n) {
  double *env = make_bufn(n*2);
  bset_rendern(b,env,n);
  mult(env,src,n);
  free(env);
}

//exponential attack decay envelope
void expad(double a,double d, double *env, int env_n) {
  int i=0;
  double t =0;
  double lambda = 1.0/(a/6);
  if (a < 1e-9) lambda = 1e9;
  
  int ac = 0;

  while (i<a*SAMPLE_RATE && (ac+i) < env_n) {
    t = (double)i/SAMPLE_RATE;
    env[ac+i] = exp(-lambda*t) * (0-1.0) + 1.0;
    i++;
  }
  
  ac += i; i=0;
  lambda = 1.0/ (d/6);
  if (d < 1e-9) lambda = 1e9;
  
  while (i<(d)*SAMPLE_RATE && (ac+i) < env_n) {
    t = (double)i/SAMPLE_RATE;    
    env[ac+i] = exp(-lambda*t);
    i++;
  }
  
  //zero fill rest of envelope
  while (ac+i < env_n) {
    env[ac+i++] = 0;
  }
}


void bandpass(double fl,double fh,double *data,int n) {
  highpass(fl,data,n);
  lowpass(fh,data,n);
}


void normalize_stereo(double *lchan, double *rchan,int n) {
  double m = maxsample(lchan,n);
  double m2 = maxsample(rchan,n);
  if (m < m2) m = m2;
  if (m < 1e-009) { printf("Normalize: maxsample < 1e-009\n"); return; }  
  gain(1.0/m,lchan,n);
  gain(1.0/m,rchan,n);
}

void whitenoise(double gain, double *data,int n) {
  for (int i=0;i<n;i++) {
    data[i] = drand(-gain,gain);
  }
}

void msmix(double *src,double dbgain, double pos, double *dstl, double *dstr,int n) {
  pos = limit(-1.0,1.0,pos);
  double g = exp(dbgain/20 * log(10));  
  for (int i=0;i<n;i++) {
    //constant power (-3db) pan law
    dstl[i] += g*cos((pos+1.0)*0.25*M_PI)*src[i];
    dstr[i] += g*sin((pos+1.0)*0.25*M_PI)*src[i];
  }
}
void mmmix(double *src,double dbgain, double *dst,int n) {
  double g = exp(dbgain/20 * log(10));  
  for (int i=0;i<n;i++) { dst[i] += g*src[i]; }
}

void tomono(double *lchan, double *rchan, double *mono, int n) {
  lerp(1.0,lchan,mono,n);
  add(rchan,mono,n);
}

double limit(double low,double high,double x) {
  return fmin(fmax(x,low),high);
}
int ilimit(int low, int high, int x) {
  return (x > high ? high : (x < low ? low : x));
}

// BSETS

void bset_init(bset *b) {
  b->interp = INTERP_LINEAR; //linear by default
  b->n = 0;  
  b->max_n = 2048*32*4;
  b->sr = SAMPLE_RATE;
  b->t = (double *)calloc(b->max_n,sizeof(double));
  b->y = (double *)calloc(b->max_n,sizeof(double));
  b->iparam = (double *)calloc(b->max_n,sizeof(double));     //param (gamma, tau);
  b->itype = (int *)calloc(b->max_n,sizeof(int));     //interpolation type
}

double bset_length(bset *b) {if (b->n==0) return 0;return b->t[b->n-1];
}

double bset_last(bset *b) {if (b->n==0) return 0;return b->y[b->n-1];}

void bset_destroy(bset *b) {
  if (b->t) free(b->t);
  if (b->y) free(b->y);
  if (b->iparam) free(b->iparam);
  if (b->itype) free(b->itype);    
  b->t =0;
  b->y =0;
}


//returns index of left bound of interval containing t
int bset_search(bset *b, double t) {
  if (b->n==0) return 0;
  int l = 0; int h = b->n;
  int i=0;
  do {
    i = (l+h)/2;
    //    printf("l: %d, h: %d, i: %d\n",l,h,i);    
    if (i==0) return t < b->t[0] ? 0 : 1;
    if (i >= b->n) return b->n;
    if (t >=  b->t[i-1] && t < b->t[i]) return i;
    if (l==h) return b->n;
    if (t >= b->t[i]) l=i+1; else h=i;
  } while (1);
  
  return i;
}

//returns index of left bound of interval containing t
int bset_search_bak(bset *b, double t) {
  int i=0;

  for (i=0;i<b->n;i++) {
    if (t < b->t[i]) break;
  }

  return i;
}


//interpolate breakpoint set at time t, linear interpolation
double bset_interp_linear(bset *b, double t) {
  int i = bset_search(b,t) - 1;
  if (i >= 0 && i < b->n-1) {
      //linear interpolation
      double a = (t - b->t[i])/ (b->t[i+1]-b->t[i]);
      return (1-a)*b->y[i] + a*b->y[i+1];
  }
  return 0;
}


//interpolate breakpoint set at time t, sample/hold
double bset_interp_samplehold(bset *b, double t) {
  int i = bset_search(b,t) - 1;
  if (i >= 0 && i < b->n-1) {
      //sample and hold      
      return b->y[i];    
  }
  return 0;
}

//interpolate breakpoint set, using the indicated interpolation method for that segment
double bset_interp(bset *b, double t) {
  int i = bset_search(b,t) - 1;
  if (i >= 0 && i < b->n-1) {
    switch(b->itype[i]) {
    case INTERP_SAMPLEHOLD:
      return bset_interp_samplehold(b,t);
      break;
    case INTERP_LINEAR:
      return bset_interp_linear(b,t);
      break;      
    default:
      return bset_interp_samplehold(b,t); // sample hold by default
      break;      
    }
  }
  return 0;
}
  
//render breakpoint set to an array of doubles
double bset_render(bset *b, double *data,int n) {
  double t = 0;
  int k=0;
  while (t < b->t[b->n-1] && k<n) {
    data[k++] = bset_interp(b,t);
    t+=1.0/b->sr;
  }
}



void bset_add(bset *b, double t, double y) {
  if (b->n >= b->max_n) { printf("more than max_n %d entries breakpoint set",b->max_n); return; }

  //insertion sort  
  int i = bset_search(b,t);
  
  for (int j=b->n;j>i;j--) {
    b->t[j] = b->t[j-1]; b->y[j] = b->y[j-1];
    b->itype[j] = b->itype[j-1];
    b->iparam[j] = b->iparam[j-1];    
  }
  b->t[i] = t; b->y[i] = y;
  b->itype[i] = b->interp; //default interpolation method for segment
  b->iparam[i] = b->param; //default interpolation param
  b->n++;
}

void bset_remove(bset *b, int j) {
  if (j >= b->n) return;

  for (int i = j;i<b->n-1;i++) {
    b->t[i] = b->t[i+1]; b->y[i] = b->y[i+1];
    b->itype[i] = b->itype[i+1];
    b->iparam[i] = b->iparam[i+1];    
  }
  b->n--;
}

//overwrite with linear segment
void bset_segment(bset *b,double t1, double t2, double y1, double y2) {
  //first remove any in the way
  for (int i=b->n-1;i>=0;i--) {
    if (b->t[i] > t1 && b->t[i] <= t2) {
      bset_remove(b,i);
    }
  }
  bset_add(b,t1,y1);
  bset_add(b,t2,y2);  
}

void bset_copy(bset *b, bset *dst) {
  bset_init(dst);
  for (int i=0;i<b->n;i++) {
    bset_add(dst,b->t[i],b->y[i]);
  }
}

void bset_clear(bset *b) { for (int i=0;i<b->n;i++) { b->t[i] = 0; b->y[i] = 0; b->itype[i] = 0; b->iparam[i] = 0;}  b->n =0; }
void bset_sort(bset *b) {
  bset b2;
  bset_init(&b2);
  
  for (int i=0;i<b->n;i++) {
    bset_add(&b2,b->t[i],b->y[i]);
  }
  bset_clear(b);
  for (int i=0;i<b2.n;i++) {
    bset_add(b,b2.t[i],b2.y[i]);
  }
}

double linear_interp(double p, double *data,int n) {
  if (p < 0) return 0;

  double a1 = p - floor(p);
  int i1 = (int)p;
  int i2 = i1 + 1;
  if (i2 >= n) return 0;

  return a1*data[i2] + (1.0-a1)*data[i1];
}


//update filter coefficients based on change in center frequency, bandwidth etc
void iir_update_params(iirfilter_params *f) {
  if (f->ftype==FTYPE_LP) {
    f->pole = 1 - 4*f->fc/SAMPLE_RATE; 
    f->gain = 2.0/(1-f->pole);
  } else if (f->ftype == FTYPE_HP) {
    f->pole = 1 - 4*f->fc/SAMPLE_RATE; 
    f->gain = 2.0/(1+f->pole);    
  } else if (f->ftype == FTYPE_RES2) {
    double R = exp(-M_PI*f->bw/SAMPLE_RATE); //pole radius in terms of bandwidth
    double A = 2*M_PI*f->fc/SAMPLE_RATE; // pole angle

    f->a1 = -2*R*cos(A);
    f->a2 = R*R;
    f->gain = 100 * (1-R*R)/2.0;
  } else if (f->ftype == FTYPE_RES2A) {
    f->rC = -1.0 * exp(-2.0*M_PI*f->bw/SAMPLE_RATE);
    f->rB = 2.0 * exp(-1.0*M_PI*f->bw/SAMPLE_RATE) * cos(2*M_PI*f->fc/SAMPLE_RATE);
    f->rA = 1.0 - f->rB - f->rC;
  }
}

void init_iirfilter(iirfilter_params *f) {
  for (int i=0;i<3;i++) { f->x[i] = 0; f->y[i] = 0; }  
  iir_update_params(f);
}

//one step of iir filter
double iir_apply(iirfilter_params *f, double s) {
  if (f->ftype==FTYPE_LP) {
    f->x[0] = f->x[1];
    f->x[1] = s/f->gain;
  
    f->y[0] = f->y[1];
    f->y[1] = (f->x[0] + f->x[1]) + (f->pole * f->y[0]);
    return f->y[1];
  } else if (f->ftype==FTYPE_HP) {
    f->x[0] = f->x[1];
    f->x[1] = s/f->gain;
  
    f->y[0] = f->y[1];
    f->y[1] = (f->x[0] - f->x[1]) + (f->pole * f->y[0]);
    return f->y[1];
  } else if (f->ftype==FTYPE_RES2) {
    double *x = f->x; double *y= f->y;
    x[0] = x[1];
    x[1] = x[2];
    x[2] = s * f->gain;

    y[0] = y[1];
    y[1] = y[2];
    y[2] = x[2] - x[0] - f->a1* y[1] -f->a2* y[0];
    return y[2];
  } else if (f->ftype==FTYPE_RES2A) {
    double *x = f->x; double *y= f->y;
    x[0] = x[1];  x[1] = x[2]; x[2] = s;
    y[0] = y[1]; y[1] = y[2];
    y[2] = f->rA * x[2] + f->rB * y[1] + f->rC *y[0];
    return y[2];
  } else {
    //unknown filter, just pass to output
    return s;
  }
  
  return s;
}

double swap(double *a, double *b) {
  double tmp = *a;
  *a = *b;
  *b = tmp;
}

//simple delay using delay buffer, and optional hp and lp filters
void ffdelay(double offset, double dur, double feedback, double fl, double fh, double *src, double*dst, int n) {
  feedback = exp(feedback/20 * log(10));
  
  int d = (dur * SAMPLE_RATE) + 1;
  
  double *buf = (double *)malloc(d * sizeof(double));
  for (int i=0;i<d;i++) buf[i] = 0;

  iirfilter_params filt_hp;
  filt_hp.ftype = FTYPE_HP;
  filt_hp.fc = fl;
    
  iirfilter_params filt_lp;
  filt_lp.ftype = FTYPE_LP;
  filt_lp.fc = fh;
  
  init_iirfilter(&filt_hp);
  init_iirfilter(&filt_lp);

  int p = 0;

  for (int i=0;i<n;i++) {
    //double x = buf[p];
    
    //apply any iir filters to buf[p] here
    if (fl > .001) buf[p] = iir_apply(&filt_hp, buf[p]);
    if (fh > .001) buf[p] = iir_apply(&filt_lp, buf[p]);      
    buf[p] = src[i] + feedback*buf[p];
    //data[i] += gain * x;
    dst[i] += buf[p];
    
    if (++p >= d) p = 0;
  }

  delay(offset*SAMPLE_RATE,dst,n); //offset delay
  
  free(buf);
}

//generators

void init_gen(gen_data *g) {
  g->freq = 440;

  g->phase = 0;

  g->pw = 0.50;
  g->amp = 1.0;
  g->i = 0;

  //for DSF
  g->freq2 = 440;  
  g->phase2 = 0;
  g->w = 0.5;
}

void sin_tone(int freq,double start,double dur,double amp,double *data) {
  int s;
  for (int i=0;i<dur*SAMPLE_RATE;i++) {
    data[i + (int)(start*SAMPLE_RATE)]  = amp * sinlu(freq * 2 * M_PI * (double)i/SAMPLE_RATE);
  }
}

double sin_gen(gen_data *g) {
  g->phase += 2*M_PI*g->freq/SAMPLE_RATE; 
  if (g->phase > 2*M_PI) g->phase -= 2*M_PI;
  g->phase = dmod(g->phase,2*M_PI);
  double x  = g->amp * sinlu(g->phase);
  return x;
}    

// sin lookup table
//int sinl_table_n=2048*256; //4mb lookup table gives about 70db SNR to math sin()
int sinl_table_n=2048*4; //4mb lookup table gives about 70db SNR to math sin()
double *sinl_table=0;
double sinlu(double w) {
  if (sinl_table==0) {
    printf("Building %d sin lookup table.. ",sinl_table_n);
    sinl_table = (double*)malloc(sinl_table_n * sizeof(double));
    for (int i=0;i<sinl_table_n;i++) {
      sinl_table[i] = sin((double)i*2.0*M_PI/(double)sinl_table_n);
    }
    printf("done.\n");
  }

  w = dmod(w, 2.0*M_PI);
  if (w < 0) w += 2.0*M_PI;
  
  double p = ((double)w*sinl_table_n/(2.0*M_PI));
  //return sinl_table[(int)round(p)];
  //optional linear interpolation of sin table

  double a1 = p - floor(p);
  int i1 = (int)p;
  int i2 = i1 + 1;
  if (i2 >= sinl_table_n) i2 = 0;

  return a1*sinl_table[i2] + (1.0-a1)*sinl_table[i1];
}

void nfadein(int d, double *src,int n) {
  double dur = (double)d/SAMPLE_RATE;  
  if (d > n) d = n;
  double t = 0;
  for (int i=0;i<d;i++) {
    double a = t/dur;
    src[i] *= a;
    t += 1.0/SAMPLE_RATE;    
  }
}
void fadein(double dur, double *src,int n) {
  int d = dur*SAMPLE_RATE;
  nfadein(d,src,n);
}

void nfadeout(int d, double *src,int n) {    
  double dur = (double)d/SAMPLE_RATE;
  if (d > n) d = n;
  double t = 0;
  for (int i=n-d;i<n;i++) {
    double a = (1-t/dur);
    src[i] *= a;
    t += 1.0/SAMPLE_RATE;    
  }
}
void fadeout(double dur, double *src,int n) {
  int d = dur*SAMPLE_RATE;
  if (d > n) d = n;
  double t = 0;
  for (int i=n-d;i<n;i++) {
    double a = (1-t/dur);
    src[i] *= a;
    t += 1.0/SAMPLE_RATE;    
  }
}

//visual code

#define CHAR_WIDTH 8 
#define CHAR_HEIGHT 12 

typedef unsigned int int32;
typedef short int16;
typedef unsigned char byte;

typedef struct img_t {
  int h;
  int w;
  int n; //number of values. 1 for grayscale, 3 for rgb
  double *p; //pixel data
} img_t;


typedef struct {
  int use_lines; //for '/-|\' type lines
  double c; //otherwise the character to use;
} drawlinesettings_t;

int ind(int x, int y,int w);

void val_to_img_ds(int w,int h, int os, double *val,img_t *im);
void val_to_img(int w,int h, double *val,img_t *im);

void write_frame(int w,int h,double os,double *val,img_t *img,int *fn,double *t);
void writebmp(img_t *i,char *filename);
void init_img(img_t *i,int n, int w, int h);
double getpixel(img_t *i, int x, int y, int n);
double *getpixelp(img_t *i, int x, int y);


//end header

//index into 2d array
int ind(int x, int y,int w) {
  return y*w + x;
}

//reverse index
void dni(int i, int *x, int *y,int w) {
  *y = i/w;
  *x = i - *y*w;
}

double getpixel(img_t *i, int x, int y, int n) {
  int index = x*i->h*i->n + y*i->n;
  double *p = &(i->p[index]);
  return p[n];
}

double *getpixelp(img_t *i, int x, int y) {
  int index = x*i->h*i->n + y*i->n;
  double *p = &(i->p[index]);
  return p;
}

void val_to_img(int w,int h, double *val,img_t *im) {
  for (int x=0;x<w;x++) { for (int y=0;y<h;y++) {
      double *v = &val[ind(x,y,w)];
      int index = x*im->h*im->n + y*im->n;
      double *p = &(im->p[index]);
      double f = *v;
      if (f < 0) f = 0; if (f > 1.0) f = 1.0;
      p[0] = f; p[1] = f; p[2] = f;

    }
  }
}

//val to img with downsampling
void val_to_img_ds(int w,int h, int os, double *val,img_t *im) {
  double os_v = 1.0/(double)(os*os);
  double *val_ds = make_bufn(((double)(w*h)/(os*os)));
  int w2 = w/os;
  int h2 = h/os;  

  for (int x=0;x<w;x++) {
    for (int y=0;y<h;y++) {
      int x2 = x/os;
      int y2 = y/os;
      val_ds[ind(x2,y2,w2)] += os_v * val[ind(x,y,w)];
    }
  }
  //  normalize(val_ds,w2*h2);  
  val_to_img(w2,h2,val_ds,im);
  free(val_ds);
}

void writebmp(img_t *i,char *filename) {
  int WIDTH = i->w;
  int HEIGHT= i->h;
  
  unsigned int headers[13];
  FILE * outfile;
  int extrabytes;
  int paddedsize;
  int x; int y; int n;
  int red, green, blue;

  extrabytes = 4 - ((WIDTH * 3) % 4);
  if (extrabytes == 4)
    extrabytes = 0;

  paddedsize = ((WIDTH * 3) + extrabytes) * HEIGHT;

  headers[0]  = paddedsize + 54;      
  headers[1]  = 0;                    
  headers[2]  = 54;                   
  headers[3]  = 40;                   
  headers[4]  = WIDTH; 
  headers[5]  = HEIGHT;
  headers[7]  = 0;                    
  headers[8]  = paddedsize;           
  headers[9]  = 0;                    
  headers[10] = 0;                    
  headers[11] = 0;                    
  headers[12] = 0;                    

  outfile = fopen(filename, "wb");
  fprintf(outfile, "BM");

  for (n = 0; n <= 5; n++) {
    fprintf(outfile, "%c", headers[n] & 0x000000FF);
    fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
    fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
    fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
  }
  fprintf(outfile, "%c", 1);
  fprintf(outfile, "%c", 0);
  fprintf(outfile, "%c", 24);
  fprintf(outfile, "%c", 0);

  for (n = 7; n <= 12; n++) {
    fprintf(outfile, "%c", headers[n] & 0x000000FF);
    fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
    fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
    fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
  }

  int stride = (i->h*i->n);
  unsigned char *data = (unsigned char*)malloc(WIDTH*3);

  for (y = HEIGHT - 1; y >= 0; y--) {
    double *p = getpixelp(i,0,y);

    for (x = 0; x <= WIDTH - 1; x++)  {
      data[x*3] = p[2]*255; 
      data[x*3+1] = p[1]*255;
      data[x*3+2] = p[0]*255;
      p+= stride;
    }
    fwrite(data,1,WIDTH*3,outfile);      
    if (extrabytes) {
      for (n = 1; n <= extrabytes; n++) {
	fprintf(outfile, "%c", 0);
      }
    }
  }

  fclose(outfile);
  return;
}

void write_frame(int w,int h,double os,double *val,img_t *img,int *fn,double *t) {
  val_to_img_ds(w,h,os,val,img);    
  char filename[128];
  sprintf(filename,"frames/frame%05d.bmp",*fn+1);
  writebmp(img,filename);
  
  if (*fn % 10==0) {
    printf("Frame %d, seed is %llu\n",*fn+1,rs[1]);
  }
  
  *t += 1/24.0;
  *fn = *fn + 1;
}

void write_pic(int w,int h,double os,double *val,img_t *img,int fn) {
  val_to_img_ds(w,h,os,val,img);    
  char f[128];
  sprintf(f,"frames/pic%d.bmp",fn);
  writebmp(img,f);
}

drawlinesettings_t dls;
void init_drawline() {
  dls.use_lines = 1;
  dls.c = '.';
}
int line_mode;
void val_line(int x0,int y0,
	      int x1,int y1,
	      double v,
	      double *val,int w, int h) {

  int dx = abs(x1-x0);
  int sx = x0 < x1 ? 1 : -1;
  int dy = -1.0*abs(y1-y0);
  int sy = y0<y1 ? 1 : -1;
  int err = dx+dy;
  int x = x0;
  int y = y0;
  int k=0;

  while (1) {
    if (x >= 0 && x < w && y >= 0 && y < h){
      if (line_mode==0)
	val[ind(x,y,w)] += v;
      else
	val[ind(x,y,w)] = v;	
    }

    if (k >= 4096) { 
      return; }
    if (floor(x)==floor(x1) && floor(y)==floor(y1)) break;
    int e2 = 2*err;
    if (e2 >= dy) {
      err += dy;
      x += sx;
    }
    if (e2 <= dx) {
      err += dx;
      y += sy;
    }
  }
}

int clampi(int f, int a, int b) {
  if (f < a) return a;
  if (f > b) return b;
  return f;
}

double clampd(double f, double a, double b) {
  if (f < a) return a;
  if (f > b) return b;
  return f;
}

void init_img(img_t *i,int n, int w, int h) {
  i->h = h; i->w = w; i->n = n;
  i->p = (double *)malloc(n*w*h*sizeof(double));
}

double sino(double f, double t) {
  return 0.5+0.5*sin(2*M_PI*f*t);
}


//functions for writing ASCII characters

//character map
unsigned int font[768] = {0x00000000,0x00000000,0xff00ff00,0x00000000,0x007e7e00,0x00181800,0xff00ff3e,0x3c3f7f18,0x0081ff6c,0x103c3c00,0xff00ff0e,0x66336318,0x00a5dbfe,0x383c7e00,0xff3cc31a,0x663f7fdb,0x0081fffe,0x7ce7ff18,0xe7669932,0x6630633c,0x0081fffe,0xfee7ff3c,0xc342bd78,0x3c3063e7,0x00bdc3fe,0x7ce77e3c,0xc342bdcc,0x1830633c,0x0099e77c,0x38181818,0xe76699cc,0x7e7067db,0x0081ff38,0x10181800,0xff3cc3cc,0x18f0e718,0x007e7e10,0x007e7e00,0xff00ff78,0x18e0e618,0x00000000,0x00000000,0xff00ff00,0x0000c000,0x00000000,0x00000000,0xff00ff00,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x80021866,0x7f7c0018,0x18180000,0x00000000,0xc0063c66,0xdbc6003c,0x3c180000,0x000010fe,0xe00e7e66,0xdb60007e,0x7e181830,0x002838fe,0xf83e1866,0xdb380018,0x18180c60,0xc06c387c,0xfefe1866,0x7b6c0018,0x1818fefe,0xc0fe7c7c,0xf83e1866,0x1bc60018,0x18180c60,0xc06c7c38,0xe00e7e00,0x1b6cfe7e,0x187e1830,0xfe28fe38,0xc0063c66,0x1b38fe3c,0x183c0000,0x0000fe10,0x80021866,0x1b0cfe18,0x18180000,0x00000000,0x00000000,0x00c6007e,0x00000000,0x00000000,0x00000000,0x007c0000,0x00000000,0x00000000,0x00006600,0x18000030,0x00000000,0x00000000,0x0018666c,0x18003830,0x0c300000,0x00000002,0x003c666c,0x7c006c30,0x18180000,0x00000006,0x003c24fe,0xc6c26c60,0x300c6618,0x0000000c,0x003c006c,0xc0c63800,0x300c3c18,0x00000018,0x0018006c,0x7c0c7600,0x300cff7e,0x00fe0030,0x0018006c,0x0618dc00,0x300c3c18,0x00000060,0x000000fe,0xc630cc00,0x300c6618,0x180000c0,0x0018006c,0x7c66cc00,0x18180000,0x18001880,0x0018006c,0x18c67600,0x0c300000,0x18001800,0x00000000,0x18000000,0x00000000,0x30000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x38187c7c,0x0cfe38fe,0x7c7c0000,0x0c00607c,0x6c38c6c6,0x1cc060c6,0xc6c61818,0x180030c6,0xc6780606,0x3cc0c006,0xc6c61818,0x300018c6,0xc6180c06,0x6cc0c00c,0xc6c60000,0x607e0c0c,0xd618183c,0xccfcfc18,0x7c7e0000,0xc0000618,0xc6183006,0xfe06c630,0xc6060000,0x60000c18,0xc6186006,0x0c06c630,0xc6061818,0x307e1800,0x6c18c6c6,0x0cc6c630,0xc60c1818,0x18003018,0x387efe7c,0x1e7c7c30,0x7c780030,0x0c006018,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x7c10fc3c,0xf8fefe3c,0xc63c1ee6,0xf0c6c638,0xc6386666,0x6c666666,0xc6180c66,0x60eee66c,0xc66c66c2,0x666262c2,0xc6180c6c,0x60fef6c6,0xdec666c0,0x666868c0,0xc6180c6c,0x60fefec6,0xdec67cc0,0x667878c0,0xfe180c78,0x60d6dec6,0xdefe66c0,0x666868de,0xc6180c6c,0x60c6cec6,0xdcc666c2,0x666260c6,0xc618cc6c,0x62c6c6c6,0xc0c66666,0x6c666066,0xc618cc66,0x66c6c66c,0x7cc6fc3c,0xf8fef03a,0xc63c78e6,0xfec6c638,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00001000,0xfc7cfc7c,0x7ec6c6c6,0xc666fe3c,0x803c3800,0x66c666c6,0x7ec6c6c6,0xc666c630,0xc00c6c00,0x66c666c6,0x5ac6c6c6,0x6c668c30,0xe00cc600,0x66c66660,0x18c6c6c6,0x38661830,0x700c0000,0x7cc67c38,0x18c6c6d6,0x383c3030,0x380c0000,0x60d66c0c,0x18c6c6d6,0x38186030,0x1c0c0000,0x60de66c6,0x18c66cfe,0x6c18c230,0x0e0c0000,0x607c66c6,0x18c6387c,0xc618c630,0x060c0000,0xf00ce67c,0x3c7c106c,0xc63cfe3c,0x023c0000,0x000e0000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x000000fe,0x30000000,0x00000000,0x00000000,0x00000000,0x3000e000,0x1c003800,0xe01806e0,0x38000000,0x18006000,0x0c006c00,0x60180660,0x18000000,0x00006000,0x0c006400,0x60000060,0x18000000,0x0078787c,0x3c7c6076,0x6c380e66,0x18ecdc7c,0x000c6cc6,0x6cc6f0cc,0x7618066c,0x18fe66c6,0x007c66c0,0xccfe60cc,0x66180678,0x18d666c6,0x00cc66c0,0xccc060cc,0x6618066c,0x18d666c6,0x00cc66c6,0xccc6607c,0x66180666,0x18d666c6,0x00767c7c,0x767cf00c,0xe63c66e6,0x3cc6667c,0x00000000,0x000000cc,0x00006600,0x00000000,0x00000000,0x00000078,0x00003c00,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x10000000,0x0000000e,0x18707600,0x00000000,0x30000000,0x00000018,0x1818dc00,0x00000000,0x30000000,0x00000018,0x18180010,0xdc76dc7c,0xfccc66c6,0xc6c6fe18,0x18180038,0x66cc76c6,0x30cc66c6,0x6cc6cc70,0x000e006c,0x66cc6670,0x30cc66d6,0x38c61818,0x181800c6,0x66cc601c,0x30cc66d6,0x38c63018,0x181800c6,0x7c7c60c6,0x36cc3cfe,0x6c7e6618,0x181800fe,0x600cf07c,0x1c76186c,0xc606fe0e,0x18700000,0x600c0000,0x00000000,0x000c0000,0x00000000,0xf01e0000,0x00000000,0x00f80000,0x00000000,0x3c000c10,0x00603800,0x10006000,0x1860c638,0x66cc1838,0xcc306c00,0x38cc3066,0x3c30c66c,0xc2cc306c,0xcc183800,0x6ccc1866,0x66181038,0xc0000000,0x0000003c,0x00000000,0x00003838,0xc0cc7c78,0x78787866,0x7c7c7c38,0x38386c6c,0xc0ccc60c,0x0c0c0c60,0xc6c6c618,0x1818c6c6,0xc2ccfe7c,0x7c7c7c66,0xfefefe18,0x1818c6c6,0x66ccc0cc,0xcccccc3c,0xc0c0c018,0x1818fefe,0x3cccc6cc,0xcccccc0c,0xc6c6c618,0x1818c6c6,0x0c767c76,0x76767606,0x7c7c7c3c,0x3c3cc6c6,0x06000000,0x0000003c,0x00000000,0x00000000,0x7c000000,0x00000000,0x00000000,0x00000000,0x18000010,0x00603060,0x00c6c618,0x3800f80e,0x30003e38,0xc6307830,0xc6c6c618,0x6c66cc1b,0x60006c6c,0xc618cc18,0xc638003c,0x6466cc18,0xfe00cc00,0x00000000,0x006cc666,0x603cf818,0x66eccc7c,0x7c7ccccc,0xc6c6c660,0xf018c47e,0x6036fec6,0xc6c6cccc,0xc6c6c660,0x607ecc18,0x787eccc6,0xc6c6cccc,0xc6c6c666,0x6018de18,0x60d8ccc6,0xc6c6cccc,0xc6c6c63c,0x607ecc18,0x66d8ccc6,0xc6c6cccc,0x7e6cc618,0xe618cc18,0xfe6ece7c,0x7c7c7676,0x06387c18,0xfc18c6d8,0x00000000,0x00000000,0x0c000000,0x00000070,0x00000000,0x00000000,0xf8000000,0x00000000,0x180c1818,0x00763c38,0x00000000,0x00000000,0x30183030,0x76dc6c6c,0x30000040,0x40180000,0x60306060,0xdc006c6c,0x300000c6,0xc6180000,0x00000000,0x00c63e38,0x000000cc,0xcc0036d8,0x78387ccc,0xdce60000,0x300000d8,0xd8186c6c,0x0c18c6cc,0x66f67e7c,0x30fefe30,0x3018d836,0x7c18c6cc,0x66de0000,0x60c00660,0x663c6c6c,0xcc18c6cc,0x66ce0000,0xc6c006dc,0xce3c36d8,0xcc18c6cc,0x66c60000,0xc6c00686,0x963c0000,0x763c7c76,0x66c60000,0x7c00000c,0x3e180000,0x00000000,0x00000000,0x00000018,0x06000000,0x00000000,0x00000000,0x0000003e,0x06000000,0x1155dd18,0x18183600,0x00363600,0x36361800,0x44aa7718,0x18183600,0x00363600,0x36361800,0x1155dd18,0x18183600,0x00363600,0x36361800,0x44aa7718,0x18183600,0x00363600,0x36361800,0x1155dd18,0x18f83600,0xf8f636fe,0xf636f800,0x44aa7718,0x18183600,0x18063606,0x06361800,0x1155dd18,0xf8f8f6fe,0xf8f636f6,0xfefef8f8,0x44aa7718,0x18183636,0x18363636,0x00000018,0x1155dd18,0x18183636,0x18363636,0x00000018,0x44aa7718,0x18183636,0x18363636,0x00000018,0x1155dd18,0x18183636,0x18363636,0x00000018,0x44aa7718,0x18183636,0x18363636,0x00000018,0x18180018,0x00181836,0x36003600,0x36003618,0x18180018,0x00181836,0x36003600,0x36003618,0x18180018,0x00181836,0x36003600,0x36003618,0x18180018,0x00181836,0x36003600,0x36003618,0x18180018,0x00181f36,0x373ff7ff,0x37fff7ff,0x18180018,0x00181836,0x30300000,0x30000000,0x1fffff1f,0xffff1f37,0x3f37fff7,0x37fff7ff,0x00001818,0x00181836,0x00360036,0x36003600,0x00001818,0x00181836,0x00360036,0x36003600,0x00001818,0x00181836,0x00360036,0x36003600,0x00001818,0x00181836,0x00360036,0x36003600,0x00001818,0x00181836,0x00360036,0x36003600,0x36000036,0x18000036,0x181800ff,0x00f00fff,0x36000036,0x18000036,0x181800ff,0x00f00fff,0x36000036,0x18000036,0x181800ff,0x00f00fff,0x36000036,0x18000036,0x181800ff,0x00f00fff,0x36ff0036,0x1f1f0036,0xff1800ff,0x00f00fff,0x36000036,0x18180036,0x181800ff,0x00f00fff,0xffffff3f,0x1f1f3fff,0xfff81fff,0xfff00f00,0x00183600,0x00183636,0x180018ff,0xfff00f00,0x00183600,0x00183636,0x180018ff,0xfff00f00,0x00183600,0x00183636,0x180018ff,0xfff00f00,0x00183600,0x00183636,0x180018ff,0xfff00f00,0x00183600,0x00183636,0x180018ff,0xfff00f00,0x00000000,0x00000000,0x00000000,0x00000000,0x0000fe00,0xfe000000,0x7e38381e,0x00031c00,0x0000c600,0xc6000000,0x186c6c30,0x0006307c,0x007cc6fe,0x60006676,0x3cc6c618,0x007e60c6,0x76c6c06c,0x307e66dc,0x66c6c60c,0x7edb60c6,0xdcfcc06c,0x18d86618,0x66fec63e,0xdbdb7cc6,0xd8c6c06c,0x30d86618,0x66c66c66,0xdbf360c6,0xd8c6c06c,0x60d87c18,0x3cc66c66,0x7e7e60c6,0xdcfcc06c,0xc6d86018,0x186c6c66,0x006030c6,0x76c0c06c,0xfe706018,0x7e38ee3c,0x00c01cc6,0x00c00000,0x0000c000,0x00000000,0x00000000,0x00400000,0x00000000,0x00000000,0x00000000,0x00000000,0x00180000,0x3800000f,0xd8700000,0x0000300c,0x0e180000,0x6c00000c,0x6cd80000,0xfe181818,0x1b181800,0x6c00000c,0x6c300000,0x00180c30,0x1b181876,0x3800000c,0x6c607c00,0x007e0660,0x181800dc,0x0000000c,0x6cc87c00,0xfe180c30,0x18187e00,0x0018000c,0x6cf87c00,0x00181818,0x18180076,0x001818ec,0x00007c00,0x0000300c,0x18d818dc,0x0000006c,0x00007c00,0xfe000000,0x18d81800,0x0000003c,0x00007c00,0x00ff7e7e,0x18700000,0x0000001c,0x00000000,0x00000000,0x18000000,0x00000000,0x00000000,0x00000000,0x18000000,0x00000000,0x00000000};

//how much to scale the characters by
#define CHAR_SCALE 1

int getfontpixel(int x, int y) {
  x /= CHAR_SCALE;
  y /= CHAR_SCALE;
  int bit = y*CHAR_WIDTH*16 + x;
  return (font[bit >> 5] & (1 << (31-(bit & 0x1f))));
}
int getfontpixel_noscale(int x, int y) {
  int bit = y*CHAR_WIDTH*16 + x;
  return (font[bit >> 5] & (1 << (31-(bit & 0x1f))));
}


double tex_lookup(double *src,double xxx,double yyy, int w, int h) {
  double x = xxx * w;
  double y = yyy * h;

  double xx = clampd(x-0.5,0,(double)(w-1));
  double yy = clampd(y-0.5,0,(double)(h-1));
	
  int x1 = clampi((int)xx,0,w-1);
  int x2 = clampi((int)xx + 1,0,w-1);

  int y1 = clampi((int)yy,0,h-1);
  int y2 = clampi((int)yy + 1,0,h-1);

  double b1 = src[ind(x1,y1,w)];
  double b2 = src[ind(x2,y1,w)] - src[ind(x1,y1,w)]; 
  double b3 = src[ind(x1,y2,w)] - src[ind(x1,y1,w)]; 
  double b4 = src[ind(x1,y1,w)] - src[ind(x2,y1,w)] - src[ind(x1,y2,w)] + src[ind(x2,y2,w)]; 
	
  double dx = xx - (double)x1;
  double dy = yy - (double)y1;

  double tot =  b1 + b2*dx + b3*dy + b4*dx*dy;

  return tot;  
  
}

//non space varying box blur
void val_boxblur(double rhox,double rhoy, double *val, int w, int h) {
  double a = 0;
  double mx = rhox*2 + 1;
  double my = rhoy*2 + 1;  
  double *dst = make_bufn(w*h);
  
  for (int i=0;i<w;i++) {
    a = 0;
    for (int j=0;j<h;j++) {
      int xx = i;
      int yy= j-my;
      
      a += val[ind(i,j,w)];
      if (!(xx < 0 || xx >= w || yy < 0 || yy >= h)) {
	a -= val[ind(xx,yy,w)];
      }

      yy = j - rhoy;
      if (!(xx < 0 || xx >= w || yy < 0 || yy >= h)) {
	dst[ind(xx,yy,w)] = a;
      }
    }
  }
  lerp(1.0,dst,val,w*h);
  fill(0,dst,w*h);


  for (int j=0;j<h;j++) {    
    a = 0;
    for (int i=0;i<w;i++) {
      int xx = i - mx;
      int yy= j;
      
      a += val[ind(i,j,w)];
      if (!(xx < 0 || xx >= w || yy < 0 || yy >= h)) {
	a -= val[ind(xx,yy,w)];
      }

      xx = i-rhox;
      if (!(xx < 0 || xx >= w || yy < 0 || yy >= h)) {
	dst[ind(xx,yy,w)] = a;
      }
    }
  }
  lerp(1.0,dst,val,w*h);
  fill(0,dst,w*h);
  double fact = 1.0/(mx*my);
  for (int i=0;i<w*h;i++) {
    val[i] *= fact;
  }
    
  free(dst);
}


void blit_rot(int mode, double am, double *v1, double *v2, int sx, int sy, int sw, int sh,
	      int dx, int dy, int dw, int dh, int w, int h, int w2, int h2,
	      double a, double b, double c, double d) {
  
  for (int x=dx;x<dx+dw;x++) {
    for (int y=dy;y<dy+dh;y++) {
      //normalized coordinates
      
      double nx = (double)(x-dx)/dw;
      double ny = (double)(y-dy)/dh;

      nx = nx*sw + sx;
      ny = ny*sh + sy;
      
      //optionally apply rotation/translation matrix
      nx /= w;
      ny /= h;
      nx -= 0.5;
      ny -= 0.5;
      double nnx = a*nx + b*ny + 0.5;
      double nny = c*nx + d*ny + 0.5;

      double av;
      if (am >=0) { 
	av = am*tex_lookup(v1,nnx,nny,w,h);
      } else { //invert if less than 0
	av = am * (tex_lookup(v1,nnx,nny,w,h) -1.0);
      }

      if (mode==1) {
      v2[ind(x,y,w2)] = av;
      } else {
      v2[ind(x,y,w2)] += av;	
      }
    }
  }
}

int instring(char *s, int c) {
  for (int i=0;i<strlen(s);i++) {
    if (s[i]==c) return 1;
  }
  return 0;
}



int map_touches(int *map,int x,int y,int gw,int c) {
  for (int ii=x-1;ii <= x+1;ii++) {
    for (int jj=y-1;jj <= y+1;jj++) {
      int i = ii; int j= jj;
      
      if (ii<0) i = 0;
      if (ii >= gw) i = gw -1;

      if (jj<0) j = 0;
      if (jj >= gw) j = gw -1;

      if (map[ind(i,j,gw)]==c) return 1;
    }
  }
  return 0;
}

int map_clear(int *map,int x,int y,int gw,char *stuff) {
  int num = 0;
  for (int ii=x-1;ii <= x+1;ii++) {
    for (int jj=y-1;jj <= y+1;jj++) {
      int i = ii; int j= jj;
      if (ii<0 || ii >= gw) continue;
      if (jj<0 || jj >= gw) continue;      
      int ch = map[ind(i,j,gw)];
      
      if (ch != 0 && strchr(stuff,ch)) num++;
    }
  }
  return num;
}

void filli(int a, int *env,int n) {
  for (int i=0;i<n;i++) env[i] = a;
}

typedef struct level_t {
  int p[2];
  int du[2]; //door up
  int dd[2]; //door down
  int st; //state
  int goal[2];
  int gw;
  double ts; //timestamp of last state change
  double lts[2048];
  int ln;
} level_t;


int view_radius;
void add_seen(int *seen, level_t *l) {
  int gw = l->gw;
  int radius = view_radius;
  for (int i = l->p[0]-radius;i<= l->p[0]+radius;i++) {
    for (int j = l->p[1]-radius;j<= l->p[1]+radius;j++) {
      if (i < 0 || i >= gw || j < 0 || j >= gw) continue;
      int dx = i-l->p[0];
      int dy = j-l->p[1];
      if (sqrt(dx*dx+dy*dy) >= radius) continue;
      seen[ind(i,j,gw)] = 1;
    }
  }
}
  
int find_dir(int *dir, int *gmap, int x, int y, int gw) {
  dir[0] = 0; dir[1] = 0;
  if (x < 0 || x >= gw || y < 0 || y>= gw) {printf("x: %d y: %d gw: %d",x,y,gw); return -1; }
  
  int c = gmap[ind(x,y,gw)];
  if (c==-1) c = 10000;
  int choice_x[9];
  int choice_y[9];  
  int ci = 0;
  for (int i=x-1;i<=x+1;i++) {
    for (int j=y-1;j<=y+1;j++) {
      if (x < 0 || x >= gw || y < 0 || y>= gw) continue;
      int d = gmap[ind(i,j,gw)];
      if (d == -1) continue;
      //      printf("d: %d c: %d\n",d,c);
      if (d < c) {
	choice_x[ci] = i;
	choice_y[ci++] = j;
      }
    }
  }
  
  if (ci==0) { printf("couldn't find direction\n"); return -1; }

  //when there are multiple choices, pick one randomly
  int i = randp()%ci;
  dir[0] = choice_x[i];
  dir[1] = choice_y[i];
}

void calc_gmap(int *gmap,int *map,int x, int y, int gw) {
  filli(-1,gmap,gw*gw);
  gmap[ind(x,y,gw)] = 0;
  printf("doing gmap at %d %d\n",x,y);
  
  int m = 0;
  while (m < 1000) {
    int flag = 0;
    
    for (int i=0;i<gw;i++) {
      for (int j=0;j<gw;j++) {
	int c = gmap[ind(i,j,gw)];
	if (c != m) continue;
	for (int ii=i-1;ii<=i+1;ii++) {
	  for (int jj=j-1;jj<=j+1;jj++) {
	    if (ii < 0 || ii >= gw) continue;
	    if (jj < 0 || jj >= gw) continue;
	    if (ii==i && jj==j) continue;
	    if (!instring(".><",map[ind(ii,jj,gw)])) continue;
	    int c2 = gmap[ind(ii,jj,gw)];
	    if (c2==-1 || (c+1 < c2)) {
	      gmap[ind(ii,jj,gw)] = c+1;
	      flag=1;	      
	    }

	  }
	}
      }
    }

    if (flag==0) {
      printf("calc_gmap finished with max %d path length\n",m);
      //no changes made
      return;
    }
    m++;
  }
}

double *grid;
int *map,*seen,*inview;
int *gmap[3];

double *font2;
img_t font_img;


level_t l;

void etch_char(unsigned char ch, double fg, int dx, int dy, double os, double *val, int vw,int vh) {
  int fw = CHAR_WIDTH*16;
  int gw = l.gw;
  
  int npr = 16;  //number of characters per row

  int c = (int)ch;

  int cwidth = CHAR_WIDTH;
  int cheight = CHAR_HEIGHT;

  int cy = c/npr;
  int cx = c - cy*npr;

  int sx = cwidth * cx;
  int sy = cheight * cy;

  double dw = (double)vw/gw * 0.9 * (double)cwidth/cheight;
  double dh = (double)vh/gw * 0.9;

  /*
  dx += drand(-1,1)*(double)cwidth * 1e-1;
  dy += drand(-1,1)*(double)cheight * 1e-1;

  dx = limit(0,vw-1,dx);
  dy = limit(0,vh-1,dy);
  */
  
  blit_rot(1,fg,font2,val,sx,sy,cwidth,cheight,
	   dx,dy,dw,dh, fw,193, vw,vh,
	   1,0,0,1);
  
  return;
}


void dragon_curve(int num_levels, double *v, int w) {
  int bytesize = (4 << num_levels)*20;
  char *buf = (char *)malloc(bytesize);
  char *buf2 = (char *)malloc(bytesize);
  
  strcpy(buf,"F");
  buf[1] = 0;

  double rp = 1.979;
  double lp = 1.998;  

  for (int k = 0;k<num_levels;k++) {
    int j =0;
    for (int i=0;i<strlen(buf);i++) {
      if (buf[i]=='F' && drand(0,1.0) < rp)  {
	strcpy(&buf2[j],"F+G");
	j+=3;
      } else if (buf[i]=='G' && drand(0,1.0) < lp) {
	strcpy(&buf2[j],"F-G");
	j+=3;
      } else {
	buf2[j++] = buf[i];
      }
    }
    buf2[j] = 0;
    
    //printf("1. %s\n2. %s\n",buf,buf2);
    strcpy(buf,buf2);
  }

  double ang0 = drand(0,M_PI*2);
  double ang = ang0;
  double dd = 1.0;
  double x = 0,y=0;
  double x2=0,y2=0;

  double minx=x,miny=y,maxx=x,maxy=y;

  int jnum = 1;
  double jang=2*M_PI;
  double jsign[] = {1.0,1.0,1.0,1.0};
  int jw[] = {1,1,1,1};
  int mode = randp()%3;

  printf("mode: %d\n",mode);
  double fl = drand(0,100);
  if (fl < 50) {
    mode =0;
  } else if (fl < 85) {
    mode = 2;
  } else if (fl < 100) {
    mode = 1;
  }
  
  if (mode==1) {
    jang = 2*M_PI*0.5;
    jnum = 4;
  } else  if (mode==2) {
    jang = 2*M_PI*0.25;
    jnum = 4;
    //    jsign[1] = -1; jsign[3]= -1;
  }

  
  for (int j=0;j<jnum;j++) {
    for (int i=0;i<strlen(buf);i++) {
      if (buf[i]=='F' || buf[i]=='G') {
	x2 = x + sin(ang)*dd;
	y2 = y + cos(ang)*dd;
      } else if (buf[i]=='-') {
	ang -= M_PI/2 * jsign[j%4];
      } else if (buf[i]=='+') {
	ang += M_PI/2 * jsign[j%4];
      }
      if (1 || jw[j%4]) {
	if (x2 < minx) minx = x2;
	if (y2 < miny) miny = y2;
    
	if (x2 > maxx) maxx = x2;
	if (y2 > maxy) maxy = y2;
      }
      x = x2; y = y2;
    }
    ang += jang;
  
  }

  double scf = fmax(fabs(maxx-minx), fabs(maxy - miny));
  scf *= 1+1e-4;
  printf("minmax: %.3f %.3f, %.3f %.3f: %.3f\n",minx, maxx, miny, maxy,scf);
  double mm = fmax(minx,miny);
  
  x = (-minx)* (1.0/scf)*w;// + 2;
  y = (-miny)*(1.0/scf)*w;// + 2;

  if (fabs(maxy-miny) > fabs(maxx-minx)) {
    x += ((maxy-miny) - (maxx-minx))/2.0 * (1.0/scf)*w;
  } else {
    y += ((maxx-minx) - (maxy-miny))/2.0 * (1.0/scf)*w;    
  }
  
  dd = (1.0/scf) * w;
  ang = ang0;

  minx = miny = maxx = maxy = 0;
  for (int j=0;j<jnum;j++) {
    for (int i=0;i<strlen(buf);i++) {
      if (buf[i]=='F' || buf[i]=='G') {
	x2 = x + sin(ang)*dd;
	y2 = y + cos(ang)*dd;
      } else if (buf[i]=='-') {
	ang -= M_PI/2* jsign[j%4];
      } else if (buf[i]=='+') {
	ang += M_PI/2* jsign[j%4];
      }
      //    val_line(x,y,x2,y2,1.0,val[0],w,h);
      if (x2 < -5) 
	printf("x2: %.3f y2: %.3f\n",x2,y2);
      if (1 || jw[j%4]) {      
	val_line(x,y,x2,y2,1.0,v,w,w);

	if (x2 < minx) minx = x2;
	if (y2 < miny) miny = y2;
    
	if (x2 > maxx) maxx = x2;
	if (y2 > maxy) maxy = y2;

      }
      x = x2; y = y2;
    }
    ang += jang;
  
  }
  printf("after minmax: %.3f %.3f, %.3f %.3f: %.3f\n",minx, maxx, miny, maxy,scf);
  
  free(buf);
  free(buf2);
}

void do_water(int *map, int ch, int num) {
  int choice[l.gw*l.gw];
  int ci = 0;
  int i = 0;
  int k = 0;

  char stuff[3];
  stuff[0] = '.';
  stuff[1] = ch;
  stuff[2] = 0;
  
  while(k++ < 1000) {
    i = randp()%(l.gw*l.gw);
    int x,y;
    dni(i,&x,&y,l.gw);
    if (map_clear(map,x,y,l.gw,stuff)==9) break;
  }
  if (k >= 1000) return;

  map[i] = ch;
  choice[ci++] = i;

  for (int j=0;j<num;j++) {
    i = choice[randp()%ci];
    int x,y;
    dni(i,&x,&y,l.gw);
    int x2 = x;
    int y2 = y;
    
    do {
      if (randp()%2==0) { x2 = x + (randp()%2==0? 1 : -1); }
      else { y2 = y + (randp()%2==0 ? 1: -1);}
    } while (x2 < 0 || x2 >= l.gw || y2 < 0 || y2 >= l.gw);

    
    if (map_clear(map,x2,y2,l.gw,stuff) != 9) continue;
    map[ind(x2,y2,l.gw)] = ch;
    if (ci > l.gw*l.gw-1) break;;
    choice[ci++] = ind(x2,y2,l.gw);
  }
}

void newlevel(int gw, int num_levels) {
  l.gw = gw;
  printf("Making newlevel\n");
  
  if (grid) free(grid);
  if (map) free(map);
  if (seen) free(seen);
  if (inview) free(inview);
  
  grid = make_bufn(gw*gw);
  
  map = (int *)calloc(gw*gw,sizeof(int));
  int *map2 = (int *)calloc(gw*gw,sizeof(int));  
  seen = (int *)calloc(gw*gw,sizeof(int));
  inview = (int *)calloc(gw*gw,sizeof(int));  

  for (int i=0;i<3;i++) {
    if (!gmap[i]) free(gmap[i]);
    gmap[i] = (int *)calloc(gw*gw,sizeof(int));  
  }

  int goff = 4;
  
  dragon_curve(num_levels,grid,gw-goff);
  printf("Done Dragon curve: %d %d\n",num_levels,gw-goff);
  //translate fractal grid to map
  for (int x = 0;x<gw-goff;x++) {
    for (int y = 0;y<gw-goff;y++) {
      if (grid[ind(x,y,gw-goff)] > 0) {
	map[ind(x+goff/2,y+goff/2,gw)] = '.';
      }
    }
  }

  for (int x = 0;x<gw;x++) {
    for (int y = 0;y<gw;y++) {  
      if (map[ind(x,y,gw)]==0 && map_touches(map,x,y,gw,'.')) {
	map[ind(x,y,gw)] = '#';
      }
    }
  }

  //2nd pass
  int pass2 = randp()%2;
  
  if (pass2) {
    int dx = 0,dy=0,xflip = 0,yflip = 0,xyflip = 0;

    if (randp()%2==0)
      dx = randp()%3-1;
    else
      dy = randp()%3-1;
    
    xflip = randp()%2;
    yflip = randp()%2;
    xyflip = randp()%2;

    xflip = 0;
    yflip = 0;
    xyflip = 0;
    
    for (int x = 0;x<gw;x++) {
      for (int y = 0;y<gw;y++) {
	int x2=x,y2=y;
	if (xflip) x2 = gw -x2;
	if (yflip) y2 = gw -y2;
	x2 += dx;
	y2 += dy;
	if (x2 < 0 || x2 >= gw) continue;
	if (y2 < 0 || y2 >= gw) continue;

	if (xyflip) {
	  int tmp = x2;
	  x2 = y2;
	  y2 = tmp;
	}
	map2[ind(x2,y2,gw)] = map[ind(x,y,gw)];
      }
    }
    
    for (int x = 0;x<gw;x++) {
      for (int y = 0;y<gw;y++) {
	if (map2[ind(x,y,gw)] != 0 && strchr("><",map2[ind(x,y,gw)])) continue;
	    
	if (map2[ind(x,y,gw)]=='.') {
	  map[ind(x,y,gw)]='.';
	} else if (map2[ind(x,y,gw)]=='#' && map[ind(x,y,gw)]==0) {
	  map[ind(x,y,gw)] = choose3('/','-','\\');
	}
      }
    }

    for (int x = 0;x<gw;x++) {
      for (int y = 0;y<gw;y++) {  
	if (map[ind(x,y,gw)]==0 && map_touches(map,x,y,gw,'.')) {
	  map[ind(x,y,gw)] = choose3('/','-','\\');	
	}
      }
    }
  }

  int pass3 = 0;
  double proll = drand(0,100);
  if (proll < 25) {
    pass3 = 2;
  } else if (proll < 55) {
    pass3 = 1;
  } else {
    pass3 = 0;
  }
  
  if (pass3) {
    for (int x = 0;x<gw;x++) {
      for (int y = 0;y<gw;y++) {  
	if (map[ind(x,y,gw)]=='#') {
	  switch(pass3) {
	  case 0:
	    break;
	  case 1:
	    map[ind(x,y,gw)] = choose2(choose2('-','|'),choose2('/','\\'));	    
	    break;
	  case 2:
	    map[ind(x,y,gw)] = choose2(choose2('*','+'),choose2('O','*'));
	    break;
	  }
	}
      }
    }
  }

  int num_water = 2+randp()%8;
  num_water = ceil(num_water*(double)gw/100.0);

  for (int i=0;i<num_water;i++) {
    do_water(map,'~',gw*gw*(5000/100));
  }

  for (int i=0;i<2*(gw*gw)/250;i++)
    do_water(map,'^',5);

  for (int i=0;i<2*(gw*gw)/250;i++)
    do_water(map,'v',5);

  //Place door down at random location
  int flag = 0;
  do {
    int x = randp()%gw;
    int y = randp()%gw;
    if (map[ind(x,y,gw)] =='.') {
      map[ind(x,y,gw)] = '>';
      l.dd[0] = x; l.dd[1] = y;
      flag = 1;
    }
  } while (!flag);

  //Place door up at random location
  flag = 0;
  do {
    int x = randp()%gw;
    int y = randp()%gw;
    if (map[ind(x,y,gw)] =='.') {
      map[ind(x,y,gw)] = '<';
      l.p[0] = x; l.p[1] = y; //player starts at door up location
      l.goal[0] = x; l.goal[1] = y;
      l.du[0] = x; l.du[1] = y;      
      flag = 1;
    }
  } while (!flag);

  add_seen(seen,&l);
  free(map2);
}

double calc_move_dur(int _gw) {
  double gw = (double)_gw;
  double r = 0;
  if (gw < 40)  r =  gw/20;
  else if (gw < 60)  r =  gw/22;  
  else if (gw < 80) r =  gw/25;
  else if (gw < 100) r =  gw/29;  
  else if (gw < 120) r =  gw/34;
  else if (gw < 140) r =  gw/40;
  else r = gw/47;
  r *= 6;
  return r;
}

int audio_only=0;
int pic_only = 0;

void visual(double dur) {
  int w2 = 1080;
  int h2 = 1080;
  
  randp_setstream(1);
  
  printf("Generating visual, %d x %d..\n",w2,h2);

  int cw = 24;
  int ch = 16;
  int cn = cw*ch;
  
  double *cha[4];
  double *cval[4];
  
  for (int i=0;i<4;i++) {
    cha[i] = make_bufn(cn);
    cval[i] = make_bufn(cn);    
  }

  printf("w: %d, h: %d\n",w2,h2);
  
  double os=1; //oversampling rate during processing
  double dos=1; //sampling rate before rendering  (leave at one for downsample from 'os' value)
  int w = w2*os;
  int h = h2*os;
  
  img_t img;
  //  init_img(&img,3,w2,h2);
  init_img(&img,3,w2*dos,h2*dos);

  int n = w*h;
  double *val[6]; //buffer storing images
  for (int i=0;i<6;i++) {
    val[i] = make_bufn(n); 
  }


  //----------------------------------------

  font2 = make_bufn(16*16*CHAR_WIDTH*CHAR_HEIGHT);
  for (int x = 0;x < 16*CHAR_WIDTH;x++)
    for (int y = 0;y < 16*CHAR_HEIGHT;y++)
      font2[ind(x,y,16*CHAR_WIDTH)] = getfontpixel(x,y) ? 1.0: 0.0;
  
  double t = 0;
  
  int fn = 0;
  
  /*  write_frame(16*CHAR_WIDTH,16*CHAR_HEIGHT,os/dos,font2,&img,&fn,&t);
  exit(0);  
  */
  grid = 0;
  map = seen = inview = 0;
  for (int i=0;i<3;i++) gmap[i] = 0;
    
  //initialize new level
  int _nl = 11;
  int _gw = 35;
  _nl = 15; _gw = 35;
    //_nl = 10; _gw = 45;  
  //17, 180
  
  view_radius = 7;

  //initialize level_t struct
  l.ln = 0;

  newlevel(_gw,_nl);
    
  /*
  write_frame(_gw-4,_gw-4,os/dos,grid,&img,&fn,&t);
  exit(0);
  */
  l.st = 0; //coming up from door
  l.ts = 0;
  
  /*
  dragon_curve(7,val[0],w);
  write_frame(w,h,os/dos,val[0],&img,&fn,&t);  
  exit(0);
  */
  
  
  double move_dur = (double)_gw/100.0 * 2;
  move_dur = calc_move_dur(_gw);
  double next_ts = move_dur;

  double move_counter = 0;

  double fi_length = 1.2;
  double fo_length = 1.3;

  int door_anim_ctr =0;
  
  while (t <= dur) {
    fill(0,val[0],n);
    double *v = val[1];

    if (l.st==0) {
      //fade in
      if (t < l.ts+fi_length)	goto dodraw;
      printf("st: 1 t: %.3f\n",t);            
      l.st =1;
    } else if (l.st==1) {
      //we do regular moving
    } else if (l.st==2) {
      //fade out

      if (t < l.ts+fo_length) goto dodraw;

      fill(0,val[1],n);
      fill(0,val[2],n);            
      fill(0,val[3],n);
      
      printf("st: 0 t: %.3f\n",t);      
      l.st=0;
      l.ts = t;
      l.lts[l.ln++]=t;
	
      _gw *= drand(1.2,1.4);
      _nl += 1;
      if (_nl > 15) _nl = 15;
      if (_gw > 150) _gw = 35;
	
      newlevel(_gw,_nl);
      move_dur = calc_move_dur(_gw);
      goto dodraw;
    }
    
    move_counter += (move_dur/24.0  * (1 + drand(0,0.25)));
    
    while (move_counter > 1.0) {
      move_counter-= 1.0;
      add_seen(seen,&l);

      if (l.p[0] == l.goal[0] && l.p[1] == l.goal[1] &&
	  l.p[0] == l.dd[0] && l.p[1] == l.dd[1]) {
	//we've arrived at  the door, transition to fade-out state
	printf("st: 2 t: %.3f\n",t);            	
	l.st=2;
	l.ts = t;
	
	goto dodraw;
      }

      /*
      if (seen[ind(l.dd[0],l.dd[1],gw)])  {
	//if we've discovered the door, make it new goal
	calc_gmap(gmap[0],map,l.dd[0],l.dd[1],gw);
      }
      */
      
      if (l.p[0] == l.goal[0] && l.p[1] == l.goal[1]) {
	printf("we've arrived at goal, calculate a new one\n");
	
	int choice[l.gw*l.gw];
	int ci = 0;
	for (int i=0;i<l.gw*l.gw;i++) {
	  if ((map[i] != 0 && strchr(".#",map[i])) && seen[i]==0) {
	    choice[ci++] = i;
	  }
	}

	if (ci==0) {
	  //we've explored everywhere, go to the door
	  printf("proceeding to the door\n");
	  l.goal[0] = l.dd[0]; l.goal[1] = l.dd[1];
	} else {
	  long m = 1e10;
	  int x=0; int y =0;
	  
	  for (int k=0;k<100;k++) {
	    dni(choice[randp()%ci],&x,&y,l.gw);
	    int dx = x - l.p[0];
	    int dy = y - l.p[1];
	    if (dx*dx + dy*dy < m) {
	      m = dx*dx + dy*dy;
	      l.goal[0] = x; l.goal[1] = y;
	    }
	  }
	}
	calc_gmap(gmap[0],map,l.goal[0],l.goal[1],l.gw);
      }
      
      // perform movement            
      int dir[2];
      //  printf("player beginging at: %d %d\n",l.p[0],l.p[1]);
      
      if (find_dir(dir, gmap[0],l.p[0],l.p[1],l.gw) >= 0) {
	if (map[ind(dir[0],dir[1],l.gw)] =='#') {
	  l.goal[0] = l.p[0];
	  l.goal[1] = l.p[1];	  
	} else {
	  l.p[0] = dir[0];
	  l.p[1] = dir[1];
	}
      } else {
	//there was a problem, maybe abandon goal?
	printf("abandoning goal\n");
	seen[ind(l.goal[0],l.goal[1],l.gw)] = 1;
	if (l.p[0] >= l.gw) l.p[0] = l.gw -1;
	if (l.p[1] >= l.gw) l.p[1] = l.gw -1;
	
	l.goal[0] = l.p[0];
	l.goal[1] = l.p[1];
	
      }
      // add discovery
      add_seen(seen,&l);

      // calculate what's in view
      filli(0,inview,l.gw*l.gw);
      add_seen(inview,&l);
    } //done move


    //now draw
  dodraw:

    //    calc_gmap(gmap[0],map,20,12,l.gw);
	
    for (int x = 0;x<l.gw;x++) {
      for (int y =0;y<l.gw;y++) {
	double xx = (double)x/l.gw * w;
	double yy = (double)y/l.gw * h;	

	int ch =map[ind(x,y,l.gw)];

	//draw player
	if (x == l.p[0] && y == l.p[1]) {
	  if (l.st==2 || l.st==0)  {
	    char buf[] = "\\/-H";
	    door_anim_ctr++;
	    etch_char(buf[(door_anim_ctr>>1)%4], 1.0,xx,yy,os,v,w,h);
	    continue;
	  }
	  etch_char('H', 1.0,xx,yy,os,v,w,h);
	  continue;
	}

	
	if (!pic_only && seen[ind(x,y,l.gw)]==0) continue; //only draw parts of map that have been seend
	int iv = inview[ind(x,y,l.gw)]; //check if in view (TODO upgrade this to line of sight)
	double br = iv ? 0.8 : 0.1;

	double dx = x - l.p[0];
	double dy = y - l.p[1];
	if (iv) br = limit(0.1,0.8, pow(1.0 - limit(0,1.0,sqrt(dx*dx + dy*dy)/(double)view_radius),1.0));


	/*
	if (ch =='.') {
	  ch = gmap[0][ind(x,y,l.gw)];	  
	  etch_char('a'+ch, 1.0,xx,yy,os,v,w,h);
	  continue;
	}
	*/
	
	if (ch==0) {
	  //	  etch_char(&font_img, '#', 0.4,xx,yy,os,v,w,h);	  
	} else if (ch=='.') {
	  etch_char('.', br + drand(0.1,0.15),xx,yy,os,v,w,h);

	  /*
	  square3( xx,yy,
		   1.0*1.0/(_nl*0.8), //size factor
		   t,dur,v,w,h);
	  */

	}else if (ch=='>' || ch=='<') {
	  etch_char(ch, br + drand(0.2,0.4),xx,yy,os,v,w,h);
	} else if (ch=='H') {
	  etch_char(ch, br,xx,yy,os,v,w,h);
	} else if (strchr("#oO0*",ch)) {
	  etch_char(ch, br + drand(0.15,0.3),xx,yy,os,v,w,h);
	} else if (strchr("/|-\\",ch)) {
	  etch_char(ch, br + drand(0.3,0.45),xx,yy,os,v,w,h);
	} else if (strchr("^v~",ch)) {
	  etch_char(ch, br + drand(0.15,0.20),xx,yy,os,v,w,h);	  
	}
      }
    }

    
    //    scale(1.0/nib,0,v,n);
    
    add(v,val[2],n);
    double sa = 0.91;
    
    //blur radius
    double ba = drand(0.01,0.03);
    ba = 0.035;
    
    if (0 && randp()%18==0) { //TODO change back to 12
      //skip boxblur scaling on 1/12 of frames
    } else {
      scale(sa,0,val[2],n);
      val_boxblur(w*ba,h*ba,val[2],w,h);
    }

    //noise1(0.2,0.9,val[2],w,h);    
    lerp(1.0,val[2],val[3],n);
    add(v,val[3],n);
    normalize(val[3],n);    

    softclip(1.8,0.0,val[3],n);
    
    //    noise1(0.1,0.8,val[3],w,h);
    
    //val_boxblur(w*ba*0.1,h*ba*0.1,val[3],w,h);
    //    lerp(1.0,v,val[3],n);
    
    //    val_lens2(-0.05,val[3],w,h);
    //    val_boxblur_radial(w*0.02,h*0.02,val[3],w,h);
    
    normalize(val[3],n);

    double a = 1.0;



    if (l.st==0) {
      //fade in
      a = (t - l.ts)/fi_length;
    } else if (l.st==2) {
      //fade out
      a = 1.0 - (t - l.ts)/fo_length;
    }

    if (!pic_only)
      scale(a,0,val[3],n);      

    int num_exp = 1;
    
    for (int i=0;i<num_exp;i++) {
      write_frame(w,h,os/dos,val[3],&img,&fn,&t);
      if (pic_only) exit(0);
    }
  }

  //list of timestamps
  l.lts[l.ln++] = t;
  for (int i=0;i<l.ln;i++) {
    printf("%d: %.3f\n",i,l.lts[i]);
  }

}

//res2a with time dependent coefficients, bwenv is in ocatves
void res2a_bset(bset *fenv,
		bset *bwenv,
		double fmult,
		double *data,int n) {

  bset fenv2,bwenv2;
  bset_init(&fenv2);
  bset_add(&fenv2,0,440);
  bset_add(&fenv2,100000,440);

  bset_init(&bwenv2);
  bset_add(&bwenv2,0,1.0);
  bset_add(&bwenv2,100000,1.0);    

  if (fenv==0) fenv = &fenv2;
  if (bwenv==0) bwenv = &bwenv2;

  double x[3] = {0,0,0}; double y[3] = {0,0,0};

  for (int i=0;i<n;i++) {
    double t = (double)i/SAMPLE_RATE;
    double fc = bset_interp(fenv,t) * fmult;        
    if (t > bset_length(fenv)) {
      fc = bset_last(fenv) * fmult;
    }

    double bw = bset_interp(bwenv,t) * fc; //bwenv is in octaves

    if (fc <= 10) {
      // printf("fc is: %.3f\n",fc);  //TODO
    }
    double C = -1.0 * exp(-2.0*M_PI*bw/SAMPLE_RATE);
    double B = 2.0 * exp(-1.0*M_PI*bw/SAMPLE_RATE) * cos(2*M_PI*fc/SAMPLE_RATE);
    double A = 1.0 - B - C;

    x[0] = x[1];  x[1] = x[2]; x[2] = data[i];
    y[0] = y[1]; y[1] = y[2];
    y[2] = A * x[2] + B * y[1] + C *y[0];
    data[i] = y[2];
  }

  bset_destroy(&fenv2);
  bset_destroy(&bwenv2);
  return;
  
  iirfilter_params rf;
  rf.ftype = FTYPE_RES2A;
  rf.fc = bset_interp(fenv,0);
  rf.bw = bset_interp(bwenv,0);

  init_iirfilter(&rf);
  
  for (int i=0;i<n;i++) {
    double t = (double)i/n;
    
    double x = data[i];
    
    rf.fc = bset_interp(fenv,t);
    rf.bw = bset_interp(bwenv,t);
    //rf.bw = t*bwenv[ib] + (1.0-t)*rf.fc;
    iir_update_params(&rf);

    x = iir_apply(&rf,x);

    data[i] = x;
  }
}

int only_metadata;


void chick(double *tone,int n) {
  double *env = make_bufn(n);
  double *t = make_bufn(n);  

  fill(0,tone,n);

  double f = 1300;
  
  whitenoise(1.0,t,n);

  expad(.015,0.14,env,n);
  mult(env,t,n);
  res2a(f,2000,t,n);
  res2a(f*1.5,2000,t,n);
  //  res2a(f*2.3,1000,t,n);  
  bandpass(800,5000,t,n);
  bandpass(800,5000,t,n);  
  //  res2a(f*3,300,t,n);    
  add(t,tone,n);

  whitenoise(1.0,t,n);
  bandpass(600,1100,t,n);
  bandpass(700,900,t,n);
  expad(.01,0.15,env,n);
  mult(env,t,n);
  dbgain(-5,t,n);
  res2a(f,3000,t,n);
  res2a(f*2,3000,t,n);    
  add(t,&tone[(int)(n*0.22)],n);

  highpass(2800,tone,n);
  highpass(1000,tone,n);  
  free(env);
  free(t);  
}  



typedef struct {
  double *Xr2;
  double *Xi2;
  double *cos2pidN;
  double *sin2pidN;
  int *ind;
  int N;
} fft_plan;
void release_fft_mem(fft_plan *f) {
  if (f->N != 0) {
    if (f->Xr2) free(f->Xr2);
    if (f->Xi2) free(f->Xi2);

    if (f->cos2pidN) free(f->cos2pidN);
    if (f->sin2pidN) free(f->sin2pidN);
    if (f->ind) free(f->ind);
  }  
}
void prepare_fft(fft_plan *f,int N) {
  if (f->N == N) {
    return;
  }
  if (f->N != 0) {
    if (f->Xr2) free(f->Xr2);
    if (f->Xi2) free(f->Xi2);

    if (f->cos2pidN) free(f->cos2pidN);
    if (f->sin2pidN) free(f->sin2pidN);
    if (f->ind) free(f->ind);
  }
  printf("Preparing fft %d\n",N);
  f->N = N;
  f->Xr2 = (double *)malloc(N * sizeof(double));
  f->Xi2 = (double *)malloc(N * sizeof(double));

  f->cos2pidN = (double *)malloc(N * sizeof(double));
  f->sin2pidN = (double *)malloc(N * sizeof(double));
  f->ind = (int *)malloc(N * sizeof(int));
  
  //precompute cos/sin table
  for (int i=0;i<N;i++) {
    f->cos2pidN[i] = cos(2*M_PI*(double)i/N);
    f->sin2pidN[i] = -sin(2*M_PI*(double)i/N);    
  }
  
  //build index
  f->ind[0] = 0;
  for (int n= N; n>1; n /=2) {
    for (int i=0;i<N;i+=n) {
      f->ind[i+n/2] = f->ind[i] + N/n;
    }
  }
}

void fft(double *x, double *Xr, double *Xi,fft_plan *f,int N) {
  double *X = Xr;
  if (N != f->N) {
    prepare_fft(f,N);
  }
  
  double *Xr2 = f->Xr2;
  double *Xi2 = f->Xi2;

  //compute first layer
  for (int i=0;i<N;i+=2) {
    Xr[i] = x[f->ind[i]] + x[f->ind[i+1]];
    Xr[i+1] = x[f->ind[i]] - x[f->ind[i+1]];
    
    Xi[i] = 0; 
    Xi[i+1] = 0;
  }
  
  double *Er,*Or,*Ei,*Oi;
	     
  for (int n = 2;n<N;n*=2) {
    int CI = (N/(n*2));
    for (int j = 0;j<N;j+=n*2) {
      Er = &Xr[j];
      Or = &Xr[j+n];

      Ei = &Xi[j];
      Oi = &Xi[j+n];

      double *Xr2ij=&Xr2[j];
      double *Xi2ij=&Xi2[j];
      double *Xr2ijn=&Xr2[j+n];
      double *Xi2ijn=&Xi2[j+n];

      double *sint = f->sin2pidN;
      double *cost = f->cos2pidN;
      for (int i=0;i<n;i++,sint+=CI,cost+=CI) {
	//compute n*2 point DFT
	//1st half

	*Xr2ij++ = *Er + (*Or*(*cost) - *Oi*(*sint));
	*Xi2ij++ = *Ei + (*Or*(*sint) + *Oi*(*cost));	

	//2nd half (TODO: omit in last level)

	*Xr2ijn++ = *Er++ - (*Or*(*cost) - *Oi*(*sint));
	*Xi2ijn++ = *Ei++ - (*Or++ *(*sint) + *Oi++ *(*cost));
      }
    }
    
    //now swap pointers (TODO: omit in last level)
    double *tmp = Xr;Xr = Xr2;Xr2 = tmp;
    tmp = Xi;Xi = Xi2; Xi2 = tmp;
  }

  // copy to passed array if required
  if (X != Xr) {
    double *tmp = Xr;Xr = Xr2;Xr2 = tmp;
    tmp = Xi;Xi = Xi2; Xi2 = tmp;
    //copy Xr2
    for (int j=0;j<N;j++) { Xr[j] = Xr2[j]; Xi[j] = Xi2[j];
    }
  }
}

void ifft(double *x, double *Xr, double *Xi,fft_plan *f,int N) {
  double *X = Xr;
  if (N != f->N) {
    prepare_fft(f,N);
  }
  double *Xr2 = f->Xr2;
  double *Xi2 = f->Xi2;
  
  //compute first layer
  for (int i=0;i<N;i+=2) {
    Xr2[i] = Xr[f->ind[i]] + Xr[f->ind[i+1]];
    Xr2[i+1] = Xr[f->ind[i]] - Xr[f->ind[i+1]];

    Xi2[i] = Xi[f->ind[i]] + Xi[f->ind[i+1]];
    Xi2[i+1] = Xi[f->ind[i]] - Xi[f->ind[i+1]];
  }

  //now swap pointers (TODO: omit in last level)
  double *tmp = Xr;Xr = Xr2;Xr2 = tmp;
  tmp = Xi;Xi = Xi2; Xi2 = tmp;
    
  double *Er,*Or,*Ei,*Oi;
	     
  for (int n = 2;n<N;n*=2) {
    int CI = (N/(n*2));
    for (int j = 0;j<N;j+=n*2) {
      Er = &Xr[j];
      Or = &Xr[j+n];

      Ei = &Xi[j];
      Oi = &Xi[j+n];

      double *Xr2ij=&Xr2[j];
      double *Xi2ij=&Xi2[j];
      double *Xr2ijn=&Xr2[j+n];
      double *Xi2ijn=&Xi2[j+n];

      double *sint = f->sin2pidN;
      double *cost = f->cos2pidN;
      for (int i=0;i<n;i++,sint+=CI,cost+=CI) {
	//compute n*2 point DFT
	//1st half

	*Xr2ij++ = *Er + (*Or*(*cost) + *Oi*(*sint));
	*Xi2ij++ = *Ei + (*Or*(-*sint) + *Oi*(*cost));	

	//2nd half (TODO: omit in last level)

	*Xr2ijn++ = *Er++ - (*Or*(*cost) + *Oi*(*sint));
	*Xi2ijn++ = *Ei++ - (*Or++ *(-*sint) + *Oi++ *(*cost));
      }
    }

    if (n == N/2 && X == Xr) {
      //don't swap
    }else {
      //now swap pointers (TODO: omit in last level)
      double *tmp = Xr;Xr = Xr2;Xr2 = tmp;
      tmp = Xi;Xi = Xi2; Xi2 = tmp;
    }
  }

  // copy to passed array if required
  if (X != Xr) {
    double *tmp = Xr;Xr = Xr2;Xr2 = tmp;
    tmp = Xi;Xi = Xi2; Xi2 = tmp;
    //copy Xr2
    for (int j=0;j<N;j++) { Xr[j] = Xr2[j]; Xi[j] = Xi2[j]; }
  }
  
  for (int j=0;j<N;j++) { x[j] = Xr[j]/(double)N;}
}

//convolution using fft
void convfft(double *src,int src_n,double *dst,int dst_n) {
  int N = 1;
  while (src_n+dst_n > N) N*=2;
  N*=2;
  fft_plan f;
  f.N = 0;  

  printf("Performing FFT with N = %d\n",N);
  printf("Src: %d, DST: %d\n",src_n,dst_n);
  
  double *srcx = (double*)malloc(N * sizeof(double));
  double *SRCr = (double*)malloc(N * sizeof(double));
  double *SRCi = (double*)malloc(N * sizeof(double));

  fill(0,srcx,N);
  add(src,srcx,src_n);
  
  fft(srcx,SRCr,SRCi,&f,N);
  printf("Completed FFT\n");
  double *dstx = (double*)malloc(N * sizeof(double));
  double *DSTr = (double*)malloc(N * sizeof(double));
  double *DSTi = (double*)malloc(N * sizeof(double));

  fill(0,dstx,N);
  add(dst,dstx,dst_n);

  fft(dstx,DSTr,DSTi,&f,N);

  for (int i=0;i<N;i++) {
    double Xr = DSTr[i];
    double Xi = DSTi[i];

    DSTr[i] = Xr * SRCr[i] - Xi*SRCi[i];
    DSTi[i] = Xr * SRCi[i] + Xi*SRCr[i];    
  }

  ifft(dstx,DSTr,DSTi,&f,N);
  for (int i=0;i<dst_n;i++) {
    dst[i] = dstx[i];
  }
  release_fft_mem(&f);
  free(srcx); free(SRCr); free(SRCi);
  free(dstx); free(DSTr); free(DSTi);
}

//filtered white noise exponential trail
//g gain in db
void reverb(double dur, double g, double predelay, double *src,double *dst,int n) {
  int dn;
  double *buf = make_buf(dur*2.5 + predelay,&dn);
  double *env = make_bufn(dn);

  expenv(dur/6,env,dn);
  if (predelay > 1e-6) delay((int)(predelay*SAMPLE_RATE),env,dn);

  whitenoise(1.0,buf,dn);
  //origianl 100hz to 1500hz
  lowpass(2800,buf,dn); //for brighter reverb
  highpass(200,buf,dn);
  highpass(200,buf,dn);
  mult(env,buf,dn);
  
  lerp(1.0,src,dst,n);
  //  add(src,dst,n);
  convfft(buf,dn,dst,n);
  dbgain(g,dst,n);
  
  free(buf);
  free(env);
}

//perform convolution reverb in windowed chunks
void stereochunkreverb(double dur, double g, double predelay, double *src0,double *src1, double *dst0,double *dst1,
			double chunk_dur,int n) {
  int cn = chunk_dur * SAMPLE_RATE;
  int ov = 1.0 * SAMPLE_RATE; //1 second overlap  
  double *buf = make_bufn(cn+ov*2+dur*2.5);
  double *buf2 = make_bufn(cn+ov*2+dur*2.5);  
  
  int num_chunks = n/cn + 1;
  for (int i=0;i<num_chunks;i++) {
    printf("FFT reverb chunk %d of %d\n",i,num_chunks);

    for (int j=0;j<2;j++) {
      double *src,*dst;
      if (j==0) { src = src0; dst = dst0; }
      if (j==1) { src = src1; dst = dst1; }      
      fill(0,buf,cn+ov*2+dur*2.5);
      fill(0,buf2,cn+ov*2+dur*2.5);      
      
      if (i==0) {
	//first chunk
	lerp(1.0,src,buf,cn+ov);
	nfadeout(ov*2,buf,cn+ov);
	reverb(dur,g,predelay,buf,buf2,cn+ov*2+dur*2.5);
	add(buf2,dst,cn+ov*2+dur*2.5);
      } else if (i==num_chunks-1) {
	//last chunk

	int max_remaining = n - (i*cn-ov);
	max_remaining--;
	printf("last chunk: insert point is %.3f\n",(double)(i*cn-ov)/SAMPLE_RATE);
	printf("last chunk: max dur is %.3f\n",(double)max_remaining/SAMPLE_RATE);
	
	lerp(1.0,&src[i*cn-ov],buf,limit(0,max_remaining,cn+ov));
	nfadein(ov*2,buf,cn+ov);
	reverb(dur,g,predelay,buf,buf2,cn +ov*2+dur*2.5);
	add(buf2,&dst[i*cn-ov],limit(0,max_remaining,cn+ov*2+dur*2.5));
      } else {
	//middle chunk      
	lerp(1.0,&src[i*cn-ov],buf,cn+ov*2);
	nfadein(ov*2,buf,cn+ov*2);
	nfadeout(ov*2,buf,cn+ov*2);
	reverb(dur,g,predelay,buf,buf2,cn+ov*2+dur*2.5);
	add(buf2,&dst[i*cn-ov],cn+ov*2+dur*2.5);

      }
    }
  }
  printf("Freeing chunk reverb window memory.\n");
  free(buf);
  free(buf2);
  printf("Done freeing chunk reverb window memory.\n");  
}

void thwack(double dur, double tf, double *dst) {
  int p = SAMPLE_RATE*dur;
  double *env = make_bufn(p);
  whitenoise(1.0,dst,p);
  res2a(tf, tf*0.8,dst,p);
  expad(.002,dur/6,env,p);
  mult(env,dst,p);
  free(env);
}

//sin tone with time varying frequency envelope
void sin_tone_bset(bset *f,double *data,int n) {
  gen_data g;

  init_gen(&g);
  g.amp = 1.0;
  double t =0;
  double y = bset_interp(f,t);
  //  double slew = .03;
  
  for (int i=0;i<n;i++) {
    double dy = bset_interp(f,t) - y;
    //    if (dy > slew) dy = slew;
    //    if (dy < -slew) dy = -slew;
    y += dy;
    g.freq = y;
    
    data[i] = sin_gen(&g); //TODO: change this to a function pointer
    t += 1.0/SAMPLE_RATE;
  }
}


void audio(double dur) {
  printf("Generating audio..\n");
  randp_setstream(2);  
  
  int tone_n = (dur+20.0) * SAMPLE_RATE;  

  double *tone[5];
  double *bar[5];  
  for (int i=0;i<5;i++) {
    tone[i] = make_bufn(tone_n);
    fill(0,tone[i],tone_n);

    bar[i] = make_bufn(tone_n);
    fill(0,bar[i],tone_n);    
  }

  int env_n = (dur+20.0) * SAMPLE_RATE;
  double *env[5];
  for (int i=0;i<5;i++) {
    env[i] = make_bufn(env_n);
    fill(0,env[i],env_n);
  }

  int chan_n;
  double *chan[2];
  double *al[2];  

  chan[0] = make_buf(dur*2.1,&chan_n);
  chan[1] = make_bufn(chan_n);
  al[0] = make_bufn(chan_n);
  al[1] = make_bufn(chan_n);
  
  double *efx[2];
  efx[0] = make_bufn(chan_n);
  efx[1] = make_bufn(chan_n);
    
  int p =0;

  double t =0;
  p = dur*SAMPLE_RATE;



  //------------------------------------------------------------------------


  
  bset ff,amp;
  bset ffb;
  bset_init(&ff);
  bset_init(&amp);  
  bset_init(&ffb);

  double tempo = 4;
  double f = drand(100,2000);
  f = 900;
  double hf = 3200;

  double bw = 0.004;
  bset bwenv;
  bset_init(&bwenv);
  bset_segment(&bwenv,0,10000,bw,bw);
  
  double bdur = 0;
  drand(0,1.0);
  
  for (int i=0;i<10;i++) {
    
    bset_clear(&ff);
    bset_clear(&amp);
    f = drand(200,1000);
    hf = limit(1200,20000,f*drand(2.0,3.7));
    t = drand(0,3);
    bset_add(&ff,0,f);

    bdur = drand(0.1,0.13*1000.0/f);
    
    fill(0,tone[0],p);
    fill(0,tone[1],p);
    fill(0,tone[2],p);
    double ddur = drand(4.0,8);
    
  while (t < dur) {

    bset_add(&ff,t,f*(1+0.03*drand(-1,1)));
    double aa = 0.6 + drand(-1,1)*0.05;
    bset_add(&ff,t+bdur*drand(0.25,0.35),aa*hf+(1-aa)*f);
    
    bset_add(&ff,t+bdur*0.9,hf);
    bset_add(&ff,t+ddur*0.99,hf);        

    bset_add(&amp,t,0);

    bset_add(&amp,t+drand(0.008,0.012),1.0);
    

    bset_add(&amp,t+drand(0.037,0.049),drand(0.18,0.23));

    bset_add(&amp,t+bdur*0.5,0.0);

    int ci = t*SAMPLE_RATE;
    
    //thwack sound
    thwack(bdur,drand(88,92),&tone[1][ci]);
    //    thwack(bdur,drand(2800,3200),&tone[2][ci]);
    thwack(bdur,drand(6800,7200),&tone[2][ci]);        
  
    t += ddur;
  }

  p = dur*SAMPLE_RATE;
  sin_tone_bset(&ff,tone[0],p);


  whitenoise(1.0,tone[0],p);
  res2a_bset(&ff,&bwenv,1,tone[0],p);

  /*
  lowpass(2000,tone[0],p);
  lowpass(4000,tone[0],p);
  lowpass(6000,tone[0],p);  
  */
  apply_amp(&amp,tone[0],p);

  dbgain(drand(-2,2),tone[1],p);
  add(tone[1],tone[0],p);
  
  dbgain(drand(-17,-10),tone[2],p);
  add(tone[2],tone[0],p);  

  highpass(300,tone[0],p);
  highpass(300,tone[0],p);  
  
  double panf = drand(-1,1.0);

  double gg = limit(-30,0,-3*i);
  gg += drand(-3,3);
  
  dbgain(gg,tone[0],p);
  
  msmix(tone[0],0,panf,chan[0],chan[1],p);

  }

  whitenoise(1.0,tone[0],p);
  res2a(80,100,tone[0],p);
  highpass(200,tone[0],p);
  lowpass(5000,tone[0],p);
  
  //  smoothlines_lfo(1.0/80,1.0/1000,1.0,0,tone[0],p);
  msmix(tone[0],-4,0,chan[0],chan[1],p);
  
  double dtime = 0.3;
  ffdelay(0,dtime,-10,200,2800,chan[0],efx[0],p);
  ffdelay(0,dtime,-10,200,3000,chan[1],efx[1],p);
  int ddtime = dtime*0.5*SAMPLE_RATE;
  delay(ddtime,efx[1],p-ddtime);
  dbgain(-10,efx[0],p);
  dbgain(-10,efx[1],p);  
  add(efx[0],chan[0],p);
  add(efx[1],chan[1],p);


  fill(0,efx[0],p);
  fill(0,efx[1],p);
  
  //calcualte reverb (optional)
  int use_reverb = 1;
  
  if (use_reverb) {
    if (dur > 51) {
      //windowed convolution reverb when the audio track is very long
      double window_size = 50.0;
      printf("Calculating FFT convolution reverb in windows of 50.0 s...\n");
      stereochunkreverb(3.0,-43,.02, chan[0],chan[1],efx[0],efx[1],50.0,p);
      normalize_stereo(efx[0],efx[1],p);
      dbgain(13,efx[0],p);
      dbgain(13,efx[1],p);     
    } else {
      printf("Calculating FFT convolution reverb...\n");        
      //for 50s or less, crunch reverb all at once
      reverb(3.0,-30, .02, chan[0],efx[0],p);
      reverb(3.0,-30, .02, chan[1],efx[1],p);
    }
  
    add(efx[0],chan[0],p);
    add(efx[1],chan[1],p);
  }

  normalize_stereo(chan[0],chan[1],p);
  stereomixdown(chan[0],chan[1],p);

  free(efx[0]);
  free(efx[1]);
  for (int i=0;i<5;i++) {free(env[i]); }
}

unsigned long long seed2long(char *seed,int i) {
  char sbuf[9];
  memcpy( sbuf, &seed[2+i*8], 8 );
  sbuf[8] = '\0';
  unsigned long long v = strtol(sbuf,NULL,16);
  return v;
}


int defaultSeed() {
  if (strcmp(seed,"0x2c68ca657f166e1e7a5db2266ee1edca9daf2da5ecc689c39c1d8cfe0b9aad2e")==0) {
    return 1;
  }
  return 0;
}

#define rand() myrand()

int myrand() {
  for (int i=0;i<20;i++) {
    printf("WARNING!!!! do not use rand();\n");
  }
  exit(0);
  return 0;
}

int main(int argc, char *argv[]) {
  init_token_params();

  //duration of 7 minutes by default. Can be changed with command line option '-d'
  double dur = 60*7;
  
  int visual_only=0;
  int ex = 0;
  
  if (argc > 1 && strcmp(argv[1],"-a")==0) audio_only = 1;
  if (argc > 1 && strcmp(argv[1],"-v")==0) {printf("Option: visual only\n"); visual_only = 1; }
  
  if (argc > 1 && strcmp(argv[1],"-a")==0) {
    audio_only = 1;
    if (argc > 2) {
      ex = atoi(argv[2]);
      printf("using extra seed ex: %d\n",ex);            
    }
    
    printf("Option: Audio only\n");    
  }

  if (argc > 1 && strcmp(argv[1],"-pi")==0) {
    
    if (argc > 2) {
      ex = atoi(argv[2]) * 13847234;
      printf("using extra seed ex: %d\n",ex);      
    }

    pic_only = 1;
    printf("Option: export still image only");
  }
  
  if (argc > 2 && strcmp(argv[1],"-d")==0) {
    dur = atof(argv[2]);
    dur = limit(10,60*10,dur);
    printf("Option: Duration of %.3f seconds\n",dur);
  }

  if (argc > 3 && strcmp(argv[2],"-d")==0) {
    dur = atof(argv[3]);
    dur = limit(10,60*10,dur);
    printf("Option: Duration of %.3f seconds\n",dur);
  }

  if (argc > 4 && strcmp(argv[3],"-d")==0) {
    dur = atof(argv[4]);
    dur = limit(10,60*10,dur);
    printf("Option: Duration of %.3f seconds\n",dur);
  }
  
  if (argc > 4 && strcmp(argv[3],"-d")==0) {
    dur = atof(argv[4]);
    dur = limit(10,60*10,dur);
    printf("Option: Audio and visual, with duration of %.3f seconds\n",dur);
  }
  
  // seed 8 PRNG streams using token_param[0]
  for (int i=0;i<8;i++) {
    unsigned long long s = 0;
    s = 298375823498234 + ex + i + token_param[0];
    randp_seed(i,s); 
  }

  if (!audio_only)
    visual(dur); //compute visual

  //audio uses PRNG stream 2
  if (!visual_only)
    audio(dur); //compute audio
  
  return 0;
}