// wandersweep.c -- Take wander7 output and turn it into sound.
//                  Sample rate starts at 1KHz and sweeps upward.
// Probably run like: gzcat wander7out.gz | wandersweep >sweep.wav
// wav because lame likes wav?

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "sound.h"

#define MINFREQ (1000.0/SECOND)
#define MINWAVELEN (1.0/MINFREQ)

#define N_TIN (1.0 * (1<<27))
#define N_TOUT (1.0 * MINUTE)


double tinf_e, tinf_b;

// Return sample edges within input file so that output is an exponential
// sweep starting at 1000 input samples/sec and ending at N_TOUT output
// samples having sampled all N_TIN input samples.
double tin_f( double tout ) {
  return( tinf_e * ( exp( tinf_b * tout ) - 1 ) );
}


// Calculate the constants that tin_f() uses.
//
prep_tin_f() {
  int i;

  tinf_e = 1.0;
  // I did this math on the mistaken impression that [see notebook 2008.04.13]
  //     tin_f( tout ) = exp( tinf_a + tinf_b * tout ) - 1   instead of
  //                     tinf_e * ( exp( tinf_b * tout ) - 1 ).
  // ...which would have allowed calculating both constants directly.
  // Then ran it & saw that it didn't work.  Then I thought to try calculating
  // tinf_b based on a guess of tinf_e, then tinf_e based on tinf_b. 
  // Works well enough (13 iterations) so I kept it.
  // 
  // tin_f( N_TOUT ) = tinf_e * ( exp( tinf_b * N_TOUT ) - 1 ) = N_TIN
  //                   exp( tinf_b * N_TOUT ) = ( N_TIN / tinf_e ) + 1
  //                   tinf_b = log( ( N_TIN / tinf_e ) + 1 ) / N_TOUT
  //
  // tin_f( 1 ) = tinf_e * ( exp( tinf_b ) - 1 ) = MINFREQ
  //              tinf_e = MINFREQ / ( exp( tinf_b ) - 1 )
  //
  for( i = 1;  i < 20;  i++ ) {
    tinf_b = log( ( N_TIN / tinf_e ) + 1.0 ) / N_TOUT;
    tinf_e = MINFREQ / ( exp( tinf_b ) - 1 );
#ifdef BALLYHOO
    printf( "\ni = %d,  tinf_e = %f,  tinf_b = %g\n", i, tinf_e, tinf_b );
    printf( "tin_f( 0.0 ) = %f\n", tin_f( 0.0 ) );
    printf( "1 / tin_f( 1.0 ) = %f\n", 1 / tin_f( 1.0 ) );
    printf( "tin_f( MINUTE ) - N_TIN = %g\n", tin_f( MINUTE ) - N_TIN );
    printf( "=====================\n" );
#endif
  }
#ifdef BALLYHA
  printf( "tin_f( 0.0 ) = %f\n", tin_f( 0.0 ) );
  printf( "tin_f( 1.0 ) = %f\n", tin_f( 1.0 ) );
  printf( "tin_f( 2.0 ) = %f\n", tin_f( 2.0 ) );
  printf( "  ...\n" );
  printf( "tin_f( N_TOUT - 2 ) = %f\n", tin_f( N_TOUT - 2 ) );
  printf( "tin_f( N_TOUT - 1 ) = %f\n", tin_f( N_TOUT - 1 ) );
  printf( "tin_f( N_TOUT     ) = %f\n", tin_f( N_TOUT ) );
#endif
}


main( ) {
  double x = 0;
  double y = 0;
  double z = 0;
  double total;
  double mn = 0.0;
  double mx = 0.0;
  double mnt = 0.0;
  double mxt = 0.0;
  double d = 0.0;
  double tin, dt, tin0, tin9;  // 0 and 9 are short for start and end.
  int c;
  SoundStream *s;
  short sample;
  long tout;

  prep_tin_f();

  s = sound_open_stream( stdout, "w", "wav" );
  mn = mx = 0.0;
  tin0 = 0.0;

  c = 0;  
  tin9 = tin_f( 0.0 );
  for( tout = 0;  tout < N_TOUT  &&  c >= 0;  tout++ ) {
    total = 0.0;
    tin0 = tin9;
    tin9 = tin_f( tout + 1 );
    for( tin = tin0;  tin < tin9;  tin = floor( tin + 1 ) ) {
      if( tin == (long) tin ) {
        c = getchar();
	if( c < 0 )  break;

      }
      dt = fmin( floor(tin)+1, tin9 ) - tin;
      total += ( c == '='?  -80 : 80 ) * dt;
    }
    z = z * ( 1.0 - 1.0 / 4 ) + total / pow( (tin9-tin0), .0625 );
    y = y * ( 1.0 - 1.0 / 16 ) + z;
    x = x * ( 1.0 - 1.0 / 64 ) + y;
    if( x < mn ) {
      mn = x;
      mnt = ( tin0 + tin9 ) * .5;
    }
    if( x > mx ) {
      mx = x;
      mxt = ( tin0 + tin9 ) * .5;
    }
    sample = x;
    write_LR_short( sample, sample, s );
  }
  fprintf( stderr, "min = %f  at input sample ~%ld\n",
	   mn, (long) round( mnt ) );

  fprintf( stderr, "max = %f  at input sample ~%ld\n",
	   mx, (long) round( mxt ) );

  fprintf( stderr, "finally undersampling by %f\n", tin9 - tin0 );

  sound_close( s );
}