/*
CHRONOPHOTOGRAPH

Chronophotograph is the companion collection to Vol. 2, Series 1 - Noumenon. This work is made in support of The Los Angeles County Museum of Art, presented in dialogue with the pioneering chronophotopgraphic work of Edweard Muybridge from the LACMA archives.

On a time locked schedule, a blockchain transaction can be triggered to metaphorically "capture" an observation of a Noumenon, minting a new ERC721 token in this companion series. The blocknumber serves as the seed for a deterministic program to generate an image in the style of Muybridge's chronophotographs, representing the attempt to transcend sensory limitations through the use of technology.

Viewers are encouraged to reflect on the perception of time, consensus reality, photography and blockchain as markers of time, sources of objective truth, limitations thereof.

by 0xDEAFBEEF
May 2023

This self contained C file will compile and run to generate 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. Fill in the values in init_token_params() with parameters retrieved
from Ethereum smart contract at 0xda1bf9b5de160cecde3f9304b187a2f5f5b83707

2. Compile and run:
# gcc -Ofast main.c -lm && ./a.out

This will produce:
a. PGM(Portable gray map) file named 'image.pgm'

16bit grayscale, linear (no gamma)

*/

typedef struct {
long long blocknumber;
int noumenon_token_id;
int chronophotograph_token_id;
int epoch;
int res;
int cat;
} tokenparams_struct_t;
tokenparams_struct_t token_params;

// Fill in the following parameters with return values from Etheereum smart contract
// 0xda1bf9b5de160cecde3f9304b187a2f5f5b83707 function getTokenParams(tokenID)

void init_token_params() {
// insert token ID...
token_params.chronophotograph_token_id = 0;

// ...and the parameters returned from smart
// contract function getTokenParams(tokenID)

//blocknumber the chronophotograph was minted
token_params.blocknumber = 0;

//ID of parent noumenon token. 'parent_id' var in getTokenParams return vals
token_params.noumenon_token_id = 0;

//epoch
token_params.epoch = 0;
}


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

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


//indexing
uint64 ind(uint64 x, uint64 y,uint64 w);
void dni(int i, int *x, int *y,int w);

//for pixel rendering
void val_line(int x0,int y0, int x1,int y1, double v, double *val,int w, int h);
void val_line_list(int x0,int y0,
int x1,int y1,
double v,
int *pa,int *pk, int w, int h);


void val_radial(int w,int h,double p, double ang2, double *val);
void val_rect(double x0, double y0, double x1, double y1, double v, double *val, int w, int h);
void val_lines(double *p,int n,double *val,int w, int h);
void val_arc(double x, double y, double r, double ang1, double ang2,double *val,int w,int h);
void val_poly(double *p,int n,double *val,int w, int h);
void val_invert(double *val,int w,int h);


//lines and drawing
void init_drawline();
void drawline(int x0,int y0,
int x1,int y1,
double v,
double *ch,double *val,int w, int h);
void drawpoly(double *p,int n,double *ch,double *val,int w, int h);
void drawlines(double *p,int n,double *ch,double *val,int w, int h);
void drawarc(double x, double y, double r, double ang1, double ang2,double *ch,double *val,int w,int h);

//masks
void replace_with_mask(double *src, double *dst, double *mask, int n);


void noisefield(double *v, int w, int h);

int param_only = 0;

#define SAMPLE_RATE 44100

typedef struct abuf {
int n; // number of samples
int sr; //sample rate
double *data; //pointer to samples
} abuf;


//breakpoint set
typedef struct bset {
int destroy; //flag to indicate it should be destroyed at end of function
int n; // number of samples used
int max_n; //max allocated samples
int sr; //sample rate
int interp; //default interpolation
double param; //default param
double *t; //pointer to samples
double *y; //pointer to samples
double *iparam; //parameter for interpolation technique (gamma, tau, etc);

double *p; //parameters, 8 per index
int *pi; //index into parameters

int *itype; //type of interpolation for next segment
} bset;

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

//collection of notes, 0 is middle c
typedef struct ncol {
int n;
int d[128];
} ncol;

typedef struct {
const char *name;
int n; //number of pitches/steps
int type; //type of scale, either pitches or semitones
double const * const p; //list pitches per octave relative to a fundamental, maximum 100
int s[12]; //relative list of intervals p steps, evenly tempered scale
} scale_pattern_t;

#define SP_STEP 0
#define SP_PITCH 1


typedef struct {
double freq;
double phase;
double amp;
double pw; //pulse width for pulse generator
int i; // sample counter

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

//for noisetone
double hp_filt;
double t;
double tt;
} gen_data;


typedef struct {
int ftype;
double fc; //center freq
double bw; //bandwidth (for resonent filters)

double pole;
double a1;
double a2;
double rA; // for RES2A filter
double rB; // for RES2A filter
double rC; // for RES2A filter
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


typedef struct {
const char *name;
double F[6]; //formants
double BW[6]; //formant bandwidths
double A[6]; // noise filter params
double AB; //wide band fricative noise amplitude

double AV; //voicing amplitude
double AVS; //lowpassed voicing amplitude

double AF; //fricative noise amplitude
double AH; //aspiration noice amplitude
double FNZ; //nasal zero
double FNP; //nasal pole
} speech_params_t;


//FIR filter of max length 1024
#define FIR_BUF_MAX 1024

//memory allocation
double *make_buf(double dur,int*n);
double *dup_buf(double *src,int n);
double *make_bufn(uint64 n);


//functions
void fill(double a, double *env,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);
double maxsample(double *data,int n);
void normalize(double *data,int n);
void normalize_stereo(double *lchan, double *rchan,int n);
void removedc(double *data,int n);
void dbgain(double db,double*data,int n);
void gain(double g,double*data,int n);
void reverse(double *data,int n);
void scale(double a,double b, double *dst, int n);
void delay(int b, double *data,int n);
void shift(int b, double *data,int n);
double linear_interp(double p, double *data,int n);
double linear_interp_circular(double p, double *data,int n);
double l2norm(double *data, int n);
void slewlimit(double slew, double *data,int n);

double triangle_naive(double w);
void scrub(double *pos,double *src, int n);


//compression, gating, envelope detection
void envdetect(double attack, double release, double *src,double *dst,int n);
void compress(double ratio, double threshold, double *cs, double *src,double *dst,int n);
void compress_inplace(double attack, double release, double db, double *src,int n);
void compress_inplace_stereo(double attack, double release, double db, double *src0, double *src1, int n);
void trigger_seq(double *env, int env_n, double *trig, int trig_n, double *chan);
void trigger_bset(double *env, int env_n, bset *b,double *chan);
int itrigger(int p, double *data, int n);
void beatdetect(double min_dur, double db_over_rms, double *src,bset *b, int n);

//distort
void softclip(double gain,double offset,double *data,int n);
void hardclip(double gain,double *data,int n);
void softclip_t(double *env, int env_n, double offset,double *data,int n);

//fm modulation
void modulate(double A, double freq, double *src, double src_n, double *dst);
void modulate_t(double A, double *freq, double *src, double src_n, double *dst);
void modulate_t2(double *A, double *freq, double *src, double src_n, double *dst);

//random
double drand ( double low, double high );
double exprand(double lambda);
double nrand();
int irand ( int low, int high );
double choose2(double c1,double c2);
double choose3(double c1,double c2,double c3);
double exprand(double lambda);
double tnrand(double mean, double sd, double low, double high);
double tnrands(double mean,double sd, int n);
long randp();
void randp_setstream(int i);
unsigned long long randp_getseed(int i);
int randp_seed(int i,unsigned long long entropy);
void randp_reset(int i);

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

//filters
void lowpass(double fc,double *data,int n);
void reslp(double fc,double q, double *data,int n);
void reslp_t(double *fenv, double env_n,double q, double *data,int n);
void reslp_bset(bset *fenv,double q, double *data,int n);

void lpg_t(double *fenv, double env_n,double *data,int n);
void highpass(double fc,double *data,int n);
void bandpass(double fl,double fh,double *data,int n);
void res2a(double fc,double bw, double *data,int n);
void vowel(int v, double *data,int n);
void res2(double fc,double q, double *data,int n);
void nres2a(double fc,double bw, double *data,int n);
void res2a_t(double *fenv, double env_n,
double *bwenv, double bwenv_n,
double *data,int n);

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

//envelopes
void expad(double a,double d, double *env, int env_n);
void expadsr(double a,double d, double s, double r, double dur, double *env, int env_n);
void expenv(double tc, double *env, int env_n);
void smoothlines_lfo(double dur1,double dur2,double amp,double offset,double *data,int n);
void sawtooth_lfo(double freq,double amp,double offset,double *data,int n);
double pulse_naive(double w,double pw);
double sawtooth_naive(double w);
void sin_lfo(double freq,double amp,double offset,double *data,int n);
void sin_lfo_ph(double freq,double phase,double amp,double offset,double *data,int n);
void pulse_lfo(double freq,double pw, double amp,double offset,double *data,int n);
void triangle_lfo(double freq,double amp,double offset,double *data,int n);
double sin_gen(gen_data *g);
double triangle_gen(gen_data *g);
double pulse_gen(gen_data *g);
double saw_gen(gen_data *g);



void itrapz(int a, double *env, int n);
void trapz(int a, double *env, int n);
void line(double a,double b,double dur,double *env,int env_n);
void linen(double a,double b,double *env,int n);
void logline(double a,double b,double dur,double *env,int env_n);

void fadeinout(double dur1, double dur2, double *src, int n);
void fadein(double dur, double *src,int n);
void fadeout(double dur, double *src,int n);


//generators
void init_gen(gen_data *g);
double sinlu(double w);
void sin_tone(int freq,double start,double dur,double amp,double *data);
void saw_tone(int freq,double start,double dur,double amp,double *data);
void pulse_tone(int freq,double pw, double start,double dur,double amp,double *data,int data_n);

void sin_tone_t(double *fenv, double env_n, double start, double dur,double amp, double *data,int n);
void sin_tone_bset(bset *f,double *data,int n);
void cos_tone_bset(bset *f,double *data,int n);
void saw_tone_bset(bset *f,double *data,int n);
void triangle_tone_bset(bset *f,double *data,int n);
void pulse_tone_bset(bset *f,bset *pw_bset,double *data,int n);

void saw_tone_t(double *fenv, double env_n, double start, double dur,double amp, double *data,int n);
void pulse_tone_t(double *fenv, double env_n,
double *pwenv,int pwenv_n,
double start, double dur,double amp, double *data,int n);
void triangle_tone_t(double *fenv, double env_n,
double start, double dur,double amp, double *data,int n);

void whitenoise(double gain, double *data,int n);
void bassdrum(double freq1,double freq2, double dur,double *dst, int dst_n);
void snare(double freq1, double *dst, int dst_n);
void circular_membrane_drum(double freq,double *dst,int dst_n );
void dsf_tone_bset(bset *fenv,
bset *fmenv,
bset *wenv,
double *data, int n);

void dsf_bass(double f, double dur,double *chan[], double t);

//resonators
void guitar_harmonics(double freq,double *data, int n);
void guitar_harmonics_t(double *fenv,double *data, int n);
void resbank(double freq,double *bw, double *amp, int bn, double *data, int n);


//bsets
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_cubic(bset *b, double t);
double bset_interp_samplehold(bset *b, double t);
double bset_interp_exp(bset *b, double t,double tau);
double bset_interp_gamma(bset *b, double t,double gamma);
double bset_render(bset *b, double *data);
void bset_dump(bset *b);
void bset_write(bset *b,char *filename);
void bset_read(bset *b,char *filename);
void bset_add(bset *b, double t, double y);
void bset_remove(bset *b, int j);
void bset_clear(bset *b);

void bset_inertia(bset *b, double a);
void bset_parse(bset *b, char *s);
void bset_parse2(bset *b, char *s);

void bset_stretch(bset *b,double f);
void bset_clip(bset *b,double t1,double dur);
void bset_shift(bset *b,double f);
void bset_delay(bset *b,double f);
double bset_length(bset *b);
double bset_last(bset *b);
void bset_dup(bset *b, double times, double gap);
void bset_replacem(bset *b, char *p);
bset *bp(char *s);
bset *bramp(double t0, double t1, double v0, double v1);
void bset_blit(bset *src, bset *dst, double t1, double t2, double dur1, double dur2);
void bset_scale(bset *b, double s, double t);

double bset_rendern(bset *b, double *data,int n);
void bset_addm(bset *b, char *p);
void bset_clean(bset *b);
void bset_erase(bset *b, double t1, double t2);
void bset_splice(bset *b, double t1, double t2);
void bset_place(bset *b, bset *b2, double t);
void bset_flatten(bset *b);

void bset_copy(bset *b, bset *dst);


//time stretching
int playback_t(double *speed,int fn,double *src,int src_n,double *dst,int dst_n);
int playback_bset(bset *speed,double *src,int n);
int playback_bset_stereo(bset *speed,double *src1, double *src2,int n);

//delay,reverb
void ffdelay(double dur, double feedback, double fl, double fh, double *src, double*dst, int n);
void tapedelay(bset *speed, double dur, double feedback, double fl, double fh, double *src, double *dst,int n);
void reverb(double dur, double g, double predelay, double *src,double *dst,int n);

//mixdown, stereo/mono, panning
void tomono(double *lchan, double *rchan, double *mono, int n);
void mmmix(double *src,double dbgain, double *dst,int n);
void msmixdown(double*mono,double *l, double *r, int n);
void msmixp(double *src,double dbgain, bset *pos_bset, double *dstl, double *dstr,int n);
void msmix(double *src,double dbgain, double pos, double *dstl, double *dstr,int n);
void monomixdown(double *chan,int n);
void stereomixdown(double *lchan,double*rchan,int n);
void pan_t(double *src,double *pos,int pn, double *dstl, double *dstr,int n);

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

//plotting
void plot(double *data,int n);
void spectrum(double *data,int data_n);

//file read/writing
int read_stereo_wav_file(char *filename,double **lchannel,double **rchannel,int max_samples);
void write_raw(double *data, char *filename, int n);
void read_raw(double *data, char *filename, int n);


//index into 2d array
uint64 ind(uint64 x, uint64 y,uint64 w) {
return y*w + x;
}
//reverse index
void dni(int i, int *x, int *y,int w) {
*y = i/w;
*x = i - *y*w;
}

//chooses random character from a string
unsigned char choosechar(char *s) {
int n = strlen(s);
return s[rand()%n];
}



void val_invert(double *val,int w,int h) {
for (int x=0;x<w;x++) { for (int y=0;y<h;y++) {
double *v = &val[ind(x,y,w)];
*v = 1.0 - *v;
}
}
}

void writepgm(double *v,char *filename, uint64 w, uint64 h) {
FILE *outfile = fopen(filename, "wb");
fprintf(outfile, "P5\n");
fprintf(outfile, "# pixel intensities stored as 16 bit linear (no gamma)\n");
fprintf(outfile, "%llu %llu\n65535\n",w,h);

for (uint64 i=0;i<w*h;i++) {
double a = v[i];
if (a < 0) a = 0;
if (a > 1.0) a = 1.0;
unsigned int x = a * 65535;
fprintf(outfile, "%c%c", (x & 0xFF00) >> 8, x & 0xFF);
}
fclose(outfile);
}

void writebmp2(double *i,char *filename, int w, int h) {
unsigned int WIDTH = w;
unsigned int HEIGHT= h;

unsigned int headers[13];
FILE * outfile;
unsigned int extrabytes;
unsigned int paddedsize;
int x; int y; int n;
int red, green, blue;

extrabytes = 4 - ((WIDTH * 3) % 4); // How many bytes of padding to add to each
// horizontal line - the size of which must
// be a multiple of 4 bytes.
if (extrabytes == 4)
extrabytes = 0;

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

// Headers...
// Note that the "BM" identifier in bytes 0 and 1 is NOT included in these "headers".

headers[0] = paddedsize + 54; // bfSize (whole file size)
headers[1] = 0; // bfReserved (both)
headers[2] = 54; // bfOffbits
headers[3] = 40; // biSize
headers[4] = WIDTH; // biWidth
headers[5] = HEIGHT; // biHeight

// Would have biPlanes and biBitCount in position 6, but they're shorts.
// It's easier to write them out separately (see below) than pretend
// they're a single int, especially with endian issues...

headers[7] = 0; // biCompression
headers[8] = paddedsize; // biSizeImage
headers[9] = 0; // biXPelsPerMeter
headers[10] = 0; // biYPelsPerMeter
headers[11] = 0; // biClrUsed
headers[12] = 0; // biClrImportant

outfile = fopen(filename, "wb");

//
// Headers begin...
// When printing ints and shorts, we write out 1 character at a time to avoid endian issues.
//

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);
}

// These next 4 characters are for the biPlanes and biBitCount fields.

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);
}

//
// Headers done, now write the data...
//


unsigned char *data = (unsigned char*)malloc(WIDTH*3);

for (y = HEIGHT - 1; y >= 0; y--) // BMP image format is written from bottom to top...
{
for (x = 0; x <= WIDTH - 1; x++)
{
double a = i[ind(x,y,w)];
if (a < 0) a = 0;
if (a > 1.0) a = 1.0;

data[x*3] = a*255;
data[x*3+1] = a*255;
data[x*3+2] = a*255;
}
fwrite(data,1,WIDTH*3,outfile);

if (extrabytes) // See above - BMP lines must be of lengths divisible by 4.
{
for (n = 1; n <= extrabytes; n++)
{
fprintf(outfile, "%c", 0);
}
}
}

fclose(outfile);
return;
}

//PARSE
typedef enum {BEGIN, COMMAND, COMMENT} parse_state;
#define PARSE_CB_LENGTH 1024

typedef struct {
parse_state state;
double tvals[1024*2];
double yvals[1024*2];
char buf[1024];
char cbuf[PARSE_CB_LENGTH]; //circular buffer for past characters
int cbp;
int bp;
int voice;
int selected_voice; //voice we wish to filter for
int n;
double t;
double v;
double cur_dur;
double cur_octave;
double transpose;
double beat_dur;
double barline_t;
int verbose;
int beats_per_bar;
int last_was_whitespace;
} parse_t;


//functions
void parse_init(parse_t *p);

void parse_str(parse_t *p,char *c);
void parse_str2(parse_t *p,char *c);

#define PARSE_LOG(p,...) sprintf(logbuf,__VA_ARGS__); parse_log(p,logbuf)
double *make_bufn(uint64 n);

char logbuf[1024*4];

void parse_init(parse_t *p) {
p->verbose = 1;
p->transpose = 0; //number of semitones to transpose
p->t = 0;
p->v = 0;
p->voice = 0;
p->selected_voice = 0;
p->state=0;
p->n = 0;
p->bp = 0;
p->cbp = 0;
p->last_was_whitespace = 1;
for (int i=0;i<PARSE_CB_LENGTH;i++) p->cbuf[i] = 0; // zero circular delay buffer
p->barline_t = 0;
p->beats_per_bar = 4;
p->beat_dur = 0.25;
p->cur_dur = p->beat_dur;
p->cur_octave = 4;
}

int whitespace(char c) {
return (c==' ' || c=='\n' || c=='\t') ? 1 : 0;
}

void parse_log(parse_t *p,char *s) {
if (p->verbose) {
printf("%s",s);
}
}
//read character that occured n characters in the past
char pastchar(parse_t *p, int n) {
int i = p->cbp - n;
if (i < 0) i += PARSE_CB_LENGTH;
return p->cbuf[i];
}
void parse_putc(parse_t *p,char c) { p->buf[p->bp++] = c;}
void parse_rbuf(parse_t *p) { p->bp = 0;}
void parse_cr(parse_t *p) { parse_putc(p,0); parse_rbuf(p); }

//
void parse_command(parse_t *p) {
char *argv[10];
char *s = strtok(p->buf," ");
PARSE_LOG(p,"# command: %s\n",p->buf);
if (s != NULL) {
s++;
if (!strcmp(s,"ts")) {
//time signature command
char *a1 = strtok(NULL," ");
if (a1==NULL) { PARSE_LOG(p,"# error: ts command\n"); return; }
int b = atoi(a1);
if (b > 0) {p->beats_per_bar = b; } else {
PARSE_LOG(p,"# error: ts command is 0\n"); return; }
PARSE_LOG(p,"# time signature, beats_per_bar=%d\n",p->beats_per_bar);
return;
}

if (!strcmp(s,"beat_dur")) {
//time signature command
char *a1 = strtok(NULL," ");
if (a1==NULL) { PARSE_LOG(p,"# error: beat_dur command\n"); return; }
double b = atof(a1);
if (b > 1e-6) { p->beat_dur = b; } else {
PARSE_LOG(p,"# error: beat_dur command is 0\n"); return; }
PARSE_LOG(p,"# beat_dur=%.3f\n",p->beat_dur);
return;
}

if (!strcmp(s,"voice")) {
//voice command
char *a1 = strtok(NULL," ");
if (a1==NULL) { PARSE_LOG(p,"# error: voice command\n"); return; }
int b = atoi(a1);
p->voice = atoi(a1);
PARSE_LOG(p,"# voice=%d\n",p->voice);
return;
}

if (!strcmp(s,"tr")) {
//transpose command
char *a1 = strtok(NULL," ");
if (a1==NULL) { PARSE_LOG(p,"# error: tr\n"); return; }
int b = atoi(a1);
p->transpose = atoi(a1);
PARSE_LOG(p,"# transpose=%.3f\n",p->transpose);
return;
}
}
PARSE_LOG(p,"# unknown command: %s\n",p->buf);
}

//after white space detected, parse token
void parse_token(parse_t *p) {
if (p->voice != p->selected_voice) return; // skip tokens for voices other than the one we're processing
double dur = 0;
double val = 0;
int rest_flag = 0;
int move_flag = 0;
int absolute_flag = 0;
p->bp = 0;
if (strlen(p->buf)==0) return;
char v_str[256];
char t_str[256];
t_str[0] = 0;
char *token = strtok(p->buf,":");
if (token==NULL) {
// no value!
return;
}

if (token != NULL) {
strcpy(v_str,token);
int l = strlen(v_str);
char c = v_str[0];
int notnote_flag = 0;
switch(c) {
case 'a': val = 220; break;
case 'b': val = 247; break;
case 'c': val = 262; break;
case 'd': val = 294; break;
case 'e': val = 330; break;
case 'f': val = 349; break;
case 'g': val = 392; break;
case 'r': rest_flag = 1; break;
case 'm': move_flag = 1; break;
default: notnote_flag = 1;break;
}
if (notnote_flag) {
val = atof(v_str);
}
if (!notnote_flag && l > 1) {
c = v_str[1];
if (c=='#') { val *= pow(2,1/12.0); }
if (c=='&') { val /= pow(2,1/12.0); }
if (strchr("1234567890",c)) p->cur_octave = c - 48;

if (l > 2) {
c = v_str[2];
if (strchr("1234567890",c)) p->cur_octave = c - 48;
}
}
val *= pow(2,p->cur_octave - 4);
val *= pow(2,p->transpose/12.0);
}

token = strtok(NULL,":");
if (token==NULL) {
//no duration info, use saved duration
dur = p->cur_dur;
} else if (token != NULL) {
int l = strlen(token);
double units;
if (l >= 2 && token[l-1]=='b') {
//measured in bars
units = p->beat_dur * p->beats_per_bar;
} else if (l >=3 && token[l-2]=='m' && token[l-1]=='s') {
//measured in milliseconds
units = 1e-3;
} else if (l >=2 && token[l-1]=='s') {
//measured in seconds
units = 1.0;
} else {
//measured in beats
units = p->beat_dur;
}

//duration was supplied, parse it
dur = atof(token) * units;
// if (dur < 1e-10) dur = p->cur_dur;
}

//check for optional 3rd part of token to indicate other options
token = strtok(NULL,":");
if (token != NULL) {
if (token[0]=='a') {
//absolute flag
absolute_flag = 1;
}
}

PARSE_LOG(p,"# token: %s: %s_%s val: %.3f dur: %.3f at time %.3f\n",p->buf,v_str,t_str,val,dur,p->t);
if (move_flag==1) {
if (absolute_flag==1) {
p->t = dur;
p->barline_t = p->t; //when we move, reset barline position
PARSE_LOG(p,"# move to %.3f\n",p->t);
} else {
p->t += dur;
p->barline_t = p->t; //when we move, reset barline position
PARSE_LOG(p,"# move %s %.3f to %.3f\n",(dur > 0 ? "forward" : "back"),dur,p->t);
}
return;
}

if (rest_flag) val = 0;

if (absolute_flag) {
//move to absolute position before printing token
p->t = dur;
PARSE_LOG(p,"# absolute pos %.3f\n",p->t);
PARSE_LOG(p,"%.3f\t%.3f\n",p->t,val);
p->tvals[p->n] = p->t;
p->yvals[p->n] = val;
p->n++;
return;
}

PARSE_LOG(p,"%.3f\t%.3f\n",p->t,val);
p->tvals[p->n] = p->t;
p->yvals[p->n] = val;
p->n++;

p->v = val;
p->cur_dur = dur;
p->t += dur;
}

void parse_1char(parse_t *p, char c) {
switch(p->state) {
case BEGIN: //beginning state
if (whitespace(c)) { parse_cr(p); parse_token(p); break; }
if (c=='|') {
//bar line
parse_cr(p); parse_token(p);
if (p->voice != p->selected_voice) break;
//set time to next bar
p->t = p->barline_t + p->beats_per_bar * p->beat_dur;
p->barline_t = p->t; //time at beginning of this bar
PARSE_LOG(p,"# barline at %.3f\n",p->t);
break;
}

if (c=='/' && p->last_was_whitespace) { p->state = COMMAND; parse_putc(p,c); break; }

if (c=='#' && p->last_was_whitespace) { p->state = COMMENT; parse_putc(p,c); break; }
// if (c=='/' && pastchar(p,1)=='/' && whitespace(pastchar(p,2))) { p->state = COMMENT; parse_putc(p,c); break; }
parse_putc(p,c);
break;
case COMMAND:
// command beginning with :
if (c=='\n') { p->state = BEGIN; parse_cr(p); parse_command(p); break; }
parse_putc(p,c);
break;
case COMMENT:
if (c=='\n') {
p->state = BEGIN;
parse_cr(p);
PARSE_LOG(p,"# comment: %s\n",p->buf);
break; }
parse_putc(p,c);
break;
default:
break;
}
//other state flags
p->last_was_whitespace = whitespace(c);

//buffer of past characters
p->cbuf[p->cbp++] = c;
if (p->cbp >= PARSE_CB_LENGTH) p->cbp = 0;
}

void parse_str(parse_t *p,char *c) {
do {
parse_1char(p,*c++);
} while (*c != 0);
parse_1char(p,'\n');
}

void parse_str2(parse_t *p,char *s) {
//stack for storing sequence lengths
int stack[256];
int st_ptr = 0;
double durstack[1024];
int dst_ptr = 0;

int max_notes = 1024;

double *dur_list =make_bufn(max_notes);
double *val_list =make_bufn(max_notes);


int ml_ptr = 0;
int seq_length = 0;
double seq_dur = 0;

int octave=4;
int sn = strlen(s);

for (int i=0;i<sn;i++) {
char c = s[i];
if (c==' ') continue;

if (strchr("abcdefgr0123456789-",c)) {
int notnote_flag = 0;
int rest_flag = 0;

double val = 0;
switch(c) {
case 'a': val = 220; break;
case 'b': val = 247; break;
case 'c': val = 262; break;
case 'd': val = 294; break;
case 'e': val = 330; break;
case 'f': val = 349; break;
case 'g': val = 392; break;
case 'r': rest_flag = 1; break;
default: notnote_flag = 1;break;
}
if (notnote_flag) {
//this is a numeric token
val = atof(&s[i]);
while (i<sn && !strchr(": )",s[i])) i++; i--; //skip to next part of token
} else {
//note token
if (i< sn-1 && strchr("#&",s[i+1])) {
if (s[i+1]=='#') { val *= pow(2,1/12.0); }
if (s[i+1]=='&') { val /= pow(2,1/12.0); }
i++;
}
if (i< sn-1 && strchr("1234567890",s[i+1])) {
octave = s[i+1] - 48;
i++;
}
while (i<sn && !strchr(": )",s[i])) i++; i--; //skip to next part of token
val *= pow(2,octave - 4);
}

double dur = 1.0;
if (i < sn-2 && s[i+1]==':') {
i+=2;
dur = atof(&s[i]);
while (i<sn && !strchr(": )",s[i])) i++; i--; //skip to next part of token
}
dur_list[ml_ptr] = dur;
val_list[ml_ptr++] = val;
seq_length++;
seq_dur+= dur;
if (ml_ptr >= max_notes) goto pend;
continue;
}

if (c=='(') {
//new stack frame
if (st_ptr==255) { printf("parse_str2: Stack overflow\n"); return; }
stack[st_ptr++] = seq_length; seq_length = 0;
durstack[dst_ptr++] = seq_dur; seq_dur = 0;
}

if (c==')') {
double dur = 1.0;
int mult = 1;
if (i < sn-2 && s[i+1]==':') {
i+=2;
dur = atof(&s[i]);
//advance whitespace
while (i < sn && s[i] != '*' && s[i] != ' ' && s[i] != '\n' && s[i] != ')') i++;
if (s[i]=='*' || s[i]==')') i--;
}
if (i < sn-2 && s[i+1]=='*') {
// printf("processing mult\n");
i+=2;
mult = atoi(&s[i]);
//advance whitespace
while (i < sn && s[i] != ' ' && s[i] != '\n' && s[i] != ')') i++;
if (s[i]==')') i--;
}

//divide previous sequence by the sequence duration(each has equal values)
for (int j=0;j<seq_length;j++) dur_list[ml_ptr-1-j] *= dur/seq_dur;

//copy multiples if necessary
for (int k=0;k<mult-1;k++) {
for (int j=0;j<seq_length;j++) {
dur_list[ml_ptr] = dur_list[ml_ptr-seq_length];
val_list[ml_ptr] = val_list[ml_ptr-seq_length];
ml_ptr++;
if (ml_ptr >= max_notes) goto kdone;
}
}

kdone:
//pop previous sequence length from stack
if (st_ptr==0) { printf("parse_str2: Stack underflow\n"); return; }

int this_seq_length = seq_length;
seq_length = stack[--st_ptr];
seq_dur = durstack[--dst_ptr];
seq_length += (this_seq_length * mult);
seq_dur += dur*mult;
if (ml_ptr >= max_notes) goto pend;
}
}

pend:
if (0) 0;

double t =0;
for (int i=0;i<ml_ptr;i++) {
p->tvals[p->n] = t;
p->yvals[p->n] = val_list[i];
p->n++;
t += dur_list[i];
}

free(val_list);
free(dur_list);

}


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(uint64 n) {
double *a = (double *)calloc(n, sizeof(double));
if (a==NULL) { printf("CALLOC FAIL!\n");
return NULL;
}
return a;
}


//clamp high and low values
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));
}
//RANDOM

#define RAND_MAX_P 2147483647

//LCD PRNG with 16 streams
unsigned long long rs[16] = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};
unsigned long long rs_initial[16] = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};
int seed_index = -1;
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;
}

unsigned long long randp_getseed(int i) {
return rs[i];
}

int randp_getstream() { return seed_index; }
void randp_setstream(int i) {
if (i>=0 && i< 16) 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() {
if (seed_index < 0) return rand() % RAND_MAX_P;
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;
}

double choose2(double c1,double c2) { return (randp() % 2 ==0 ? c1 : c2); }

double choose3(double c1,double c2,double c3) {
int x = randp() % 3;
if (x==0) return c1;
if (x==1) return c2;
return c3;
}

//sample from exponential distribution
double exprand(double lambda) {
return -1.0 * log(drand(0,1))/lambda;
}


//sample normal distribution via central limit theorem
double nrand() {
double x = 0;
for (int i=0;i<16;i++) {
x += (double)rand()/RAND_MAX;
}
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;
}
}
}

//symmetric truncated normal distribution, bounded by n standard deviations
double tnrands(double mean,double sd, int n) {
return tnrand(mean,sd,mean-sd*n, mean+sd*n);
}



void bset_init(bset *b) {
b->interp = INTERP_LINEAR; //linear by default
b->destroy = 0;
b->n = 0;
b->max_n = 2048*16;
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);

//parameters
b->p = (double *)calloc(b->max_n*8,sizeof(double)); //param (gamma, tau);
b->pi = (int *)calloc(b->max_n,sizeof(int)); //parameter indices

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;
}

//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, exponential interpolation
double bset_interp_exp(bset *b, double t, double tau) {
int i = bset_search(b,t) - 1;

if (i >= 0 && i < b->n-1) {
//exponential interpolation
double a = (t - b->t[i])/ (b->t[i+1]-b->t[i]);
double trans_val = 0.75;
if (a >= trans_val) {
double exp_v = (b->y[i] - b->y[i+1])*exp(-1.0*a*tau) + b->y[i+1];
double lin1 = (b->y[i] - b->y[i+1])*exp(-1.0*trans_val*tau) + b->y[i+1];
double c = (a-trans_val)/(1-trans_val);
double lin_v = (1-c)*lin1 + c*b->y[i+1];
return (1-c)*exp_v + c*lin_v;
} else {
return (b->y[i] - b->y[i+1])*exp(-1.0*a*tau) + b->y[i+1];
}
}
return 0;
}

//interpolate breakpoint set at time t, power law
double bset_interp_gamma(bset *b, double t,double gamma) {
int i = bset_search(b,t) - 1;

if (i >= 0 && i < b->n-1) {
//power law interpolation
double a = (t - b->t[i])/ (b->t[i+1]-b->t[i]);
a = pow(a,1.0/gamma);
return (1-a)*b->y[i] + a*b->y[i+1];
}
return 0;
}

//interpolate breakpoint set at time t, sample/hold
double bset_interp_cubic(bset *b, double t) {
int i = bset_search(b,t) - 1;
if (i >= 0 && i < b->n-1) {
//cubic interpolation
double a = -2*(b->y[i+1] - b->y[i]);
double c = 3 * (b->y[i+1] - b->y[i]);

double d = (t - b->t[i])/(b->t[i+1] - b->t[i]);
return a*d*d*d + c*d*d + b->y[i];
}
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;
case INTERP_CUBIC:
return bset_interp_cubic(b,t);
break;
case INTERP_GAMMA:
return bset_interp_gamma(b,t,b->iparam[i]); //default gamma param 5.0
break;
case INTERP_EXP:
return bset_interp_exp(b,t,b->iparam[i]); //default exp param 6.0
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) {
double t = 0;
int k=0;
while (t < b->t[b->n-1]) {
data[k++] = bset_interp(b,t);
t+=1.0/b->sr;
}
}

void bset_dump(bset *b) {
for (int i=0;i<b->n;i++) {
printf("bset[%f] = %f\n",b->t[i],b->y[i]);
}
}

void bset_write(bset *b,char *filename) {
FILE *f = fopen(filename,"w");
for (int i=0;i<b->n;i++) {
fprintf(f,"%lf\t%lf\n",b->t[i],b->y[i]);
}
fclose(f);
}

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--;
}

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; }

//make a copy and append, with optional time gap between
void bset_dup(bset *b, double times, double gap) {
double dur = b->t[b->n-1];

int n = b->n;
for (int i=0;i<n;i++) {
for (int j=0;j<(times-1);j++) {
bset_add(b,b->t[i] + (dur + gap)*(j+1), b->y[i]);
}
}
}

void bset_parse(bset *b, char *s) {
parse_t p;
parse_init(&p);
parse_str(&p,s);
for (int i=0;i<p.n;i++)
bset_add(b,p.tvals[i],p.yvals[i]);
}

//uses alternate syntax for nested rhythms
void bset_parse2(bset *b, char *s) {
parse_t p;
parse_init(&p);
parse_str2(&p,s);
bset_clear(b);

for (int i=0;i<p.n;i++)
bset_add(b,p.tvals[i],p.yvals[i]);
}




void bset_stretch(bset *b,double f) {
for (int i=0;i<b->n;i++) {
b->t[i] *= f;
}
}


void bset_shift(bset *b,double f) {
//traverse backwards so points can be safely removed
for (int i=b->n-1;i>=0;i--) {
b->t[i] -= f;
if (b->t[i] < 0) bset_remove(b,i);
}
}

void bset_delay(bset *b,double f) {
for (int i=b->n-1;i>=0;i--) {
b->t[i] += f;
}
}


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_place(bset *b, bset *b2, double t) {
for (int i=0;i<b2->n;i++) {
bset_add(b,b2->t[i]+t,b2->y[i]);
}
}

void bset_splice(bset *b, double t1, double t2) {
if (t2 <= t1) {
printf("bset_splice: invalid time interval\n");
return;
}

double y1 = bset_interp(b,t1);
double y2 = bset_interp(b,t2);
// printf("t1: %.3f y1: %.3f, t2: %.3f y2: %.3f\n",t1,y1,t2,y2);
double si = bset_search(b,t1);

bset c;
bset_init(&c);

for (int i=0;i<b->n;i++) {
double t = b->t[i];
if (t >=t1 && t <= t2) continue;
bset_add(&c,b->t[i],b->y[i]);
}

// add splice points

bset_add(&c,t1-(1e-9),y1);
bset_add(&c,t1,0);
bset_add(&c,t2-(1e-9),0);
bset_add(&c,t2,y2);

bset_clear(b);
for (int i =0;i<c.n;i++) {
b->t[i] = c.t[i];
b->y[i] = c.y[i];
}
b->n = c.n;
bset_destroy(&c);
}


void bset_erase(bset *b, double t1, double t2) {
if (t2 <= t1) {
printf("bset_erase: invalid time interval\n");
return;
}

bset c;
bset_init(&c);

for (int i=0;i<b->n;i++) {
double t = b->t[i];
// if (t > t1 && t < t2) continue;
bset_add(&c,b->t[i],b->y[i]);
}

bset_copy(&c,b);
// bset_copy_fast(&c,b);
bset_destroy(&c);
}


//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;
}
}


/* apply and amplitude envelope from a bset */



void bset_scale(bset *b, double s, double t) {
for (int i=0;i<b->n;i++) {
b->y[i] *= s;
b->y[i] += t;
}
}

bset *bp(char *s) {
bset *b = (bset *)calloc(1, sizeof(bset));
bset_init(b);
b->destroy = 1;
bset_parse2(b,s);
return b;
}

void bset_replacem(bset *b, char *p) {
char buf[4096];
strcpy(buf,p);
// printf("%s\n",buf);
char *s = strtok(buf," ");
char *s2 = strtok(NULL, " ");
bset d;
bset_init(&d);

while (1) {
if (s == NULL || s2 == NULL) break;
if (s2 == NULL) break;

double t = atof(s);
double v = atof(s2);
bset_add(&d,t,v);
// printf(" adding %.3f %.3f\n",t,v);

s = strtok(NULL," ");
s2 = strtok(NULL, " ");
}
// printf("yep");
// bset_dump(&d);

double t1 = d.t[0];
double t2 = d.t[d.n-1];

double y1 = bset_interp_linear(b,t1);
double y2 = bset_interp_linear(b,t2);
double si = bset_search(b,t1);

bset c;
bset_init(&c);

for (int i=0;i<b->n;i++) {
double t = b->t[i];
if (t >=t1 && t <= t2) continue;
bset_add(&c,b->t[i],b->y[i]);
}

// add splice points
bset_add(&c,t1-1e-9,y1);

bset_place(&c,&d,0);

bset_add(&c,t2+1e-9,y2);
bset_clear(b);
for (int i =0;i<c.n;i++) {
b->t[i] = c.t[i];
b->y[i] = c.y[i];
}
b->n = c.n;

bset_destroy(&c);
bset_destroy(&d);

}

double dmod(double x, double y) {
return x - (int)(x/y) * y;
}



// 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]; }
}

// 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]; }
}


// 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; }
}

//reverse
void reverse(double *data,int n) {
for (int i=0;i<n/2;i++) {
double t = data[i];
data[i] = data[(n-1)-i];
data[(n-1)-i] = t;
}
}



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

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;
}


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


//FILTERS

// 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;
}

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; }
}



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;
}
}

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 minsample(double *data,int n) {
double m = 0;
for (int i=0;i<n;i++) {
double x = data[i];
if (x < m) m = x;
}
return m;
}

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




// 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);
}


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

double dur_func(double t, double d1,double d2) {
double d = (d2-d1)*t + d1;
return nrand() * d/3 + d;
}

double os_rate = 1;

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



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;
}


double tex_lookup(double *src,double xxx,double yyy, int w, int h) {
//lookup texture using bilinear interpolation

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;

}



double l1norm(double *val,int n) {
double t = 0;
for (int i=0;i<n;i++) {
t += fabs(val[i]);
}
return t;
}


void setnorm(double *val,int n) {
double p1 = l1norm(val,n)/(double)n;
printf("norm is: %.4f\n",p1);

if (p1 > 1e-5) {
scale(1.0/p1 * 0.5,0,val,n);
}
}


//offset sin, 0.0 to 1.0, with frequency f
double sino(double f, double t) {
return 0.5+0.5*sin(2*M_PI*f*t);
}

//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);
int imx = mx;
int imy = my;

for (int i=0;i<w;i++) {
a = 0;
for (int j=0;j<h;j++) {
int xx = i;
int yy= j-my;

// if (imx%2==0) yy+=1;

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;

// if (imy%2==0) xx+=1;

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);
}

//breakpoint set
typedef struct list {
int n; // number of samples used
int w;
int h;
int top; //start of queue
int bot; //end of queue
int *lx;
int *ly;
int *mask;
} list;



void write_frame(int w,int h,double *val,int *fn,double *t) {
char filename[128];
// if (*fn+1 == 34)

sprintf(filename,"frames/frame%04d.bmp",*fn+1);
writebmp2(val,filename,w,h);

if (*fn % 10==0) {
printf("Frame %d\n",*fn+1);
}

*t += 1/24.0;
*fn = *fn + 1;
}



double edgefunc(double *a,double *b,double *c) {
return (c[0] - a[0]) * (b[1] - a[1]) - (c[1]-a[1])*(b[0] - a[0]);
}

int drawtri_mode;

void drawtri(double *_p,double col, double *v, int w, int h) {
double p[6];
for (int i=0;i<6;i++) p[i] = _p[i];

double min_x,min_y,max_x,max_y;
min_x = w; max_x = 0;
min_y = h; max_y = 0;

//check crossproduct
/*
double v1x = p[2] - p[0];
double v1y = p[3] - p[1];
double v2x = p[4] - p[0];
double v2y = p[5] - p[1];
*/
// double cp = v1x*v2y - v1y*v2x;
double cp = (p[2]-p[0])*(p[5]-p[1]) - (p[3] - p[1])*(p[4] - p[0]);

if (cp > 0) {
swap(&p[4],&p[2]);
swap(&p[5],&p[3]);
}

for (int i=0;i<3;i++) {
min_x = p[i*2] < min_x ? p[i*2] : min_x;
min_y = p[i*2+1] < min_y ? p[i*2+1] : min_y;
max_x = p[i*2] > max_x ? p[i*2] : max_x;
max_y = p[i*2+1] > max_y ? p[i*2+1] : max_y;
}
for (int x = min_x; x< max_x; x++) {
for (int y = min_y; y< max_y; y++) {
double j[2];
j[0] = x+0.5; j[1] = y+0.5;
double w0 = edgefunc(&p[2],&p[4],j);
double w1 = edgefunc(&p[4],&p[0],j);
double w2 = edgefunc(&p[0],&p[2],j);
if (w0 >= 0 && w1 >= 0 && w2 >=0) {
if (x >= 0 && x < w && y >= 0 && y < h) {

if (drawtri_mode==0)
v[ind(x,y,w)] = col;
else
v[ind(x,y,w)] += col;

}
}
}
}
}

void drawquad(double *p,double col, double *v, int w, int h) {
double o[8];
for (int i=0;i<8;i++) o[i] = p[i];
o[2] = o[0]; o[3] = o[1];
drawtri(p,col,v,w,h);
drawtri(&p[2],col,v,w,h);
drawtri(&o[2],col,v,w,h);

}

typedef struct {
int st;
int lev;
double x0,y0,x1,y1;
} sub_t;


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



double dist(double x1,double y1,double x2,double y2) {return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));}

double mdist(int t,double x1,double y1,double x2,double y2) {

if (t==1) return fabs(x1-x2) + fabs(y1-y2); //manhattan
if (t==2) return pow(pow(x1-x2,8) + pow(y1-y2,8),1.0/8.0); //approx chebychev?
//euclidian by default
return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
}


void trsc(double x, double y, double sx,double sy, double *v, int n) {
for (int i=0;i<n;i++) {
v[i*2] =v[i*2]*sx+x;
v[i*2+1] =v[i*2+1]*sy+y;
}
}

void square3(double siz, double dx, double dy,double a,
double t, double *v, int w, int h) {
static double off = 0;
siz *= a;

double x1 = 0,y1=0;
double x = x1 - siz*w;
double y = y1 - siz*h;

double x2 = x1 + siz*w;
double y2 = y1 + drand(-1.0,1.0)*siz*h;

off += drand(-1,1)*4e-3*w;

//off is a static variable, cumulative offset.
// it wanders, and then is reset every period
// it affects the randomization of some of the parts
if (sin((1.0/3.0)*t) > 0.98) off = 0;

x += dx*w; y += dy*h;
x2 += dx*w; y2 += dy*h;
double p[8];

p[0] = x; p[1] = y;
p[2] = x2; p[3] = y;
p[4] = x2; p[5] = y2;
p[6] = x; p[7] = y2;

p[randp()%8]+= off;
p[randp()%8]+= off;
p[randp()%8]+= off;

if (randp()%2==0) {
drawquad(p,drand(0,1.0),v,w,h);
} else {
drawtri(p,drand(0,1.0),v,w,h);
}
}

void square4(double siz, double dx, double dy,double a,
double t, double *v, int w, int h) {
static double off = 0;
siz *= a;

double x1 = 0,y1=0;
double x = x1 - siz*w;
double y = y1 - siz*h;

double x2 = x1 + siz*w;
// double y2 = y1 + drand(-1.0,1.0)*siz*h;
double y2 = y1 + siz*h;

off += drand(-1,1)*4e-3*w;

//off is a static variable, cumulative offset.
// it wanders, and then is reset every period
// it affects the randomization of some of the parts
if (sin((1.0/3.0)*t) > 0.98) off = 0;

x += dx*w; y += dy*h;
x2 += dx*w; y2 += dy*h;
double p[8];

p[0] = x; p[1] = y;
p[2] = x2; p[3] = y;
p[4] = x2; p[5] = y2;
p[6] = x; p[7] = y2;

p[randp()%8]+= off;
p[randp()%8]+= off;
p[randp()%8]+= off;

if (randp()%10==0) {
drawquad(p,drand(0,1.0),v,w,h);
} else {
drawtri(p,drand(0,1.0),v,w,h);
}
}


void square5(double siz, double dx, double dy,double a,double off,
double t, double *v, int w, int h,double col) {

siz *= a;

double x1 = 0,y1=0;
double x = x1 - siz*w;
double y = y1 - siz*h;

double x2 = x1 + siz*w;
double y2 = y1 + drand(-1.0,1.0)*siz*h;

x += dx*w; y += dy*h;
x2 += dx*w; y2 += dy*h;
double p[8];

p[0] = x; p[1] = y;
p[2] = x2; p[3] = y;
p[4] = x2; p[5] = y2;
p[6] = x; p[7] = y2;

p[randp()%8]+= off;
p[randp()%8]+= off;
p[randp()%8]+= off;

if (randp()%2==0) {
drawquad(p,drand(0,1.0)*col,v,w,h);
} else {
drawtri(p,drand(0,1.0)*col,v,w,h);
}
}

void square6(double siz, double dx, double dy,double dx2,double dy2,double a,double off,
double t, double *v, int w, int h,double col) {

siz *= a;

double x1 = 0,y1=0;
double x = x1 - siz*w;
double y = y1 - siz*h;

double x2 = x1 + siz*w;
double y2 = y1 + drand(-1.0,1.0)*siz*h;

x += dx*w; y += dy*h;
x2 += dx2*w; y2 += dy2*h;
double p[8];

p[0] = x; p[1] = y;
p[2] = x2; p[3] = y;
p[4] = x2; p[5] = y2;
p[6] = x; p[7] = y2;

p[randp()%8]+= off;
p[randp()%8]+= off;
p[randp()%8]+= off;

if (randp()%2==0) {
drawquad(p,drand(0,1.0)*col,v,w,h);
} else {
drawtri(p,drand(0,1.0)*col,v,w,h);
}
}


// for pixels over a threshold, 50% chance of reducing their value
// proportionally by factor a
void noise1(double thresh, double a, double *v,int w,int h) {
int n = w*h;
double *v2 = make_bufn(n);

for (int i=0;i<n;i++) {
if (v[i] >thresh && randp()%2==0) { v2[i] = v[i]*a; }
else { v2[i] = v[i]; }
}
lerp(1.0,v2,v,n);
free(v2);
}


// for pixels over a threshold, 50% chance of reducing their value
// proportionally by factor a
void noise2(double th1, double th2, double a, double *v,int w,int h) {
int n = w*h;
double *v2 = make_bufn(n);

for (int i=0;i<n;i++) {
if (randp()%2==0) {v2[i] = v[i]; continue; }
if (v[i] >th2) { v2[i] = v[i]*a; }
else if (v[i] > th1) {
double a2 = (v[i]-th1)/(th2-th1);
v2[i] = v[i] * (a2 * a + (1-a2) * 1.0);
}
else { v2[i] = v[i]; }
}
lerp(1.0,v2,v,n);
free(v2);
}



void makefloor(double *v, int w, int h) {
double *g = make_bufn(w*h);
for (int x=0;x<w;x++) {
for (int y=0.9*h;y<h;y++) {
v[ind(x,y,w)] = 0.5;
}
}

val_boxblur(w*0.01,h*0.01,v,w,h);
val_boxblur(w*0.01,h*0.01,v,w,h);

noisefield(g,w,h);
scale(0.5,0.5,g,w*h);
mult(g,v,w*h);
free(g);
}

//bresenham
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)
val[ind(x,y,w)] += v;

if (k >= 4096) { //printf("LINE overflow: %d\n",k);
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;
}
}
}




void background(double *v, double sp, double off,
double off2,
double _am,
double ra, int w, int h) {

drawtri_mode = 1;
double p[8];

double th = w*(1.0/800);
if (th > 1) {
// printf("Using drawtri for background\n");
} else {
// printf("Using val_line for background\n");
}

for (double x =off-off2; x < off-off2 + 1.0;x+=sp) {
double xx = dmod(x,1.0);
if (xx < 0) xx += 1.0;
double am = _am * (1+drand(-1,1)*ra);


p[0] = xx*w - 0;
p[1] = 0;
p[2] = xx*w - 0;
p[3] = h;

p[4] = xx*w + th;
p[5] = h;

p[6] = xx*w + th;
p[7] = 0;

if (th >= 1) {
drawquad(p,am,v,w,h);
} else {
val_line(xx*w,0,xx*w,h,am,v,w,h);
}
}

for (double x =off; x < 1.2;x+=sp) {
double am = _am * (1+drand(-1,1)*ra);
double y = x;

p[0] = 0;
p[1] = y*h - 0;
p[2] = h;
p[3] = y*h - 0;

p[4] = h;
p[5] = y*h + th;

p[6] = 0;
p[7] = y*h + th;

if (th >= 1) {
drawquad(p,am,v,w,h);
} else {
val_line(0,y*h,w,y*h,am,v,w,h);
}

}
drawtri_mode = 0;
}

void addthresh(double th1, double th2, double *src, double *dst,int n) {
for (int i=0;i<n;i++) {
if (dst[i] > th2) continue;
if (dst[i] > th1) {
double a2 = (dst[i]-th1)/(th2-th1);
dst[i] += src[i] * a2;
} else {
dst[i] += src[i];
}
}
}

void noisefield(double *v, int w, int h) {
int n = w*h;
double wx[32];
double wy[32];
double px[32];
double py[32];
for (int i=0;i<32;i++) {
wx[i] = drand(0,1) * exp(-i*0.3);
wy[i] = drand(0,1) * exp(-i*0.3);
px[i] = drand(0,M_PI*2);
py[i] = drand(0,M_PI*2);
}
fill(0,v,n);

for (int i=0;i<n;i++) {
int x,y;
dni(i,&x,&y,w);

double xx = (double)x/w;
double yy = (double)y/h;

for (int j=0;j<3;j++) {

v[i] += wx[j]*cos(0.2*j*M_PI*2*xx + px[j]*M_PI*2);
v[i] += wy[j]*cos(0.2*j*M_PI*2*yy + py[j]*M_PI*2);
}
}
double x = minsample(v,n);
for (int i=0;i<n;i++) {
v[i] -=x;
}
normalize(v,n);
/*
for (int j=0;j<1;j++) {
val_boxblur(w*0.01,h*0.01,v,w,h);
}
*/

}


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) {

int wrap = 1;
for (int x2=dx;x2<dx+dw;x2++) {
for (int y2=dy;y2<dy+dh;y2++) {
int x = x2;
int y = y2;

while (x < 0) x += w2;
x %=w2;
while (y < 0) y += h2;
y %=h2;

if (x < 0 || x >= w2 || y < 0 || y >= h2) continue;

//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;

//wrap texture
if (wrap) {
while (nnx < 0) nnx += 1.0;
while (nny < 0) nny += 1.0;
nnx = dmod(nnx,1.0);
nny = dmod(nny,1.0);
}
double av;
if (am >=0) {
av = am*tex_lookup(v1,nnx,nny,w,h);
} else {
av = am * (tex_lookup(v1,nnx,nny,w,h) -1.0);
}

// printf("%d %d -> %.2f %.2f: %.3f\n",x,y,nx/w,ny/h,a);
if (mode==1) {
v2[ind(x,y,w2)] = av;
} else {
v2[ind(x,y,w2)] += av;
}
}
}
}


// four legged creature
void anim0(double tt, double *v, int w, int h) {
static int init=0;
static int num = 6;
static double f1[16],f2[16];
static double r1[16],r2[16];
static double ph1[16],ph2[16];
static double ph1b[16],ph2b[16];
static double px[16],py[16];
static double _px[16],_py[16];
static double psiz[16];
static double la1[16],la2[16];
static double fa1[16],fa2[16];

if (!init) {
init = 1;
for (int i=0;i<16;i++) {
px[i] = drand(0.3,0.4);
py[i] = drand(0.3,0.6);

if (i > 1) px[i] = drand(0.6,0.7);

_px[i] = px[i];
_py[i] = py[i];

r1[i] = drand(0.2,0.4);
r2[i] = drand(0.1,(1-py[i]) - r1[i]);

f1[i] = drand(0.7,1.1);

// f2[i] = drand(f1[i],f1[i]*1.3);
f2[i] = f1[i];
ph1[i] = drand(0,M_PI*2);
ph2[i] = drand(0,M_PI*2);

ph1b[i] = drand(0,M_PI*2);
ph2b[i] = drand(0,M_PI*2);
psiz[i] = drand(0.04,0.15);
la1[i] = drand(0.03,0.11);
la2[i] = drand(0.03,0.12);

fa1[i] = drand(0.1,0.3); //thickness of legs
fa2[i] = drand(0.1,0.3);
}
}


for (int i=0;i<4;i++) {
px[i] = _px[i] + sin(M_PI*2*f1[i]*tt + ph1[i])*la1[i];
py[i] = _py[i] + sin(M_PI*2*f2[i]*tt + ph2[i])*la2[i];

if (i==2) {

square5(psiz[i/2] * (1 + drand(-1,1)*0.05),
px[i],py[i],1.0, 0.1*w,
tt,v,w,h,0.3);

}
// ph1 = 0.1*sin(3*M_PI*2*tt);
// ph2 = 0.1*sin(4*M_PI*2*tt);
double ang1 = tt*M_PI*2*f1[i] + ph1[i];
double ang2 = -tt*M_PI*2*f2[i] + ph2[i];

ang1 = M_PI*0.3*cos(M_PI*2*f1[i]*tt + ph1b[i]) + M_PI*0.5;

ang2 = M_PI*0.3*cos(M_PI*2*f2[i]*tt + ph2b[i]) + M_PI*0.5 + ang1;
if (i < 2)
ang2 = -(M_PI*0.3*cos(M_PI*2*f2[i]*tt + ph2b[i]) + M_PI*0.3) + ang1;


double u[8],p[8];
double x = px[i] + r1[i] * cos(ang1);
double y = py[i] + r1[i] * sin(ang1);
p[0] = px[i];
p[1] = py[i];
p[2] = x;
p[3] = y;
p[4] = (px[i]+x)/2.0 + fa1[i]*(p[1] - p[3]);
p[5] = (py[i]+y)/2.0 + fa1[i]*(p[2] - p[0]);

for (int i=0;i<8;i++)
u[i] = w*p[i];
drawtri(u,drand(0.0,1.0)*0.3,v,w,h);

double x2 = x + r2[i]*cos(ang2);
double y2 = y + r2[i]*sin(ang2);
p[0] = x; p[1] = y;
p[2] = x2; p[3] = y2;
p[4] = (x+x2)*0.5 - fa2[i]*(p[1] - p[3]);
p[5] = (y+y2)*0.5 + fa2[i]*(p[2] - p[0]);

for (int i=0;i<8;i++)
u[i] = w*p[i];
drawtri(u,drand(0.0,1.0)*0.3,v,w,h);

square5(0.02 + drand(0,0.02),
x2,y2,1.0, 0.03*w,
tt,v,w,h,0.3);

}

//body
double u[8],p[8];
p[0] = px[0];
p[1] = py[0];
p[2] = px[1];
p[3] = py[1];
p[4] = px[2];
p[5] = py[2];
p[6] = px[3];
p[7] = py[3];
for (int i=0;i<8;i++) u[i] = w*p[i];

double off = 3e-2*w;

u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;

if (randp()%2==0)
drawtri(u,drand(0,1.0)*0.3,v,w,h);
else
drawquad(u,drand(0,1.0)*0.3,v,w,h);

}


// a bunch of shapes go from floor to in line vertically
void anim1(double tt, double *v, int w, int h) {
static int init=0;
static int num = 6;
static double siz;
if (!init) {
init = 1;
num = 6;
siz = 0.012;
}
for (int i=0;i<num;i++) {
double x = i/(double)num + 0.5/(double)num;

double a = (1+cos(0.6*M_PI*2*tt))*0.5;
double ap5 = pow(a,1);
double a2 = pow(a,2);

double dx = 0.5*(1-a2) + x*a2;
double dy = 0.8*ap5 + 0.1*(1-ap5);
double c = (1+sin(3*tt + i*M_PI*2.0/num))*0.5;

square4(siz*(c*7+1),dx,dy+0.1*c,1.0,tt,v,w,h);
}
}


// a shifting polygonal figure, vertices move
// in oscillating motion
void anim2(double tt, double *v, int w, int h) {
static int init=0;
static int num = 6;
static double siz;
static double *ph;
static double *amp;
static double *p,*f,*off;

if (!init) {
init = 1;
num = 20;
ph = make_bufn(num);
p = make_bufn(num);
amp = make_bufn(num);
f = make_bufn(num);
off = make_bufn(num);
for (int i=0;i<num;i++) {
off[i] = drand(0.3,0.7);
ph[i] = drand(0,2*M_PI);
amp[i] = drand(0.1,0.2);
if (randp()%8==0) amp[i] = drand(0.1,0.5);
f[i] = drand(0.2,0.7);
}
}

double u[8];
for (int i=0;i<num;i++) {
p[i] = off[i] + amp[i] * cos(tt*f[i] * M_PI*2 + ph[i]);
if (i%2==1) {
p[i] += 0.3*cos(0.5*tt*M_PI*2);
}
}

for (int i=0;i<num-8;i+=2) {
for (int j=0;j<8;j++) {
u[j] = w*p[i+j];
}

if (randp()%2==0)
drawtri(u,drand(0.0,1.0)*0.3,v,w,h);
else
drawquad(u,drand(0.0,1.0)*0.3,v,w,h);
}
}

// human running
void anim11(double tt, double *v, int w, int h) {
static int init=0;
static int num = 6;
static double f1[16],f2[16]; //frequency
static double r1[16],r2[16];
static double ph1[16],ph2[16];
static double ph1b[16],ph2b[16];
static double px[16],py[16];
static double _px[16],_py[16];
static double psiz[16];
static double la1[16],la2[16];
static double fa1[16],fa2[16];
static double hip_amp[16],knee_amp[16];

if (!init) {
init = 1;
for (int i=0;i<16;i++) {
px[i] = drand(0.48,0.52);
py[i] = drand(0.49,0.57);

if (i > 1) py[i] = drand(0.22,0.27);

_px[i] = px[i];
_py[i] = py[i];

r1[i] = drand(0.15,0.33);
r2[i] = drand(0.15,(0.5-r1[i]));
if (py[i] + r1[i] + r2[i] > 0.97) {
printf("TOO LONG\n");
double da = (py[i] + r1[i] + r2[i]) - 0.97;
r1[i] -= da/2;
r2[i] -= da/2;
}
f1[i] = drand(0.6,0.7); //frequency and phases of joints
ph1[i] = drand(0,M_PI*2);
ph2[i] = drand(0,M_PI*2);

ph1b[i] = drand(0,M_PI*2);
ph2b[i] = drand(0,M_PI*2);

if (i==1) {
//right leg

}

// f2[i] = drand(f1[i],f1[i]*1.3);
// frequencies are synced

if (i==1) {
//right leg, sync phases 180
// ph1[i] = ph1[0];
ph1b[i] = ph1b[0] + M_PI;

//knee joint
ph2b[i] = ph2b[0] + M_PI;

f1[i] = f1[0]; //sync freqs
}

if (i==3) {
//right arm in phase with legs
ph1b[2] = ph1b[0] + M_PI; //OPTION

//right arm, sync phases 180
ph1b[i] = ph1b[2] + M_PI;

f1[i] = f1[0]; //sync freqs
}

f2[i] = f1[i]; // freqs are same for shoulder/elbow, hip/knee

psiz[i] = drand(0.04,0.08);
la1[i] = drand(0.01,0.032); //how much it bounces
la2[i] = drand(0.03,0.032);

fa1[i] = drand(0.1,0.3); //thickness of legs
fa2[i] = drand(0.1,0.3);

if (i <= 1) {
//legs
hip_amp[i] = 0.15;
knee_amp[i] = 0.8;
} else {
//arms
hip_amp[i] = 0.3; //shoulder
knee_amp[i] = 0.8; //elbows
}

}
}

for (int i=0;i<4;i++) {
px[i] = _px[i] + sin(M_PI*2*f1[i]*tt + ph1[i])*la1[i];
py[i] = _py[i] + sin(M_PI*2*f2[i]*tt + ph2[i])*la2[i]*0.01;

if (i==2) {
square5(psiz[i/2] * (1 + drand(-1,1)*0.05),
px[i],py[i]-0.1 + pow(sin(M_PI*2*f2[i]*tt + ph2[i]),2)*0.04,
1.0, 0.1*w,
tt,v,w,h,0.3);

}
// ph1 = 0.1*sin(3*M_PI*2*tt);
// ph2 = 0.1*sin(4*M_PI*2*tt);
double ang1 = tt*M_PI*2*f1[i] + ph1[i];
double ang2 = -tt*M_PI*2*f2[i] + ph2[i];

ang1 = M_PI*hip_amp[i]*cos(M_PI*2*f1[i]*tt + ph1b[i]);

if (i <= 1) {
//legs
// ang1 += 0.8*sin(0.3*M_PI*2*tt + ph2b[i]);
}
ang1 = limit(-M_PI*0.6,M_PI*0.6,ang1);
ang1 += M_PI*0.5; //angles of hips, shoulders

ang2 = M_PI*knee_amp[i]*cos(M_PI*2*f2[i]*tt + ph2b[i]);
ang2 = limit(-M_PI*0.05,M_PI*0.6,ang2);
if (i > 1) {
//arms
ang2 = -ang2 - M_PI*0.15 + ang1; //angle of elbows is inverse
} else {
//legs
// ang1 += -1.3*(0.5 + 0.5*sin(0.3*M_PI*2*tt + ph2b[i]));


ang2 = ang2 + M_PI*0.15 + ang1; //angle of knees
}

double u[8],p[8];
double x = px[i] + r1[i] * cos(ang1);
double y = py[i] + r1[i] * sin(ang1);
p[0] = px[i];
p[1] = py[i];
p[2] = x;
p[3] = y;
p[4] = (px[i]+x)/2.0 + fa1[i]*(p[1] - p[3]);
p[5] = (py[i]+y)/2.0 + fa1[i]*(p[2] - p[0]);

for (int i=0;i<8;i++)
u[i] = w*p[i];
drawtri(u,drand(0.0,1.0)*0.3,v,w,h);

double x2 = x + r2[i]*cos(ang2);
double y2 = y + r2[i]*sin(ang2);
p[0] = x; p[1] = y;
p[2] = x2; p[3] = y2;
p[4] = (x+x2)*0.5 - fa2[i]*(p[1] - p[3]);
p[5] = (y+y2)*0.5 + fa2[i]*(p[2] - p[0]);

for (int i=0;i<8;i++)
u[i] = w*p[i];
drawtri(u,drand(0.0,1.0)*0.3,v,w,h);

square5(0.02 + drand(0,0.02),
x2,y2,1.0, 0.03*w,
tt,v,w,h,0.3);

}

//body
double u[8],p[8];
p[0] = px[0];
p[1] = py[0];
p[2] = px[1];
p[3] = py[1];
p[4] = px[2];
p[5] = py[2];
p[6] = px[3];
p[7] = py[3];
for (int i=0;i<8;i++) u[i] = w*p[i];

double off = 7e-2*w;

u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;

if (randp()%2==0)
drawtri(u,drand(0,1.0)*0.3,v,w,h);
else
drawquad(u,drand(0,1.0)*0.3,v,w,h);

}


// big dog
void anim12(double tt, double *v, int w, int h) {
static int init=0;
static int num = 6;
static double f1[16],f2[16]; //frequency
static double r1[16],r2[16];
static double ph1[16],ph2[16];
static double ph1b[16],ph2b[16];
static double px[16],py[16];
static double _px[16],_py[16];
static double psiz[16];
static double la1[16],la2[16];
static double fa1[16],fa2[16];
static double hip_amp[16],knee_amp[16];

if (!init) {
init = 1;
for (int i=0;i<16;i++) {
px[i] = 0.35;
py[i] = 0.51;

if (i > 1) px[i] = 0.75;

_px[i] = px[i];
_py[i] = py[i];

r1[i] = 0.26;
r2[i] = 0.21;

f1[i] = drand(0.6,0.7); //frequency and phases of joints
ph1[i] = drand(0,M_PI*2);
ph2[i] = drand(0,M_PI*2);

ph1b[i] = drand(0,M_PI*2);
ph2b[i] = drand(0,M_PI*2);

if (i==1) {
//right leg
}

// f2[i] = drand(f1[i],f1[i]*1.3);
// frequencies are synced

if (i==1) {
//right leg, sync phases 180
// ph1[i] = ph1[0];
ph1b[i] = ph1b[0] + M_PI + drand(-1,1)*0.05;

//knee joint
ph2b[i] = ph2b[0] + M_PI + drand(-1,1)*0.05;

f1[i] = f1[0]; //sync freqs
}

if (i==3) {
//right arm in phase with legs
ph1b[2] = ph1b[0] + M_PI*0.5 + drand(-1,1)*0.05; //OPTION

//right arm, sync phases 180
ph1b[i] = ph1b[2] + M_PI + drand(-1,1)*0.05;

f1[i] = f1[0]; //sync freqs
}

f2[i] = f1[i]; // freqs are same for shoulder/elbow, hip/knee

psiz[i] = drand(0.04,0.08);
la1[i] = drand(0.01,0.032); //how much it bounces
la2[i] = drand(0.03,0.032);

// thickness of legs here
fa1[i] = 0.035; //drand(0.01,0.03); //thickness of legs
fa2[i] = 0.25*fa1[i];

if (i <= 1) {
//back legs
hip_amp[i] = 0.12;
knee_amp[i] = 0.1;
} else {
//front legs (arms)
hip_amp[i] = 0.25; //shoulder
knee_amp[i] = 0.1; //elbows
}

}
}

for (int i=0;i<4;i++) {
/*
px[i] = _px[i] + sin(M_PI*2*f1[i]*tt + ph1[i])*la1[i];
py[i] = _py[i] + sin(M_PI*2*f2[i]*tt + ph2[i])*la2[i]*0.01;
*/

if (0 && i==2) {
square5(psiz[i/2] * (1 + drand(-1,1)*0.05),
px[i],py[i]-0.1 + pow(sin(M_PI*2*f2[i]*tt + ph2[i]),2)*0.04,
1.0, 0.1*w,
tt,v,w,h,0.3);

}
// ph1 = 0.1*sin(3*M_PI*2*tt);
// ph2 = 0.1*sin(4*M_PI*2*tt);
double ang1 = tt*M_PI*2*f1[i] + ph1[i];
double ang2 = -tt*M_PI*2*f2[i] + ph2[i];

ang1 = cos(M_PI*2*f1[i]*tt + ph1b[i]);
// ang1 = ang1*ang1*ang1;
ang1 *= M_PI*hip_amp[i];


if (i <= 1) {
//legs
//ang1 += 0.8*sin(0.3*M_PI*2*tt + ph2b[i]);
}

// ang1 = limit(-0.1*M_PI,M_PI*0.4,ang1);

// ang1 = limit(-0.4*M_PI,M_PI*0.8,ang1);

ang1 += M_PI*0.8; //angles of hips, shoulders

ang2 = M_PI*knee_amp[i]*cos(M_PI*2*f2[i]*tt + ph2b[i]);
// ang2 = limit(-M_PI*0.05,M_PI*0.6,ang2);

ang2 = -ang2 - M_PI*0.52 + ang1; //angle of knees inverse too

double u[8],p[8];
double x = px[i] + r1[i] * cos(ang1);
double y = py[i] + r1[i] * sin(ang1);
double nx = y - py[i];
double ny = px[i] - x;
double l = sqrt(nx*nx + ny*ny);
nx /= l; ny /= l;
double vx = ny;
double vy = -nx;

double ovs = 0.06;
p[0] = px[i] + vx*ovs + nx*fa1[i];
p[1] = py[i] + vy*ovs + ny*fa1[i];

p[2] = px[i] + vx*ovs - nx*fa1[i];
p[3] = py[i] + vy*ovs - ny*fa1[i];

p[4] = x - nx*fa1[i]*0.5;
p[5] = y - ny*fa1[i]*0.5;

p[6] = x + nx*fa1[i]*0.5;
p[7] = y + ny*fa1[i]*0.5;

for (int i=0;i<8;i++)
u[i] = w*p[i];
drawquad(u,drand(0.0,1.0)*0.3,v,w,h);

double x2 = x + r2[i]*cos(ang2);
double y2 = y + r2[i]*sin(ang2);


ny = x-x2;
nx = y2 - y;
l = sqrt(nx*nx + ny*ny);
nx /= l; ny /= l;

p[0] = x + nx*fa2[i];
p[1] = y + ny*fa2[i];

p[2] = x - nx*fa2[i];
p[3] = y - ny*fa2[i];

p[4] = x2 - nx*fa2[i]*0.5;
p[5] = y2 - ny*fa2[i]*0.5;

p[6] = x2 + nx*fa2[i]*0.5;
p[7] = y2 + ny*fa2[i]*0.5;


for (int i=0;i<8;i++)
u[i] = w*p[i];
drawquad(u,drand(0.0,1.0)*0.3,v,w,h);

double ang3 = M_PI*0.5;
double x3 = x2 + 0.05*cos(ang3);
double y3 = y2 + 0.05*sin(ang3);

square5(0.016 + drand(0,0.01),
x3,y3,1.0, 0.03*w,
tt,v,w,h,0.3);

}

//body
double u[8],p[8];
double ox = 0.04;

p[0] = px[0] - 3*ox;
p[1] = py[0] - ox;
p[2] = px[1] - 2.5*ox;
p[3] = py[1] + ox;

p[4] = px[2] + 1.8*ox;
p[5] = py[2] + 1.2*ox;
p[6] = px[3] + 3*ox;
p[7] = py[3] - 1.5*ox;

for (int i=0;i<8;i++) u[i] = w*p[i];

double off = 1e-2*w;

u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;

if (randp()%2==0)
drawtri(u,drand(0,1.0)*0.3,v,w,h);
else
drawquad(u,drand(0,1.0)*0.3,v,w,h);

}

typedef struct {
double p[8];
int sh; //shape, 0 is tri, 1 is quad
} member_t;


//ridgid linkage
typedef struct {
member_t m[16]; //up to 16 members
double dx[16]; //translation
double dy[16];
double tx[16]; // translated point of rotation for next member
double ty[16];
double ra[16]; // rotation angle
double maxa[16]; //maximum angle
double mina[16]; //minimum angle

double ph[3]; //amplitude and phase of each
double amp[3];

} linkage_t;


double rotx(double ra,double x, double y) {
return cos(ra)*x - sin(ra)*y;

}

double roty(double ra, double x, double y) {
return sin(ra)*x + cos(ra)*y;
}

void draw_linkage(linkage_t *l, int num, double ov, double *v, int w, int h) {
double ra=0;
double tx=0,ty=0;
for (int i=0;i<num;i++) {
tx += l->dx[i]; ty += l->dy[i];
// rotation, clamped by mina, maxa

double a = limit(-ov, 1 + ov,l->ra[i]);

ra += (1-a)*l->mina[i] + a*l->maxa[i];

//apply translation
double u[8];
double *p = l->m[i].p;

//apply translation from p to u
for (int j=0;j<4;j+=1) {
u[j*2] = cos(ra)*p[j*2] - sin(ra)*p[j*2+1] + tx;
u[j*2+1] = sin(ra)*p[j*2] + cos(ra)*p[j*2+1] + ty;
}
for (int j=0;j<8;j++) u[j] *= w;


double off = 3e-2*w;

u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;

if (randp()%2==0)
drawtri(u,drand(0,1.0)*0.3,v,w,h);
else
drawquad(u,drand(0,1.0)*0.3,v,w,h);

// modify to translate tx[i],ty[i], rotate by ra[i],
tx += rotx(ra,l->tx[i],l->ty[i]);
ty += roty(ra,l->tx[i],l->ty[i]);
}
}

void deg2rad(double *a, int n) {
for (int i=0;i<n;i++) {
a[i] = (a[i]/360)*M_PI*2;
}
}

void init_hindleg(linkage_t *l) {
double le[3] = {0.2,0.15,0.15};
double we[3] = {0.05,0.02,0.01};

double mina[3] = {-70,10,-50};
double of[3] = {-30,60,-30};
double maxa[3] = {20,120,0};

deg2rad(mina,3);
deg2rad(of,3);
deg2rad(maxa,3);

double ph[3] = {0.3,0.0,0.4};
double amp[3] = {0.5,0.5,1.0};


for (int i=0;i<3;i++) {
l->ph[i]= ph[i];
l->amp[i] = amp[i];

double *p = l->m[i].p;
p[0] = -we[i]; p[1] = 0;
p[2] = we[i]; p[3] = 0;
p[4] = we[i]; p[5] = le[i];
p[6] = -we[i],p[7] = le[i];

l->tx[i] = 0;
l->ty[i] = le[i];

l->ra[i] = (of[i] - mina[i])/(maxa[i] - mina[i]);

l->mina[i] = mina[i];
l->maxa[i] = maxa[i];
}
}

void init_frontleg(linkage_t *l) {
double le[3] = {0.18,0.13,0.13};
double we[3] = {0.023,0.02,0.01};

double mina[3] = {-30,-100,10};
double of[3] = {-30,60,-30};
double maxa[3] = {70,-20,90};

deg2rad(mina,3);
deg2rad(of,3);
deg2rad(maxa,3);

double ph[3] = {0.0,0.25,-0.05};
double amp[3] = {0.5,0.5,0.9};

for (int i=0;i<3;i++) {
l->ph[i]= ph[i];
l->amp[i] = amp[i];

double *p = l->m[i].p;
p[0] = -we[i]; p[1] = 0;
p[2] = we[i]; p[3] = 0;
p[4] = we[i]; p[5] = le[i];
p[6] = -we[i],p[7] = le[i];

l->tx[i] = 0;
l->ty[i] = le[i];

l->ra[i] = (of[i] - mina[i])/(maxa[i] - mina[i]);

l->mina[i] = mina[i];
l->maxa[i] = maxa[i];
}
}

//gallop
void anim3(double tt, double *v, int w, int h) {
static int init = 0;
static linkage_t l[4];
static double px[8], py[8];
static double _px[8],_py[8];
static double ph[4] = {0,0.2,0.4,0.6};
static double f = 0.5;

px[0] = 0.3;
py[0] = 0.5;
px[1] = 0.3;
py[1] = 0.5;

px[2] = 0.65;
py[2] = 0.5;
px[3] = 0.65;
py[3] = 0.5;

if (init==0) {
init = 1;

f = drand(0.55,0.65);

//initialize linkage
for (int j=0;j<2;j++) {

init_hindleg(&l[j]);
l[j].dx[0] = px[j];
l[j].dy[0] = py[j];
}

for (int j=2;j<4;j++) {
init_frontleg(&l[j]);
l[j].dx[0] = px[j];
l[j].dy[0] = py[j];
}

double _ph = drand(0,1);
for (int j=0;j<4;j++) ph[j] += _ph;
}



double oamp = tnrand(1,0.25,0.5,1.5);
double ov = tnrand(0,0.1,-0.2,0.2);

double phl[4] = {0.1,0.1,0.6,0.6};
double la1[4] = {0.0,0.,0.0,0.0};
double la2[4] = {0.015,0.015,0.03,0.03};

//update linkage angles
for (int j=0;j<4;j++) { //each linkage

l[j].dx[0] = px[j] + sin(M_PI*2*f*tt + phl[j]*M_PI*2)*la1[j];
l[j].dy[0] = py[j] + sin(M_PI*2*f*tt + phl[j]*M_PI*2)*la2[j];

for (int i=0;i<3;i++) {
l[j].ra[i] = 0.5 + oamp*l[j].amp[i]*sin(M_PI*2*f*tt + (l[j].ph[i] + ph[j])*M_PI*2);
}

draw_linkage(&l[j],3,ov,v,w,h);
}


//body
double u[8],p[8];
double ox = 0.04;


p[0] = l[0].dx[0] - 1*ox;
p[1] = l[0].dy[0] - ox;
p[2] = l[1].dx[0] - 0.8*ox;
p[3] = l[1].dy[0] + 1.5*ox;

p[4] = l[2].dx[0] + 1*ox;
p[5] = l[2].dy[0] + 3.0*ox;
p[6] = l[3].dx[0] + 2*ox;
p[7] = l[3].dy[0] - 1.5*ox;

for (int i=0;i<8;i++) u[i] = w*p[i];

double off = 3e-2*w;

u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;

if (randp()%2==0)
drawtri(u,drand(0,1.0)*0.3,v,w,h);
else
drawquad(u,drand(0,1.0)*0.3,v,w,h);


//shapes

double s[2][8] ={ {0.5,0.5,0.7,0.5,0.8,0.3},
{0.8,0.3,0.9,0.4,0.7,0.4}};

s[0][0] = (l[0].dx[0] + l[2].dx[0])*0.5;
s[0][1] = (l[0].dy[0] + l[2].dy[0])*0.5;

s[0][2] = p[4];
s[0][3] = p[5];


{
for (int j=0;j<2;j++) {
for (int i=0;i<8;i++) u[i] = w*s[j][i];

double off = 3e-2*w;

u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;

drawtri(u,drand(0,1.0)*0.3,v,w,h);
}
}

square5(0.055 * (1 + drand(-1,1)*0.05),
0.8,0.3,1.0, 0.1*w,
tt,v,w,h,0.3);

square5(0.04 * (1 + drand(-1,1)*0.05),
p[0]-0.05,p[1],1.0, 0.03*w,
tt,v,w,h,0.3);

}

// human walking
void anim10(double tt, double *v, int w, int h) {
static int init=0;
static int num = 6;
static double f1[16],f2[16]; //frequency
static double r1[16],r2[16];
static double ph1[16],ph2[16];
static double ph1b[16],ph2b[16];
static double px[16],py[16];
static double _px[16],_py[16];
static double psiz[16];
static double la1[16],la2[16];
static double fa1[16],fa2[16];

if (!init) {
init = 1;
for (int i=0;i<16;i++) {
px[i] = drand(0.48,0.52);
py[i] = drand(0.49,0.57);

if (i > 1) py[i] = drand(0.22,0.27);

_px[i] = px[i];
_py[i] = py[i];

r1[i] = drand(0.15,0.33);
r2[i] = drand(0.15,(0.5-r1[i]));
if (py[i] + r1[i] + r2[i] > 0.97) {
printf("TOO LONG\n");
double da = (py[i] + r1[i] + r2[i]) - 0.97;
r1[i] -= da/2;
r2[i] -= da/2;
}
f1[i] = drand(0.6,0.7);

if (i==1) f1[i] = f1[0];

// f2[i] = drand(f1[i],f1[i]*1.3);
f2[i] = f1[i];
ph1[i] = drand(0,M_PI*2);
ph2[i] = drand(0,M_PI*2);

ph1b[i] = drand(0,M_PI*2);
ph2b[i] = drand(0,M_PI*2);

if (i==1) {
// ph1[i] = ph1[0];
ph1b[i] = ph1b[0] + M_PI;
}
psiz[i] = drand(0.04,0.08);
la1[i] = drand(0.01,0.032); //how much it bounces
la2[i] = drand(0.03,0.032);

fa1[i] = drand(0.1,0.3); //thickness of legs
fa2[i] = drand(0.1,0.3);
}
}


for (int i=0;i<4;i++) {
px[i] = _px[i] + sin(M_PI*2*f1[i]*tt + ph1[i])*la1[i];
py[i] = _py[i] + sin(M_PI*2*f2[i]*tt + ph2[i])*la2[i]*0.01;

if (i==2) {
square5(psiz[i/2] * (1 + drand(-1,1)*0.05),
px[i],py[i]-0.1,1.0, 0.1*w,
tt,v,w,h,0.3);

}
// ph1 = 0.1*sin(3*M_PI*2*tt);
// ph2 = 0.1*sin(4*M_PI*2*tt);
double ang1 = tt*M_PI*2*f1[i] + ph1[i];
double ang2 = -tt*M_PI*2*f2[i] + ph2[i];

ang1 = M_PI*0.15*cos(M_PI*2*f1[i]*tt + ph1b[i]) + M_PI*0.5;

ang2 = M_PI*0.15*cos(M_PI*2*f2[i]*tt + ph2b[i]) + M_PI*0.15 + ang1;
if (i > 1)
ang2 = -(M_PI*0.15*cos(M_PI*2*f2[i]*tt + ph2b[i]) + M_PI*0.15) + ang1;

double u[8],p[8];
double x = px[i] + r1[i] * cos(ang1);
double y = py[i] + r1[i] * sin(ang1);
p[0] = px[i];
p[1] = py[i];
p[2] = x;
p[3] = y;
p[4] = (px[i]+x)/2.0 + fa1[i]*(p[1] - p[3]);
p[5] = (py[i]+y)/2.0 + fa1[i]*(p[2] - p[0]);

for (int i=0;i<8;i++)
u[i] = w*p[i];
drawtri(u,drand(0.0,1.0)*0.3,v,w,h);

double x2 = x + r2[i]*cos(ang2);
double y2 = y + r2[i]*sin(ang2);
p[0] = x; p[1] = y;
p[2] = x2; p[3] = y2;
p[4] = (x+x2)*0.5 - fa2[i]*(p[1] - p[3]);
p[5] = (y+y2)*0.5 + fa2[i]*(p[2] - p[0]);

for (int i=0;i<8;i++)
u[i] = w*p[i];
drawtri(u,drand(0.0,1.0)*0.3,v,w,h);

square5(0.02 + drand(0,0.02),
x2,y2,1.0, 0.03*w,
tt,v,w,h,0.3);

}

//body
double u[8],p[8];
p[0] = px[0];
p[1] = py[0];
p[2] = px[1];
p[3] = py[1];
p[4] = px[2];
p[5] = py[2];
p[6] = px[3];
p[7] = py[3];
for (int i=0;i<8;i++) u[i] = w*p[i];

double off = 7e-2*w;

u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;
u[randp()%8]+= drand(-1,1)*off;

if (randp()%2==0)
drawtri(u,drand(0,1.0)*0.3,v,w,h);
else
drawquad(u,drand(0,1.0)*0.3,v,w,h);

}

typedef struct {
int init;
int num;
double x;
double a1,a2;
double f1,f2;
double r1,r2;
double _ph1,_ph2,_ph3;
} anim4_struct_t;


void anim4(double tt, double *v, int w, int h, anim4_struct_t *d) {
if (!d->init) {
d->init = 1;
d->r1 = 0.3;
d->r2 = 0.2;
d->f1 = 0.4;
d->f2 = 0.7;
d->_ph1 = drand(0,M_PI*2);
d->_ph2 = drand(0,M_PI*2);
d->_ph3 = drand(0,M_PI*2);
}

double ph1 = 0.1*sin(3*M_PI*2*tt);
double ph2 = 0.1*sin(4*M_PI*2*tt);

double u[8],p[8];
p[0] = d->x;
p[1] = 0.1 +0.2 +0.3*(cos(tt*M_PI*2*0.3 + d->_ph3));
p[2] = 0.5+d->r1*cos(tt*M_PI*2*d->f1 + ph1 + d->_ph1);
p[3] = 0.5+0.1*sin(tt*M_PI*3*d->f1 + ph2 + d->_ph2);
p[4] = d->x;
p[5] = 0.9;

int ne = tnrand(1,1,1,5);

if (ne < 0) ne = 1;

for (int i=0;i<8;i++)
u[i] = w*p[i];

for (int j=0;j<ne;j++) {

drawtri(u,drand(0.0,1.0)*0.3,v,w,h);
double off = 3e-2*w;

u[randp()%8]+= off;
u[randp()%8]+= off;
u[randp()%8]+= off;
}

square5(0.05 + ((sin(3*tt*M_PI*2)+1)*0.5*.1)
,p[2],p[3],1.0,1e-2*w,
tt,v,w,h,1.0);

}


//circle
void anim5(double tt, double *v, int w, int h) {
static int init=0;
static int num;
static double a1,a2;
static double f1,f2;
static double r1,r2;
static double ph1,ph2;

if (!init) {
num = randp()%4 + 4;

init = 1;
r1 = r2 = drand(0.2,0.43);

if (randp()%4==0) {
r2 = drand(0.2,0.43);
}

f1 = 0.4;
f2 = 0.7;
}

for (int i=0;i<num;i++) {
double ang = (double)i/num*M_PI*2;

ang += 0.3*tt*M_PI*2; //OPTIONAL could vary these speeds

double x = 0.5 + r1 * cos(ang);
double y = 0.5 + r2 * sin(ang);

double c = (1+sin(1.5*M_PI*2*tt + i*M_PI*2.0/num))*0.5;
c = pow(c,3);
double siz = 0.02;
square5(siz*(c*22*r1+1),
x,y,1.0,0.05*w,
tt,v,w,h,1.0);
}

}


void anim6(double tt, double *v, int w, int h) {
static int init=0;
static int num = 6;
static double a1,a2;
static double f1,f2;
static double r1,r2;
static double ph1,ph2;
static double ph3;
static bset *xb,*yb;

if (!init) {
init = 1;
xb = bp("((0.8 0.1):0.3 (0.2 0.7):0.2):1.2*20 ");
yb = bp("((0.9 0.9):0.3 (0.9 0.75 0.68 0.75 0.9):0.2):1.2*20 ");

r1 = 0.3;
r2 = 0.2;

f1 = 0.4;
f2 = 0.7;
ph3 = drand(0,1.2);
}

double u[8],p[8];
p[0] = 0.55;
p[1] = 0.1 + 0.6*(cos(tt*M_PI*2*0.3)+1)*0.5;
p[2] = 0.55+r1*cos(tt*M_PI*2*f1 + ph1);
p[3] = 0.55+0.1*sin(tt*M_PI*3*f1 + ph1);
p[4] = bset_interp(xb,tt + ph3);
p[5] = bset_interp(yb,tt + ph3);

for (int i=0;i<8;i++)
u[i] = w*p[i];
drawtri(u,drand(0.0,1.0)*0.3,v,w,h);

square5(0.02 + ((sin(3*tt*M_PI*2)+1)*0.5*.1)
,p[4],p[5],1.0,1e-2*w,
tt,v,w,h,1.0);


p[0] = 0.45;
p[1] = 0.1 + 0.3*(cos(tt*M_PI*2*0.3)+1)*0.5;
p[2] = 0.45+r1*cos(tt*M_PI*2*f1 + ph1);
p[3] = 0.45+0.1*sin(tt*M_PI*3*f1 + ph1);

p[4] = bset_interp(xb,tt+0.55 + ph3);
p[5] = bset_interp(yb,tt+0.55 + ph3);

for (int i=0;i<8;i++)
u[i] = w*p[i];
drawtri(u,drand(0.0,1.0)*0.3,v,w,h);

square5(0.02 + ((sin(3*tt*M_PI*2)+1)*0.5*.1)
,p[4],p[5],1.0,1e-2*w,
tt,v,w,h,1.0);

}

//walking figure
void anim7(double tt, double *v, int w, int h) {
static int init=0;
static int num = 6;
static double siz;
static double *ph;
static double *amp;
static double *p,*f,*off;
static bset *xb,*yb;
static double ph2;

if (!init) {
init = 1;
xb = bp("((0.8 0.1):0.3 (0.2 0.7):0.2):1.2*20 ");
yb = bp("((0.9 0.9):0.3 (0.9 0.75 0.68 0.75 0.9):0.2):1.2*20 ");
num = 30;
ph = make_bufn(num);
p = make_bufn(num);
amp = make_bufn(num);
f = make_bufn(num);
off = make_bufn(num);
for (int i=0;i<num;i++) {
off[i] = drand(0.2,0.6);
ph[i] = drand(0,2*M_PI);
amp[i] = drand(0.1,0.2);
if (randp()%8==0) amp[i] = drand(0.1,0.5);
f[i] = drand(0.2,0.7);
}
double ph2 = drand(0,1.2);
}

double u[8];
for (int i=0;i<num;i++) {
p[i] = off[i] + amp[i] * cos(tt*f[i] * M_PI*2 + ph[i]);

if (i%2==1) {
p[i] += 0.1*cos(0.5*tt*M_PI*2);
}
}

p[0] = bset_interp(xb,tt+ph2);
p[1] = bset_interp(yb,tt+ph2);

p[8] = bset_interp(xb,tt+ph2+0.55);
p[9] = bset_interp(yb,tt+ph2+0.55);

for (int i=0;i<num-8;i+=4) {
for (int j=0;j<8;j++) {
u[j] = w*p[i+j];
}

double off = 3e-2*w;

u[randp()%8]+= off;
u[randp()%8]+= off;
u[randp()%8]+= off;

double col = choose2(drand(0,0.2), drand(0.8,1.0));
col = drand(0,1.0) * 0.3;

if (randp()%2==0)
drawtri(u,col,v,w,h);
else
drawquad(u,col,v,w,h);
}
}


void trans(double *u, double *p, int w, double dx, double dy,double sx,double sy) {
for (int i=0;i<6;i+=2) {
u[i] = w*((p[i]-0.5)*sx+0.5 + dx);
u[i+1] = w*((p[i+1]-0.5)*sy+0.5 + dy);
}
}


// bird
void anim9(double tt, double *v, int w, int h) {
static int init=0;
static int num = 6;
static double siz;
static double *ph;
static double *amp;
static double *p,*f,*off;
static bset *xb,*yb;
static double wf, wp,wa,wu;
static double sx,sy;
static double bulge_amount;
if (!init) {
init = 1;
xb = bp("((0.8 0.1):0.3 (0.2 0.7):0.2):1.2*20 ");
yb = bp("((0.9 0.9):0.3 (0.9 0.75 0.68 0.75 0.9):0.2):1.2*20 ");
num = 30;
ph = make_bufn(num);
p = make_bufn(num);
amp = make_bufn(num);
f = make_bufn(num);
off = make_bufn(num);
for (int i=0;i<num;i++) {
off[i] = drand(0.2,0.6);
ph[i] = drand(0,2*M_PI);
amp[i] = drand(0.1,0.2);
if (randp()%8==0) amp[i] = drand(0.1,0.5);
f[i] = drand(0.2,0.7);
f[0] = 0.5;
}
wf = drand(0.1,0.6); //freq of overall sinusoidal movement
wp = drand(0,M_PI*2); //phase of overall sinusoidal movement
wa = drand(0.0,0.35); //amp of overall sinusoidal movement;

wu = choose2(-1,1);
sx = drand(1.0,1.2);
sy = sx;
bulge_amount = tnrand(0.03,0.03,0.01,0.09);
}

double dx = 0;

// double dy = wu*wa*(0.5 - 1/(1 + exp(-wf*100*(tt) + 1)));

double dy = wu*wa*exp(-(tt)*wf*6);
dy += wa*0.15 * sin(wf*M_PI*2*tt + wp);

double u[8];

double a; // = (1+sin(tt*M_PI*2*f[0] + ph[0]))*0.5;
a = pow(sin(tt*M_PI*2*f[0] + ph[0]),2);

//body
p[0] = 0.2+(1-a)*0.05; p[1] = 0.51+(1-a)*0.01;
p[2] = 0.62+(1-a)*0.05; p[3] = 0.49-(1-a)*0.01;
p[4] = 0.56+(1-a)*0.025; p[5] = 0.44;

trans(u,p,w,dx,dy,sx,sy); //transform to world coordinates

double of = 2e-2*w;

u[randp()%8]+= drand(-1,1)*of;
u[randp()%8]+= drand(-1,1)*of;

drawtri(u,drand(0.2,0.8)*0.3,v,w,h);

double px[] = {0.47,0.69,0.55,0.62};
double py[] = {0.1,0.82,0.35,0.6};

double dta = 0.5*(0.5*cos(wf*2.33*tt*M_PI*2 + wf*7) + 0.5);

a = a * (1.0-dta) + 0.5*dta;

double x = px[0]*a + px[1]*(1-a);
double y = py[0]*a + py[1]*(1-a);

//a = (1+sin(tt*M_PI*2*f[0]))*0.5;
// a = pow(a,1.5);

a = pow(sin(tt*M_PI*2*f[0] + ph[0]),2);
a = a * (1.0-dta) + 0.5*dta;

double x2 = px[2]*a + px[3]*(1-a);
double y2 = py[2]*a + py[3]*(1-a);


p[0] = x2; p[1] = y2;
p[4] = 0.53; p[5] = 0.49-0.035;
p[2] = 0.35; p[3] = 0.51-0.02;



trans(u,p,w,dx,dy,sx,sy); //transform to world coordinates
drawtri(u,drand(0,1.0)*0.3,v,w,h);


p[0] = x2; p[1] = y2;
p[2] = x; p[3] = y;
p[4] = 0.35; p[5] = 0.51 - 0.02;


trans(u,p,w,dx,dy,sx,sy); //transform to world coordinates
drawtri(u,drand(0,1.0)*0.3,v,w,h);

//feathers
int maxj = 20;
for (int j=0;j<maxj;j++) {
double r[6];
r[0] = x2; r[1] = y2;

double dj = 1.0/(double)maxj;
double b = (double)j*dj;


r[2] = b*p[2] + (1-b)*p[4];
r[3] = b*p[3] + (1-b)*p[5];
r[4] = r[2];
r[5] = r[3];
double ny = r[2] - r[0];
double nx = -1*(r[3] - r[1]);

//nx, ny is parallel to butt of feather
nx = p[2] - p[4];
ny = p[3] - p[5];
double l = sqrt(nx*nx + ny*ny);
nx /= l;
ny /= l;

//tx,ty is along the feather centre line
double tx = r[2] - r[0];
double ty = r[3] - r[1];
l = sqrt(tx*tx + ty*ty);
tx /= l; ty /= l;

double bu = 2*(0.5 - fabs(0.5-b));

double da = drand(-1,1)*bulge_amount;

da *= 0.2 + 0.8*bu;
dj *= 0.2;
dj += drand(-1,1)*0.04*bu;

r[2] += da*tx + dj*nx;
r[3] += da*ty + dj*ny;

r[4] += da*tx - dj*nx;
r[5] += da*ty - dj*ny;

for (int z=2;z<6;z++) //perturb
r[z] *= 1+drand(-1,1)*0.02;

trans(u,r,w,dx,dy,sx,sy); //transform to world coordinates
drawtri(u,drand(0,1.0)*0.15,v,w,h);
}

}


// clock
void anim8(double tt, double *v, int w, int h) {
static int init=0;
static int num = 6;
static double a1,a2;
static double f1,f2;
static double r1,r2;
static double ph1,ph2;

if (!init) {
r1 = 0.45;
r2 = 0.25;
init = 1;
f1 = 0.2;
f2 = f1/12.0;
ph1 = -M_PI*0.42;
ph2 = drand(0,M_PI*2);
}

//centre
square5(0.02 + ((sin(3*tt*M_PI*2)+1)*0.5*.1)
,0.5,0.5,1.0, 0.03*w,
tt,v,w,h,0.3);

double u[8],p[8];
double x = 0.5 + r1 * cos(tt*M_PI*2*f1 + ph1);
double y = 0.5 + r1 * sin(tt*M_PI*2*f1 + ph1);

double nx = -sin(tt*M_PI*2*f1+ph1);
double ny = cos(tt*M_PI*2*f1+ph1);

p[0] = 0.5 + nx*0.047;
p[1] = 0.5 + ny*0.047;
p[2] = x;
p[3] = y;
p[4] = 0.5 - nx*0.047;
p[5] = 0.5 - ny*0.047;

for (int i=0;i<8;i++)
u[i] = w*p[i];

drawtri(u,drand(0.0,1.0)*0.3,v,w,h);



//ticks
for (int j=0;j<12;j++) {
double a = j * (M_PI*2)/12.0;

double r3 = 0.5;

double x = 0.5 + r3 * cos(a);
double y = 0.5 + r3 * sin(a);

double nx = -sin(a);
double ny = cos(a);

p[0] = 0.5 + cos(a)*0.4 + nx*0.02;
p[1] = 0.5 + sin(a)*0.4 + ny*0.02;
p[2] = x;
p[3] = y;
p[4] = 0.5 + cos(a)*0.4 - nx*0.02;
p[5] = 0.5 + sin(a)*0.4 - ny*0.02;

for (int i=0;i<8;i++)
u[i] = w*p[i];


square5(0.02,0.5+cos(a)*0.4,0.5+sin(a)*0.4,
0.25, 0.03*w,
tt,v,w,h,0.2);


drawtri(u,drand(0.0,1.0)*0.3,v,w,h);
}


//hour hand
double x2 = 0.5 + r2*cos(tt*M_PI*2*f2 + ph2);
double y2 = 0.5 + r2*sin(tt*M_PI*2*f2 + ph2);

nx = -sin(tt*M_PI*2*f2+ph2);
ny = cos(tt*M_PI*2*f2+ph2);

p[0] = 0.5 + nx*0.032;
p[1] = 0.5 + ny*0.032;
p[2] = x2;
p[3] = y2;
p[4] = 0.5 - nx*0.032;
p[5] = 0.5 - ny*0.032;

for (int i=0;i<8;i++)
u[i] = w*p[i];

drawtri(u,drand(0.0,1.0)*0.3,v,w,h);
}


//different categories of animation
void do_anim(int a, double tt, double *v, int w, int h) {
static int init = 0;
static anim4_struct_t anim4_data[2];
if (!init) {
init = 1;
anim4_data[0].init = 0;
anim4_data[0].x = 0.25;
anim4_data[1].init = 0;
anim4_data[1].x = 0.75;
}

switch(a) {
case 0: anim0(tt,v,w,h); break; //four legged creature
case 1: anim1(tt,v,w,h); break; //falling to pieces
case 2: anim2(tt,v,w,h); break; //jumping abstract
case 3: anim3(tt,v,w,h); break; //gallop
case 4: anim4(tt,v,w,h,&anim4_data[0]);
anim4(tt,v,w,h,&anim4_data[1]);
break; //pair of figures
case 5: anim5(tt,v,w,h); break; //circlular zoopraxiscope
case 6: anim6(tt,v,w,h); break; //tri walk
case 7: anim7(tt,v,w,h); break; //shifting polygonal, walking
case 8: anim8(tt,v,w,h); break; //clock
case 9: anim9(tt,v,w,h); break; //bird
case 10: anim10(tt,v,w,h); break; //human walking
case 11: anim11(tt,v,w,h); break; //human running
case 12: anim12(tt,v,w,h); break; //robotic dog

default: break;
}
}

int pick_cat() {
int freqs[] = { 40, //four legged creature
20, //1 falling to pieces
100, //2 jumping abstract
40, //3 gallop
100, //4 pair of figures
80, //5 circular zoopraxiscope
80, //6 tri walk
130, //7 shifting poygonal
70, //8 clock
100, //9 bird
50, //10 human walking
50, //11 human running
25 // robotic dog walking
};
int num_cats = 13;

int lu[num_cats*100];
int k = 0;
for (int i=0;i<num_cats;i++)
for (int j=0;j<freqs[i];j++) lu[k++] = i;

return lu[randp()%k];
}

int negative = 0;



void microscale(int r, double *v, int cat,
double sp, double off,
double off2,
double _am,
int w, int h);


double rule_spacing = 0.025;
double grid_speed = 0.2; //how fast grid moves in background
double ms_amount; //microscale amount

int chrono() {
rule_spacing = 0.025;
grid_speed = 0.2;

int w2,h2;

unsigned long epoch_res[] = {35,50,75,100,125,150,200,250,300,350,
400,500,600,700,800,900,1000,1500,2250,3375,
5000,7600,10000,20000,40000,80000,160000,320000,640000,1280000};

int ep = token_params.epoch;
if (ep <30) {
w2 = epoch_res[ep];
} else {
w2 = epoch_res[29] * pow(2,ep-29);
}
h2 = w2;

if (ep < 10) {
ms_amount = 0;
} else {
double a = ((double)ep -10.0)/(22.0-10.0);
a = pow(a,0.8);
ms_amount = limit(0,0.25,0.25* a);
}

//override with command line param if supplied
if (token_params.res > 10) w2 = h2 = token_params.res;

printf("w: %d, h: %d\n",w2,h2);
printf("epoch %d res: %d ms_amount: %.3f\n",ep,w2,ms_amount);
if (ep==12) printf("end of intial mint\n");

if (ep==15) printf("after 2 weeks\n");
if (ep==20) printf("after 1 year\n");

double os=1; //oversampling rate during processing
os_rate = os;
double dos=1; //sampling rate before rendering (leave at one for downsample from 'os' value)
int w = w2*os;
int h = h2*os;

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

double st0 = 0;
double st1 = 0;
bset b; bset_init(&b);
bset_parse2(&b,"(1 0.3 0 0.7):1.3*64");

bset bs; bset_init(&bs);
bset_parse2(&bs,"(1 0.2 0.1):0.7*64");

double t = 0;
double dur = 60*10;
int fn = 0;

bset *_siz = bp("(0.05:3 0.25):1.0*30 (0:3 0.2):20 (0.05:3 0.25):12"); //siz
int trigger=1;
int fk = 0;

int cat = pick_cat();
if (token_params.cat > -1) {
//if cat is specific on command line
cat = token_params.cat;
}

negative = randp()%3 ==0 ? 1 : 0;

//PARAMS how many plates, what order
int nx = randp()%6+3;
int ny = randp()%3+2;
if (ny > nx ) ny = nx;

if (randp()%3==0 && nx <= 4) ny = nx;

int fmax = nx*ny;

ms_amount *= tnrand(1.0,0.1,0.9,1.5);
printf("ms_amount adjusted to %.3f\n",ms_amount);

printf("PARAM res: %dx%d\nPARAM cat: %d\nPARAM nx: %d\nPARAM ny: %d\nPARAM invert: %d\n",
w,h,cat,nx,ny,negative);

if (param_only) exit(0);

//PARAMS crop image depending on aspect ratio
double cropx = 0;//0.12;
double ar= (double)nx/(double)ny;
if (ar > 1.5) cropx = drand(0.0,0.1);
if (ar > 1.75) cropx = drand(0.1,0.2);
if (cat==3 || cat==0 || cat==12 || cat==13 || cat==4) cropx = 0; //no crop for certain categories

uint64 w_plate = nx*w*(1.0-cropx*2);
uint64 h_plate = ny*h;

double mx = w_plate*0.005;
double my = h_plate*0.01;

double *plate = make_bufn(w_plate*h_plate);

printf("Final plate %llu x %llu: %llu pixels\n",w_plate,h_plate,(uint64)w_plate*h_plate);

if (negative && randp()%1==0) {
fill(0.9,plate,w_plate*h_plate);
}

double skip_dur = 2.0/fmax; //spread out over 3 seconds (normally 6)
double dt = 1/24.0;
int skip_frames = skip_dur * (1.0/dt);


int add_floor = 1;

if (cat==8) {
//clock
grid_speed = 0;
dt = 1/24.0;
skip_dur = 4.7;
skip_frames = skip_dur * (1.0/dt);
}
if (cat==9) {
//bird
grid_speed = 0.3;
add_floor = 0;
}
if (cat==11) { grid_speed = 0.25; }

// include drawing of floor
makefloor(val[8],w,h);

// reset PRNG to known state here, to match with Noumenon animation
randp_seed(1,token_params.blocknumber*32 + token_params.noumenon_token_id);
for (int i=0;i<9999+1;i++) randps(1);

randp_seed(0,token_params.blocknumber*32 + token_params.noumenon_token_id);
for (int i=0;i<9999;i++) randp();
printf("PRNG seed at start: %llu\n",randp_getseed(0));

int x = 0;
int y = 0;
int frm_method=1;
while (t <= dur) {
fill(0,val[0],n);

double *v = val[1];

int nib = 3;
for (int j=0;j<nib;j++) {
double tt = t + (double)j/nib * 1/24.0;

do_anim(cat,tt,v,w,h); //params
}

scale(1.0/nib,0,v,n);

add(v,val[2],n);
double sa = 0.8;
double ba = 0.04;

if (randp()%18==0) {
//skip boxblur scaling on 1/18 of frames
} else {
scale(sa,0,val[2],n);
val_boxblur(w*ba,h*ba,val[2],w,h);
}

if (fk % 10==0) {
printf("Frame %d\n",fk+1);
}

if (++fk % skip_frames ==0) {
// we will process this frame in chronophotography sequence
} else {
t += dt;
continue; //skip this frame
}

lerp(1.0,val[2],val[3],n);
scale(0.3,0,val[3],n);
add(v,val[3],n);
normalize(val[3],n);

int ruler = 1;


if (ruler) { // OPT ruled background like in muybridge photos
fill(0,val[4],n);
background(val[4],rule_spacing,rule_spacing*0.5, t*grid_speed, 0.13,0.05,w,h);
background(val[4],rule_spacing*10,0.1+rule_spacing*0.5, t*grid_speed, 0.36,0.02,w,h);

noisefield(val[7],w,h);
scale(0.6,0.4,val[7],n);
mult(val[7],val[4],n);

add(val[4],val[3],n);

if (add_floor) { //PARAM for some traits, do not include the floor image
add(val[8],val[3],n);
}
}


softclip(tnrand(1.5,0.5,1.2,1.8),0.0,val[3],n);
noise2(0.08,0.12,0.8,val[3],w,h);

double ba2 = w*ba*0.1;
if (ba2 < 1) ba2=1;
if (w >= 1200) ba2 = 2;
if (w >= 1600) ba2 = 1;

if (ba2 >= 1) {
val_boxblur(ba2,ba2,val[3],w,h);
}

normalize(val[3],n);
if (ruler) {
addthresh(0.2,0.3,val[4],val[3],n);
}

//add smooth noise
noisefield(val[7],w,h);
scale(0.2,0,val[7],n);
add(val[7],val[3],n);

//exposure amount
double exp_amount = tnrand(0.1,0.07, 0.02, 2);
for (int i=0;i<n;i++) {
val[3][i] += exp_amount;
}

normalize(val[3],n);
int num_exp = 1;

fn++;
printf("Shutter, %d/%d\n",fn, fmax);


//add this frame to the plate
// do this in a separate PRNG context
if (frm_method==1) {
randp_setstream(1);

//add microscale to the image.
if (ms_amount > 1e-2) {
microscale(0,val[3],cat,rule_spacing,rule_spacing*0.5, t* grid_speed * (0.025/rule_spacing),
ms_amount,w,h);
}

normalize(val[3],n);

double dx = (double)x*(double)w_plate/nx;
double dy = (double)y*(double)h_plate/ny;
double dw = (double)w_plate/nx;
double dh = (double)h_plate/ny;

//noise for imperfect placement
dx += mx * (1+drand(-1,1)*0.2);
dy += my * (1+drand(-1,1)*0.01);

dw -= 2*mx;
dh -= 2*my;

double ang = 0.5*drand(-1,1)*(2*M_PI/360);

blit_rot(1,1.0,val[3],plate,cropx*w2,0,(1.0-2*cropx)*w2,h2,
dx,dy,dw,dh,
w2,h2,w_plate,h_plate,
cos(ang),-sin(ang),sin(ang),cos(ang));

x++; if (x >= nx) { x = 0; y++; }

randp_setstream(0);
} //end if (frm_method==1)

if (fn >= fmax) break;

t += dt;
}

// whether to invert the image to a negative
if (negative) {
//PARAM exposure in negative, controls how dark the darks
double d = 0.9; //0.7
scale(d,0,plate,w_plate*h_plate);
val_invert(plate,w_plate,h_plate);
}

printf("Writing to image.pgm\n");
fn = 0;
writepgm(plate,"image.pgm",w_plate,h_plate);

return 0;
}



// r is recursion level
double *ms_val[8*10];
int ms_res[8];
int ms_init[8];
int ms_fk[8];
int ms_fn[8];

void init_microscale(int r,int res) {
printf("Initializing microscale\n");
for (int i=0;i<8;i++) {
ms_res[i] = 0;
ms_fn[i] = 0;
ms_fk[i] = 0;
ms_init[i] = 0;
}
}

void init_microscale_level(int r, int res) {
if (ms_res[r] != 0) {
printf("Microscale level %d already initialized\n",r);
return;
}
ms_res[r] = res;
for (int i=0;i<10;i++) {
//memory;
ms_val[r*10+i] = make_bufn(res*res);
}
}

void getframe(int r, double *rv, double t, int cat, int res) {
if (r > 5) return; //maxium microscale level is 5
double **val = &ms_val[r*10]; //storing images

int w = ms_res[r];
int h = ms_res[r];
int w2 = w; int h2 = h;
int n = w*h;

if (!ms_init[r]) {
ms_init[r] = 1;

makefloor(val[8],w,h);

//init seeds
int j = 2+r;
randp_seed(j,token_params.blocknumber*32 + token_params.noumenon_token_id);
for (int i=0;i<9999+j;i++) randps(j);

printf("getframe init: r: %d res: %d\n",r,w);
ms_fk[r] = 0;
ms_fn[r] = 0;
}

int tmp = randp_getstream();

randp_setstream(2+r);

double skip_dur = 0.3;
int skip_frames = skip_dur * 24;

while (1) {
fill(0,val[0],n);
double *v = val[1];

int nib = 3;
for (int j=0;j<nib;j++) {
double tt = t + (double)j/nib * 1/24.0;

do_anim(cat,tt,v,w,h); //params
}

scale(1.0/nib,0,v,n);

add(v,val[2],n);
double sa = 0.8;
//fadeout 8 is a good one
double ba = 0.04;

if (randp()%18==0) {
//skip boxblur scaling on 1/18 of frames
} else {
scale(sa,0,val[2],n);
val_boxblur(w*ba,h*ba,val[2],w,h);
}

if (ms_fk[r]++ % skip_frames ==0) {
// we will process this frame in chronophotography sequence
} else {
t += 1/24.0;
continue;
}

lerp(1.0,val[2],val[3],n);
scale(0.3,0,val[3],n);
add(v,val[3],n);
normalize(val[3],n);

int ruler = 1;
int add_floor = 1;
if (ruler) { // OPT ruled background like in muybridge photos
fill(0,val[4],n);
background(val[4],rule_spacing,rule_spacing*0.5, t*0.2, 0.13,0.05,w,h);
background(val[4],rule_spacing*10,0.1+rule_spacing*0.5, t*0.2, 0.36,0.02,w,h);

noisefield(val[7],w,h);
scale(0.6,0.4,val[7],n);
mult(val[7],val[4],n);
//addthresh(0.1,0.2,val[4],val[3],n);
add(val[4],val[3],n);

if (add_floor) {//PARAM include floor
add(val[8],val[3],n);
}
}

softclip(tnrand(1.5,0.5,1.2,1.8),0.0,val[3],n);
noise2(0.08,0.12,0.8,val[3],w,h);

double ba2 = w*ba*0.1;
if (w >= 1440) ba2 = 2;

if (ba2 >= 1) {
val_boxblur(ba2,ba2,val[3],w,h);
}

normalize(val[3],n);
if (ruler) {
addthresh(0.2,0.3,val[4],val[3],n);
}

//add smooth noise
noisefield(val[7],w,h);
scale(0.2,0,val[7],n); //PARAM how much noise
add(val[7],val[3],n);


double exp_amount = tnrand(0.1,0.07, 0.02, 2);
for (int i=0;i<n;i++) {
val[3][i] += exp_amount;
}

normalize(val[3],n);

int num_exp = 1;

ms_fn[r]++;

lerp(1.0,val[3],rv,n);

randp_setstream(tmp);
return;
}
}


void microscale(int r, double *v, int cat,
double sp, double off,
double off2,
double _am,
int w, int h) {


int recurse = 1;

int w2 = w/ ((30.0*0.025)/rule_spacing);

if (w2 < 10) {
//if resolution is too small, end recursion
return;
}else if (w2 < 60) {
w2 = 35;
if (ms_amount > 0.05) {
w2 = 60;
}
recurse = 0;
}

init_microscale_level(r,w2);

printf("Microscale recursion level %d, res: %d\n",r,w2);

int h2 = w2;

double *v3 = make_bufn(w2*h2);
double t = 0;
int k = 0;
int tmp = randp_getstream();
randp_setstream(15);

for (double y =off; y < 1.0+off;y+=sp) {
double yy = dmod(y,1.0);
if (yy<0) yy+= 1.0;
if (++k % 10==0) printf("Row %d\n",k);

for (double x =off-off2; x < off-off2 + 1.0;x+=sp) {
double xx = dmod(x,1.0);
if (xx < 0) xx += 1.0;
// if (xx < 0) continue;

getframe(r,v3,t,cat,w2);

if (recurse && r <= 7)
microscale(r+1,v3,cat,sp,off,off2,_am,w2,h2);

t += 1.3/24;
if (t > 10) t -= 10;

blit_rot(0,drand(0.75,1.25)*_am,v3,v,
0,0,w2,h2,
xx*w,yy*h,sp*w,sp*h,
w2,h2,w,h,
1,0,0,1);
}
}

randp_setstream(tmp);
free(v3);
}


int motion() {
int w2,h2;

w2 = h2 = 320; //visual resolution
printf("w: %d, h: %d\n",w2,h2);

double os=1; //oversampling rate during processing
os_rate = os;
double dos=1; //sampling rate before rendering (leave at one for downsample from 'os' value)
int w = w2*os;
int h = h2*os;

int cat =pick_cat();
if (token_params.cat > -1) {
//if cat is specific on command line
cat = token_params.cat;
}

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


double st0 = 0;
double st1 = 0;

double t = 0;
double dur = 5;
int fn = 0;


makefloor(val[8],w,h);

int add_floor = 1;

if (cat==8) {
//clock
grid_speed = 0;
}
if (cat==9) {
//bird
grid_speed = 0.6;
add_floor = 0;
}
if (cat==11) { grid_speed = 0.25; }

while (t <= dur) {
fill(0,val[0],n);


double *v = val[1];

int nib = 3;
for (int j=0;j<nib;j++) {
double tt = t + (double)j/nib * 1/24.0;
do_anim(cat, tt,v,w,h); //params
}
scale(1.0/nib,0,v,n);

add(v,val[2],n);
double sa = 0.8;
double ba = 0.04;

if (randp()%18==0) {
//skip boxblur scaling on 1/18 of frames
} else {
scale(sa,0,val[2],n);
val_boxblur(w*ba,h*ba,val[2],w,h);
}

lerp(1.0,val[2],val[3],n);
scale(0.3,0,val[3],n);
add(v,val[3],n);
normalize(val[3],n);

int ruler = 1;
if (ruler) { // ruled background like in muybridge photos
fill(0,val[4],n);
background(val[4],0.025,0.025*0.5 , t*grid_speed, 0.13,0.05,w,h);
background(val[4],0.25,0.1+0.025*0.5 , t*grid_speed, 0.36,0.02,w,h);

noisefield(val[7],w,h);
scale(0.6,0.4,val[7],n);
mult(val[7],val[4],n);

add(val[4],val[3],n);


if (add_floor) {//whether to include the floor
add(val[8],val[3],n);
}
}

softclip(tnrand(1.5,0.5,1.2,1.8),0.0,val[3],n);
noise2(0.08,0.12,0.8,val[3],w,h);
val_boxblur(w*ba*0.1,h*ba*0.1,val[3],w,h);

normalize(val[3],n);
if (ruler) {
addthresh(0.2,0.3,val[4],val[3],n);
}

//add smooth noise
noisefield(val[7],w,h);
scale(0.2,0,val[7],n); //PARAM how much noise
add(val[7],val[3],n);

double exp_amount = tnrand(0.1,0.07, 0.02, 2);
for (int i=0;i<n;i++) {
val[3][i] += exp_amount;
}

normalize(val[3],n);

int num_exp = 1;
write_frame(w,h,val[3],&fn,&t);
}

return 0;
}

int main(int argc, char *argv[]) {
init_token_params();
token_params.res = 0;
token_params.cat = -1;
drawtri_mode = 0;

int do_chrono = 1;

if (argc > 1 && (strcmp(argv[1],"-c")==0 || strcmp(argv[1],"-p")==0 || strcmp(argv[1],"-m")==0)) {
//chronophotography
do_chrono = 1;

if (strcmp(argv[1],"-m")==0) do_chrono = 2; // create animation
if (strcmp(argv[1],"-p")==0) param_only = 1; // print out parameters only

if (argc > 2) token_params.noumenon_token_id = atoi(argv[2]);
if (argc > 3) token_params.chronophotograph_token_id = atoi(argv[3]);
if (argc > 4) token_params.blocknumber = atoll(argv[4]);
if (argc > 5) token_params.epoch = atoi(argv[5]);
if (argc > 6) token_params.cat = atoi(argv[6]); //override category
}

printf("Token Params\nchronophotograph_token_id: %d\nnoumenon_token_id(parent_id): %d\nblocknumber: %llu\nepoch: %d\n\n",token_params.chronophotograph_token_id,token_params.noumenon_token_id,token_params.blocknumber, token_params.epoch);

if (token_params.blocknumber==0) token_params.blocknumber=183472384;
for (int j=0;j<16;j++) {
randp_setstream(j);
randp_seed(j,token_params.blocknumber*32 + token_params.noumenon_token_id);
for (int i=0;i<9999+j;i++) randp();
}

randp_setstream(0);
printf("PRNG seed at start: %llu\n",randp_getseed(0));

if (do_chrono==1) {
printf("Creating chronographic plate...\n");
chrono();
}

if (do_chrono==2) {
printf("Creating animation...\n");
motion();
}
}