ese 519

Dependents:   PROJECT_3D_AUDIO

Files at this revision

API Documentation at this revision

Comitter:
niv17
Date:
Tue Apr 07 21:09:22 2015 +0000
Commit message:
sonic initial april 7

Changed in this revision

AudioObj.cpp Show annotated file Show diff for this revision Revisions of this file
Mixer3D.cpp Show annotated file Show diff for this revision Revisions of this file
fft.cpp Show annotated file Show diff for this revision Revisions of this file
fhtdit.cc Show annotated file Show diff for this revision Revisions of this file
fhtdit2.cc Show annotated file Show diff for this revision Revisions of this file
location.cpp Show annotated file Show diff for this revision Revisions of this file
mit_hrtf_lib.c Show annotated file Show diff for this revision Revisions of this file
velocity.cpp Show annotated file Show diff for this revision Revisions of this file
diff -r 000000000000 -r 2076b4d80327 AudioObj.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/AudioObj.cpp	Tue Apr 07 21:09:22 2015 +0000
@@ -0,0 +1,105 @@
+#include "AudioObj.h"
+//#include <iostream>
+
+Location AudioObj::getLocation () const {
+    return this->location;
+}
+
+void AudioObj::setLocation (const Location& loc) {
+    this->location = loc;
+}
+
+void AudioObj::setLocation (float x, float y, float z) {
+    this->location = Location(x,y,z);
+}
+
+Velocity AudioObj::getVelocity () const {
+    return this->velocity;
+}
+
+void AudioObj::setVelocity (const Velocity& vel) {
+    this->velocity = vel;
+}
+
+void AudioObj::setVelocity (float dx, float dy, float dz) {
+    this->velocity = Velocity(dx,dy,dz);
+}
+
+float AudioObj::getVolume() const {
+    return this->volume;
+}
+
+void AudioObj::setVolume(float vol) {
+    if (vol > 1 || vol < 0){
+        std::cout<<"Volume not in range (0-1)"<<endl;
+        //throw invalid_argument("Volume not in range (0-1)");
+    }
+    this->volume = vol;
+}
+
+void AudioObj::setRandomVolume() {
+    float randVol = (rand() % 100 + 1) / 100.0;
+	this->setVolume(randVol);
+}
+	
+bool AudioObj::isActive() const {
+    return this->active;
+}
+
+void AudioObj::setActive(bool active){
+    this->active = active;
+}
+
+bool AudioObj::isGpsObject() const {
+    return this->gpsObject;
+}
+
+bool AudioObj::isBackgroundObject() const {
+    return this->backgroundObject;
+}
+
+void AudioObj::setRepeat(bool rep) {
+    this->repeat = rep;
+}
+
+
+bool AudioObj::fillAudioData (Complex* target, unsigned int length) {
+    if(circularBuffer.readSizeRemaining() < length) {
+                return false;
+    }
+    circularBuffer.read(target, length);
+    return true;
+}
+
+void AudioObj::loadCircularBuffer() {
+    unsigned int length = (unsigned int)circularBuffer.writeSizeRemaining();
+    if(length>16384) { // TODO: How was this number chosen???
+        //cout<<"In write circ Buff : "<<length<<endl;
+        if(!(wavObject.loadMoreData(length, repeat))) {
+            // end of file reached, no repeat
+            this->active = false;
+            this->isCompleted = true;
+            return;
+        }
+        circularBuffer.write(wavObject.complexTempData, length);
+    }
+}
+
+void AudioObj::playOnceFromBeginning() {
+    this->repeat = false;
+    this->restart();
+}
+
+void AudioObj::restart() {
+    this->circularBuffer.clear();
+    this->wavObject.seekToBeginning();
+    this->loadCircularBuffer();
+    this->active = true;
+}
+
+void AudioObj::loadFromBeginning() {
+    this->circularBuffer.clear();
+    this->wavObject.seekToBeginning();
+    this->loadCircularBuffer();
+    this->active = false;
+}
diff -r 000000000000 -r 2076b4d80327 Mixer3D.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Mixer3D.cpp	Tue Apr 07 21:09:22 2015 +0000
@@ -0,0 +1,298 @@
+#include <math.h>
+#include "Mixer3D.h"
+//#include "World.h"
+#include "fft.h"
+#include "mit_hrtf_lib.h"
+#include <sys/time.h>
+
+Mixer3D::Mixer3D(int bufSize, int smpRate, int bitD) :
+bufferSize(bufSize), sampleRate(smpRate), bitDepth(bitD))
+{
+    cout << "BUF_SIZE: " << bufSize << endl;
+    cout << "SAMPLE_RATE: " << smpRate << endl;
+    cout << "BIT_DEPTH: " << bitD << endl;
+    inputAO = new Complex[2*bufferSize];
+    overlapInput = new Complex[2 * bufferSize];
+    fInput = new Complex[2 * bufferSize];
+    fFilter = new Complex[2 * bufferSize];
+    
+    azimuths = new int[World::MAX_OBJ];
+    elevations = new int[World::MAX_OBJ];
+    prevAzimuths = new int[World::MAX_OBJ];
+    prevElevations = new int[World::MAX_OBJ];
+
+    updateAngles(); // initialize azimuths and elevations
+
+    outputLeft = new Complex*[World::MAX_OBJ];
+    outputRight = new Complex*[World::MAX_OBJ];
+    overlapLeft = new Complex*[World::MAX_OBJ];
+    overlapRight = new Complex*[World::MAX_OBJ];
+
+    leftFilter = new short[2 * bufferSize];
+    rightFilter = new short[2 * bufferSize];
+   
+      
+}
+
+Mixer3D::~Mixer3D()
+{
+    delete [] inputAO;
+    delete [] outputLeft;
+    delete [] outputLeft;
+    delete [] leftFilter;
+    delete [] rightFilter;
+    delete [] complexLeftFilter;
+    delete [] complexRightFilter;
+    delete [] overlapLeft;
+    delete [] overlapRight;
+    delete [] overlapInput;
+    delete [] prevAzimuths;
+    delete [] prevElevations;
+    delete [] azimuths;
+    delete [] elevations;
+}
+
+void Mixer3D::updateAngles() {
+    AudioObj* iAudioObj;
+    Location iLocation;
+    for (int i = 0; i < myWorld->getNumObj(); i++) {
+        iAudioObj = myWorld->getAudioObj(i);
+        iLocation = iAudioObj->getLocation();
+        
+        if (iAudioObj->isGpsObject()) {
+            azimuths[i] = player.computeAzimuthFromGpsCoordinates(&iLocation);
+            elevations[i] = player.computeElevationFromGpsCoordinates(&iLocation);
+        } else {
+            azimuths[i] = player.computeAzimuth(&iLocation);
+            elevations[i] = player.computeElevation(&iLocation);
+        }
+    }
+}
+
+bool Mixer3D::isPowerOfTwo(int x) {
+    return !(x == 0) && !(x & (x - 1));
+}
+
+int Mixer3D::loadHRTF(int* pAzimuth, int* pElevation, unsigned int samplerate, unsigned int diffused, Complex *&leftFilterIn, Complex *&rightFilterIn)
+{
+    int size = mit_hrtf_get(pAzimuth, pElevation, samplerate, diffused, leftFilter, rightFilter);
+
+    if (size == 0) {
+        // TODO: Throw MIT HRTF filter does not exist exception
+    }
+    for (int i = 0; i < size; i++)
+    {
+        leftFilterIn[i] = (double)(leftFilter[i]);
+        rightFilterIn[i] = (double)(rightFilter[i]);
+    }
+
+    return size;
+}
+void Mixer3D::convolution(Complex *input, Complex *filter, Complex *output, long nSig, long nFil, long nFFT) {
+
+	if (input == NULL || filter == NULL) {
+		throw invalid_argument("Input and Filter must be non-NULL.");
+	}
+
+	if (!isPowerOfTwo(nFFT) || nFFT < nSig || nFFT < nFil) {
+        throw invalid_argument("NFFT must be a power of two, bigger than the signal length, and bigger than the filter length.");
+	}
+
+	// Perform FFT on both input and filter.
+    // TODO: "Streamline" CFFT class?
+	CFFT::Forward(input, fInput, (unsigned int)nFFT);
+	CFFT::Forward(filter, fFilter, (unsigned int)nFFT);
+
+	for (int i = 0; i < nFFT; i++) {
+		output[i] = fInput[i] * fFilter[i];
+    }
+	CFFT::Inverse(output, (unsigned int)nFFT);
+}
+
+void Mixer3D::stereoConvolution(Complex *input, Complex *leftFilter, Complex *rightFilter, Complex *leftOutput, Complex *rightOutput, long nSIG, long nFIL, long nFFT)
+{
+    // TODO: Modify parameter name, input is a data member
+    convolution(input, leftFilter, leftOutput, nSIG, nFIL, nFFT);
+    convolution(input, rightFilter, rightOutput, nSIG, nFIL, nFFT);
+}
+
+void Mixer3D::computeValidAudioObjects(vector<AudioObj*> &validAudioObjectsOut) {
+    AudioObj* iAudioObj;
+    float iDistance;
+    Location iLocation;
+    for (int i=0; i < myWorld->getNumObj(); i++) {
+        iAudioObj = myWorld->getAudioObj(i);
+        if (iAudioObj->isActive()) {
+            iLocation = iAudioObj->getLocation();
+            if (iAudioObj->isGpsObject()) {
+                if (player.isTrackingGps()) {
+                    iDistance = player.computeDistanceFromGpsCoordinates(&iLocation);
+                } else {
+                    // not tracking gps, so skip
+                    continue;
+                }
+            } else {
+                iDistance = player.computeDistanceFrom(&iLocation);
+            }
+            
+            if (iDistance < MAX_GPS_OBJ_DISTANCE) {
+                validAudioObjectsOut.push_back(iAudioObj);
+            }
+        }
+    }
+}
+
+void Mixer3D::performMix(short *ioDataLeft, short *ioDataRight)
+{
+    double runtime;
+    struct timeval start, finish;
+    gettimeofday(&start, NULL);
+    //Zero the output arrays.
+    for(int i = 0; i < bufferSize; i++)
+    {
+       ioDataLeft[i] = 0;
+       ioDataRight[i] = 0;
+    }
+    
+    //Iterate through all audio objects, obtain input data and calculate resulting audio data for each object.
+    AudioObj* iAudioObj;
+    Location iLocation;
+    float iVolume, iDistance, iAmplitudeFactor;
+    vector<AudioObj*> validAudioObjects;
+    
+    // get valid audio objects
+    computeValidAudioObjects(validAudioObjects);
+    int nValidAudioObjects = (int)validAudioObjects.size();
+    for (int i = 0; i < nValidAudioObjects; i++) {
+        iAudioObj = validAudioObjects[i];
+        
+        // loading in input data for the iteration accordingly
+        if (!(iAudioObj->fillAudioData(inputAO, bufferSize))) {
+            cout << "Audio Object" << i << " not loading input quickly enough" << endl;
+            continue;
+        }
+        
+        iVolume = iAudioObj->getVolume();
+        if (iVolume == 0) {
+            // skip 'muted' audio objects, but not after loading their input.
+            // that way they are still playing, but silently.
+            continue;
+        }
+        
+        // handle background objects
+        if (iAudioObj->isBackgroundObject()) {
+            // copy input directly to output buffers
+            for (int j=0; j < bufferSize; j++) {
+                // scale input from [-1, 1] to proper scale depending on bit depth
+                ioDataLeft[j] += inputAO[j].real()*iVolume*pow(2, bitDepth) / nValidAudioObjects;
+                ioDataRight[j] += inputAO[j].real()*iVolume*pow(2, bitDepth) / nValidAudioObjects;
+            }
+            continue;
+        }
+        
+        iLocation = iAudioObj->getLocation();
+        
+        // handle gps objects
+        if (iAudioObj->isGpsObject()) {
+            if (player.isTrackingGps()) {
+                iDistance = player.computeDistanceFromGpsCoordinates(&iLocation);
+                if (iDistance > MAX_GPS_OBJ_DISTANCE) {
+                    // not within range, skip this audio object
+                    continue;
+                }
+            } else {
+                // gps tracking hasn't started yet
+                continue;
+            }
+        } else {
+            iDistance = player.computeDistanceFrom(&iLocation);
+        }
+
+        if (iDistance == 0) {
+            iAmplitudeFactor = iVolume;
+        } else {
+            // ensure that the amplitude factor can never be > 1 to prevent clipping
+            if (iDistance > 1) {
+                iAmplitudeFactor = iVolume / pow(iDistance, 2);
+            } else {
+                iAmplitudeFactor = iVolume;
+            }
+        }
+        
+        // dynamically calculate the Azimuth and Elevation between every object and the player
+        updateAngles();
+        
+        for(int j = 0; j < bufferSize * 2; j++) {
+            if ( j >= bufferSize ) {
+                // zero pad
+                inputAO[j] = 0;
+            } else {
+                inputAO[j] *= iAmplitudeFactor;
+            }
+        }
+       
+        if (azimuths[i] != prevAzimuths[i] ||
+           elevations[i] != prevElevations[i]) {
+            // object location relative to player has changed, so fetch a new filter
+            filterLength = loadHRTF(&azimuths[i], &elevations[i], sampleRate, 1, complexLeftFilter[i], complexRightFilter[i]);
+            
+            // zero pad
+            for (int j = filterLength; j < 2 * bufferSize; j++) {
+                complexLeftFilter[i][j] = 0;
+                complexRightFilter[i][j] = 0;
+            }
+            
+            if (azimuths[i] < 0) {
+                   azimuths[i] = -azimuths[i];
+            }
+            
+            
+            //Since the filter changed, we perform a convolution on the old input with the new filter and pull out its tail.
+            stereoConvolution(overlapInput, complexLeftFilter[i], complexRightFilter[i], outputLeft[i], outputRight[i], bufferSize, filterLength, 2 * bufferSize);
+        
+            // update the overlap part for the next iteration
+            for (int j = 0; j < bufferSize; j++) {
+                    overlapLeft[i][j] = outputLeft[i][j + bufferSize];
+                    overlapRight[i][j] = outputRight[i][j + bufferSize];
+            }
+        }
+        
+        //Perform the convolution of the current input and current filter.
+        stereoConvolution(inputAO, complexLeftFilter[i], complexRightFilter[i], outputLeft[i], outputRight[i], bufferSize, filterLength, 2 * bufferSize);
+  
+        //Output the data to the output arrays in short integer format. In addition to pushing the output of the main
+        //convolution, we also need to add the overlapped tail of the last output and divide by 2. Finally, we need
+        //to divide by the number of audio objects to ensure no clipping.
+        for (int j = 0; j < bufferSize; j++)
+        {
+            ioDataLeft[j]  += (short)( (outputLeft[i][j].real() + overlapLeft[i][j].real()) / (2*nValidAudioObjects));
+            ioDataRight[j] += (short)( (outputRight[i][j].real() + overlapRight[i][j].real()) / (2*nValidAudioObjects));
+        }
+        
+        //Updating the overlapInput for the next iteration for the correpsonding object
+        for (int j = 0; j < bufferSize * 2; j++) {
+            if (j >= bufferSize) {
+                overlapInput[j] = 0;
+            } else {
+                overlapInput[j] = inputAO[j];
+            }
+        }
+
+        // TODO: If the filter has been changed, didn't we already do this?
+        //Updating the default overlap information for the next iteration if the filter won't be changed
+        for (int j = 0 ; j < bufferSize; j++) {
+            overlapLeft[i][j] = outputLeft[i][j + bufferSize];
+            overlapRight[i][j] = outputRight[i][j + bufferSize];
+        }
+    
+        //storing the Azimuth value in this iteration for the comparison for the next iteration so that
+        //we can know that whether the filter needs to be changed in the next iteration.
+        prevAzimuths[i] = azimuths[i];
+        prevElevations[i] = elevations[i];
+    }
+    gettimeofday(&finish, NULL);
+    runtime = finish.tv_sec - start.tv_sec;
+    runtime += (finish.tv_usec - start.tv_usec) / 1000000.0;
+    
+    cout << "performMix: " << runtime << endl;
+}
diff -r 000000000000 -r 2076b4d80327 fft.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/fft.cpp	Tue Apr 07 21:09:22 2015 +0000
@@ -0,0 +1,578 @@
+//   fft.cpp - impelementation of class
+//   of fast Fourier transform - FFT
+//
+//   The code is property of LIBROW
+//   You can use it on your own
+//   When utilizing credit LIBROW site
+
+//     modified by Da Cao, Oct.1st/2014.
+//     Changed into a FHT method to compute DFT
+
+#include <fstream>
+#include <iostream>
+#include <math.h>
+#include "fft.h"
+//#include "../include/complex.h"
+
+#include <cmath>
+
+//#include "revbinpermute.h"
+#include "../include/fhtloc2.h"
+#include "../include/fht.h"
+#include "../include/cmult.h"
+#include "../include/constants.h"
+#include "../include/fxttypes.h"
+#include "../include/complextype.h"
+#include "../include/revbinpermute.h"
+#include "../Sonic.h" // getPermutedIndex needs BUF_SIZE
+
+using namespace std;
+
+//Convolution function which returns the frequency representation of the result.
+//NFFT is the FFT size, nSIG is the size for the input, NFIL is the size of the filter.
+
+//fft_dit4_core_p1  define method
+// trig data used in method
+
+unsigned int getPermutedIndex (unsigned int num) {
+    unsigned int  NO_OF_BITS = sizeof(num) * 8;
+    unsigned int reverse_num = 0, i, temp;
+ 
+    for (i = 0; i < NO_OF_BITS; i++) {
+        temp = (num & (1 << i));
+        if(temp) {
+            reverse_num |= (1 << ((NO_OF_BITS - 1) - i));
+        }
+    }
+    return (reverse_num >> (NO_OF_BITS - (int)log2(BUF_SIZE)));
+}
+
+void
+fft8_dit_core_p1(double *fr, double *fi)
+// 8-point decimation in time FFT
+// isign = +1
+// input data must be in revbin_permuted order
+//.
+// Cf. Nussbaumer p.148f
+{
+    // INPUT_RE:
+    double t1r = fr[0] + fr[1];
+    double t2r = fr[2] + fr[3];
+    double t7r = t1r + t2r;
+    double t3r = fr[4] - fr[5];
+    double t4r = fr[4] + fr[5];
+    double t5r = fr[6] + fr[7];
+    double t8r = t4r + t5r;
+    double t6r = fr[6] - fr[7];
+    
+    double m0r = t7r + t8r;
+    double m1r = t7r - t8r;
+    double m2r = t1r - t2r;
+    double m3r = fr[0] - fr[1];
+    double m4r = M_SQRT1_2 * (t3r - t6r);
+    
+#define m5i t6r
+#define m6i t7r
+#define m7i t8r
+    m7i = M_SQRT1_2 * (t3r + t6r);
+    m5i = t5r - t4r;
+    m6i = fr[3] - fr[2];
+    
+    // INPUT_IM:
+    double t1i = fi[0] + fi[1];
+    double t2i = fi[2] + fi[3];
+    double t7i = t1i + t2i;
+    double t3i = fi[4] - fi[5];
+    double t4i = fi[4] + fi[5];
+    double t5i = fi[6] + fi[7];
+    double t8i = t4i + t5i;
+    double t6i = fi[6] - fi[7];
+    
+    double m0i = t7i + t8i;
+    double m1i = t7i - t8i;
+    double m2i = t1i - t2i;
+    double m3i = fi[0] - fi[1];
+    double m4i = M_SQRT1_2 * (t3i - t6i);
+    
+#define m5r t6i
+#define m6r t7i
+#define m7r t8i
+    m7r = M_SQRT1_2 * (t3i + t6i);
+    m5r = t4i - t5i;
+    m6r = fi[2] - fi[3];
+    
+#define s1r t1r
+#define s2r t2r
+#define s3r t3r
+#define s4r t4r
+    s1r = m3r + m4r;
+    s2r = m3r - m4r;
+    s3r = m6r + m7r;
+    s4r = m6r - m7r;
+    
+    // OUTPUT_RE:
+    fr[0] = m0r;
+    fr[7] = s1r + s3r;
+    fr[6] = m2r + m5r;
+    fr[5] = s2r - s4r;
+    fr[4] = m1r;
+    fr[3] = s2r + s4r;
+    fr[2] = m2r - m5r;
+    fr[1] = s1r - s3r;
+    
+#define s1i t1r
+#define s2i t2r
+#define s3i t3r
+#define s4i t4r
+    s1i = m3i + m4i;
+    s2i = m3i - m4i;
+    s3i = m6i - m7i;
+    s4i = m6i + m7i;
+    
+    // OUTPUT_IM:
+    fi[0] = m0i;
+    fi[7] = s1i + s3i;
+    fi[6] = m2i + m5i;
+    fi[5] = s2i - s4i;
+    fi[4] = m1i;
+    fi[3] = s2i + s4i;
+    fi[2] = m2i - m5i;
+    fi[1] = s1i - s3i;
+}
+// -------------------------
+
+#undef s1r
+#undef s2r
+#undef s3r
+#undef s4r
+
+#undef s1i
+#undef s2i
+#undef s3i
+#undef s4i
+
+#undef m5r
+#undef m6r
+#undef m7r
+
+#undef m5i
+#undef m6i
+#undef m7i
+// -------------------------
+//////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+static const ulong LX = 2;
+
+double allTheC [680] = {1, 0.980785, 0.92388, 0.83147, 0.707107, 0.55557, 0.382683, 0.19509, 1, 0.998795, 0.995185, 0.989177, 0.980785, 0.970031, 0.95694, 0.941544, 0.92388, 0.903989, 0.881921, 0.857729, 0.83147, 0.803208, 0.77301, 0.740951, 0.707107, 0.671559, 0.634393, 0.595699, 0.55557, 0.514103, 0.471397, 0.427555, 0.382683, 0.33689, 0.290285, 0.24298, 0.19509, 0.14673, 0.0980171, 0.0490677, 1, 0.999925, 0.999699, 0.999322, 0.998795, 0.998118, 0.99729, 0.996313, 0.995185, 0.993907, 0.99248, 0.990903, 0.989177, 0.987301, 0.985278, 0.983105, 0.980785, 0.978317, 0.975702, 0.97294, 0.970031, 0.966976, 0.963776, 0.960431, 0.95694, 0.953306, 0.949528, 0.945607, 0.941544, 0.937339, 0.932993, 0.928506, 0.92388, 0.919114, 0.91421, 0.909168, 0.903989, 0.898674, 0.893224, 0.88764, 0.881921, 0.87607, 0.870087, 0.863973, 0.857729, 0.851355, 0.844854, 0.838225, 0.83147, 0.824589, 0.817585, 0.810457, 0.803208, 0.795837, 0.788346, 0.780737, 0.77301, 0.765167, 0.757209, 0.749136, 0.740951, 0.732654, 0.724247, 0.715731, 0.707107, 0.698376, 0.689541, 0.680601, 0.671559, 0.662416, 0.653173, 0.643832, 0.634393, 0.624859, 0.615232, 0.605511, 0.595699, 0.585798, 0.575808, 0.565732, 0.55557, 0.545325, 0.534998, 0.52459, 0.514103, 0.503538, 0.492898, 0.482184, 0.471397, 0.460539, 0.449611, 0.438616, 0.427555, 0.41643, 0.405241, 0.393992, 0.382683, 0.371317, 0.359895, 0.348419, 0.33689, 0.32531, 0.313682, 0.302006, 0.290285, 0.27852, 0.266713, 0.254866, 0.24298, 0.231058, 0.219101, 0.207111, 0.19509, 0.18304, 0.170962, 0.158858, 0.14673, 0.134581, 0.122411, 0.110222, 0.0980171, 0.0857973, 0.0735646, 0.0613207, 0.0490677, 0.0368072, 0.0245412, 0.0122715, 1, 0.999995, 0.999981, 0.999958, 0.999925, 0.999882, 0.999831, 0.999769, 0.999699, 0.999619, 0.999529, 0.999431, 0.999322, 0.999205, 0.999078, 0.998941, 0.998795, 0.99864, 0.998476, 0.998302, 0.998118, 0.997925, 0.997723, 0.997511, 0.99729, 0.99706, 0.99682, 0.996571, 0.996313, 0.996045, 0.995767, 0.995481, 0.995185, 0.994879, 0.994565, 0.99424, 0.993907, 0.993564, 0.993212, 0.99285, 0.99248, 0.992099, 0.99171, 0.991311, 0.990903, 0.990485, 0.990058, 0.989622, 0.989177, 0.988722, 0.988258, 0.987784, 0.987301, 0.986809, 0.986308, 0.985798, 0.985278, 0.984749, 0.98421, 0.983662, 0.983105, 0.982539, 0.981964, 0.981379, 0.980785, 0.980182, 0.97957, 0.978948, 0.978317, 0.977677, 0.977028, 0.97637, 0.975702, 0.975025, 0.974339, 0.973644, 0.97294, 0.972226, 0.971504, 0.970772, 0.970031, 0.969281, 0.968522, 0.967754, 0.966976, 0.96619, 0.965394, 0.96459, 0.963776, 0.962953, 0.962121, 0.96128, 0.960431, 0.959572, 0.958703, 0.957826, 0.95694, 0.956045, 0.955141, 0.954228, 0.953306, 0.952375, 0.951435, 0.950486, 0.949528, 0.948561, 0.947586, 0.946601, 0.945607, 0.944605, 0.943593, 0.942573, 0.941544, 0.940506, 0.939459, 0.938404, 0.937339, 0.936266, 0.935184, 0.934093, 0.932993, 0.931884, 0.930767, 0.929641, 0.928506, 0.927363, 0.92621, 0.925049, 0.92388, 0.922701, 0.921514, 0.920318, 0.919114, 0.917901, 0.916679, 0.915449, 0.91421, 0.912962, 0.911706, 0.910441, 0.909168, 0.907886, 0.906596, 0.905297, 0.903989, 0.902673, 0.901349, 0.900016, 0.898674, 0.897325, 0.895966, 0.894599, 0.893224, 0.891841, 0.890449, 0.889048, 0.88764, 0.886223, 0.884797, 0.883363, 0.881921, 0.880471, 0.879012, 0.877545, 0.87607, 0.874587, 0.873095, 0.871595, 0.870087, 0.868571, 0.867046, 0.865514, 0.863973, 0.862424, 0.860867, 0.859302, 0.857729, 0.856147, 0.854558, 0.852961, 0.851355, 0.849742, 0.84812, 0.846491, 0.844854, 0.843208, 0.841555, 0.839894, 0.838225, 0.836548, 0.834863, 0.83317, 0.83147, 0.829761, 0.828045, 0.826321, 0.824589, 0.82285, 0.821103, 0.819348, 0.817585, 0.815814, 0.814036, 0.812251, 0.810457, 0.808656, 0.806848, 0.805031, 0.803208, 0.801376, 0.799537, 0.797691, 0.795837, 0.793975, 0.792107, 0.79023, 0.788346, 0.786455, 0.784557, 0.782651, 0.780737, 0.778817, 0.776888, 0.774953, 0.77301, 0.771061, 0.769103, 0.767139, 0.765167, 0.763188, 0.761202, 0.759209, 0.757209, 0.755201, 0.753187, 0.751165, 0.749136, 0.747101, 0.745058, 0.743008, 0.740951, 0.738887, 0.736817, 0.734739, 0.732654, 0.730563, 0.728464, 0.726359, 0.724247, 0.722128, 0.720003, 0.71787, 0.715731, 0.713585, 0.711432, 0.709273, 0.707107, 0.704934, 0.702755, 0.700569, 0.698376, 0.696177, 0.693971, 0.691759, 0.689541, 0.687315, 0.685084, 0.682846, 0.680601, 0.67835, 0.676093, 0.673829, 0.671559, 0.669283, 0.667, 0.664711, 0.662416, 0.660114, 0.657807, 0.655493, 0.653173, 0.650847, 0.648514, 0.646176, 0.643832, 0.641481, 0.639124, 0.636762, 0.634393, 0.632019, 0.629638, 0.627252, 0.624859, 0.622461, 0.620057, 0.617647, 0.615232, 0.61281, 0.610383, 0.60795, 0.605511, 0.603067, 0.600616, 0.598161, 0.595699, 0.593232, 0.59076, 0.588282, 0.585798, 0.583309, 0.580814, 0.578314, 0.575808, 0.573297, 0.570781, 0.568259, 0.565732, 0.563199, 0.560662, 0.558119, 0.55557, 0.553017, 0.550458, 0.547894, 0.545325, 0.542751, 0.540171, 0.537587, 0.534998, 0.532403, 0.529804, 0.527199, 0.52459, 0.521975, 0.519356, 0.516732, 0.514103, 0.511469, 0.50883, 0.506187, 0.503538, 0.500885, 0.498228, 0.495565, 0.492898, 0.490226, 0.48755, 0.484869, 0.482184, 0.479494, 0.476799, 0.4741, 0.471397, 0.468689, 0.465976, 0.46326, 0.460539, 0.457813, 0.455084, 0.45235, 0.449611, 0.446869, 0.444122, 0.441371, 0.438616, 0.435857, 0.433094, 0.430326, 0.427555, 0.42478, 0.422, 0.419217, 0.41643, 0.413638, 0.410843, 0.408044, 0.405241, 0.402435, 0.399624, 0.39681, 0.393992, 0.39117, 0.388345, 0.385516, 0.382683, 0.379847, 0.377007, 0.374164, 0.371317, 0.368467, 0.365613, 0.362756, 0.359895, 0.357031, 0.354164, 0.351293, 0.348419, 0.345541, 0.342661, 0.339777, 0.33689, 0.334, 0.331106, 0.32821, 0.32531, 0.322408, 0.319502, 0.316593, 0.313682, 0.310767, 0.30785, 0.304929, 0.302006, 0.29908, 0.296151, 0.293219, 0.290285, 0.287347, 0.284408, 0.281465, 0.27852, 0.275572, 0.272621, 0.269668, 0.266713, 0.263755, 0.260794, 0.257831, 0.254866, 0.251898, 0.248928, 0.245955, 0.24298, 0.240003, 0.237024, 0.234042, 0.231058, 0.228072, 0.225084, 0.222094, 0.219101, 0.216107, 0.21311, 0.210112, 0.207111, 0.204109, 0.201105, 0.198098, 0.19509, 0.19208, 0.189069, 0.186055, 0.18304, 0.180023, 0.177004, 0.173984, 0.170962, 0.167938, 0.164913, 0.161886, 0.158858, 0.155828, 0.152797, 0.149765, 0.14673, 0.143695, 0.140658, 0.13762, 0.134581, 0.13154, 0.128498, 0.125455, 0.122411, 0.119365, 0.116319, 0.113271, 0.110222, 0.107172, 0.104122, 0.10107, 0.0980171, 0.0949635, 0.091909, 0.0888536, 0.0857973, 0.0827403, 0.0796824, 0.0766239, 0.0735646, 0.0705046, 0.0674439, 0.0643826, 0.0613207, 0.0582583, 0.0551952, 0.0521317, 0.0490677, 0.0460032, 0.0429383, 0.0398729, 0.0368072, 0.0337412, 0.0306748, 0.0276081, 0.0245412, 0.0214741, 0.0184067, 0.0153392, 0.0122715, 0.00920375, 0.00613588, 0.00306796};
+
+double allTheS [680] = {0, 0.19509, 0.382683, 0.55557, 0.707107, 0.83147, 0.92388, 0.980785, 0, 0.0490677, 0.0980171, 0.14673, 0.19509, 0.24298, 0.290285, 0.33689, 0.382683, 0.427555, 0.471397, 0.514103, 0.55557, 0.595699, 0.634393, 0.671559, 0.707107, 0.740951, 0.77301, 0.803208, 0.83147, 0.857729, 0.881921, 0.903989, 0.92388, 0.941544, 0.95694, 0.970031, 0.980785, 0.989177, 0.995185, 0.998795, 0, 0.0122715, 0.0245412, 0.0368072, 0.0490677, 0.0613207, 0.0735646, 0.0857973, 0.0980171, 0.110222, 0.122411, 0.134581, 0.14673, 0.158858, 0.170962, 0.18304, 0.19509, 0.207111, 0.219101, 0.231058, 0.24298, 0.254866, 0.266713, 0.27852, 0.290285, 0.302006, 0.313682, 0.32531, 0.33689, 0.348419, 0.359895, 0.371317, 0.382683, 0.393992, 0.405241, 0.41643, 0.427555, 0.438616, 0.449611, 0.460539, 0.471397, 0.482184, 0.492898, 0.503538, 0.514103, 0.52459, 0.534998, 0.545325, 0.55557, 0.565732, 0.575808, 0.585798, 0.595699, 0.605511, 0.615232, 0.624859, 0.634393, 0.643832, 0.653173, 0.662416, 0.671559, 0.680601, 0.689541, 0.698376, 0.707107, 0.715731, 0.724247, 0.732654, 0.740951, 0.749136, 0.757209, 0.765167, 0.77301, 0.780737, 0.788346, 0.795837, 0.803208, 0.810457, 0.817585, 0.824589, 0.83147, 0.838225, 0.844854, 0.851355, 0.857729, 0.863973, 0.870087, 0.87607, 0.881921, 0.88764, 0.893224, 0.898674, 0.903989, 0.909168, 0.91421, 0.919114, 0.92388, 0.928506, 0.932993, 0.937339, 0.941544, 0.945607, 0.949528, 0.953306, 0.95694, 0.960431, 0.963776, 0.966976, 0.970031, 0.97294, 0.975702, 0.978317, 0.980785, 0.983105, 0.985278, 0.987301, 0.989177, 0.990903, 0.99248, 0.993907, 0.995185, 0.996313, 0.99729, 0.998118, 0.998795, 0.999322, 0.999699, 0.999925, 0, 0.00306796, 0.00613588, 0.00920375, 0.0122715, 0.0153392, 0.0184067, 0.0214741, 0.0245412, 0.0276081, 0.0306748, 0.0337412, 0.0368072, 0.0398729, 0.0429383, 0.0460032, 0.0490677, 0.0521317, 0.0551952, 0.0582583, 0.0613207, 0.0643826, 0.0674439, 0.0705046, 0.0735646, 0.0766239, 0.0796824, 0.0827403, 0.0857973, 0.0888536, 0.091909, 0.0949635, 0.0980171, 0.10107, 0.104122, 0.107172, 0.110222, 0.113271, 0.116319, 0.119365, 0.122411, 0.125455, 0.128498, 0.13154, 0.134581, 0.13762, 0.140658, 0.143695, 0.14673, 0.149765, 0.152797, 0.155828, 0.158858, 0.161886, 0.164913, 0.167938, 0.170962, 0.173984, 0.177004, 0.180023, 0.18304, 0.186055, 0.189069, 0.19208, 0.19509, 0.198098, 0.201105, 0.204109, 0.207111, 0.210112, 0.21311, 0.216107, 0.219101, 0.222094, 0.225084, 0.228072, 0.231058, 0.234042, 0.237024, 0.240003, 0.24298, 0.245955, 0.248928, 0.251898, 0.254866, 0.257831, 0.260794, 0.263755, 0.266713, 0.269668, 0.272621, 0.275572, 0.27852, 0.281465, 0.284408, 0.287347, 0.290285, 0.293219, 0.296151, 0.29908, 0.302006, 0.304929, 0.30785, 0.310767, 0.313682, 0.316593, 0.319502, 0.322408, 0.32531, 0.32821, 0.331106, 0.334, 0.33689, 0.339777, 0.342661, 0.345541, 0.348419, 0.351293, 0.354164, 0.357031, 0.359895, 0.362756, 0.365613, 0.368467, 0.371317, 0.374164, 0.377007, 0.379847, 0.382683, 0.385516, 0.388345, 0.39117, 0.393992, 0.39681, 0.399624, 0.402435, 0.405241, 0.408044, 0.410843, 0.413638, 0.41643, 0.419217, 0.422, 0.42478, 0.427555, 0.430326, 0.433094, 0.435857, 0.438616, 0.441371, 0.444122, 0.446869, 0.449611, 0.45235, 0.455084, 0.457813, 0.460539, 0.46326, 0.465976, 0.468689, 0.471397, 0.4741, 0.476799, 0.479494, 0.482184, 0.484869, 0.48755, 0.490226, 0.492898, 0.495565, 0.498228, 0.500885, 0.503538, 0.506187, 0.50883, 0.511469, 0.514103, 0.516732, 0.519356, 0.521975, 0.52459, 0.527199, 0.529804, 0.532403, 0.534998, 0.537587, 0.540171, 0.542751, 0.545325, 0.547894, 0.550458, 0.553017, 0.55557, 0.558119, 0.560662, 0.563199, 0.565732, 0.568259, 0.570781, 0.573297, 0.575808, 0.578314, 0.580814, 0.583309, 0.585798, 0.588282, 0.59076, 0.593232, 0.595699, 0.598161, 0.600616, 0.603067, 0.605511, 0.60795, 0.610383, 0.61281, 0.615232, 0.617647, 0.620057, 0.622461, 0.624859, 0.627252, 0.629638, 0.632019, 0.634393, 0.636762, 0.639124, 0.641481, 0.643832, 0.646176, 0.648514, 0.650847, 0.653173, 0.655493, 0.657807, 0.660114, 0.662416, 0.664711, 0.667, 0.669283, 0.671559, 0.673829, 0.676093, 0.67835, 0.680601, 0.682846, 0.685084, 0.687315, 0.689541, 0.691759, 0.693971, 0.696177, 0.698376, 0.700569, 0.702755, 0.704934, 0.707107, 0.709273, 0.711432, 0.713585, 0.715731, 0.71787, 0.720003, 0.722128, 0.724247, 0.726359, 0.728464, 0.730563, 0.732654, 0.734739, 0.736817, 0.738887, 0.740951, 0.743008, 0.745058, 0.747101, 0.749136, 0.751165, 0.753187, 0.755201, 0.757209, 0.759209, 0.761202, 0.763188, 0.765167, 0.767139, 0.769103, 0.771061, 0.77301, 0.774953, 0.776888, 0.778817, 0.780737, 0.782651, 0.784557, 0.786455, 0.788346, 0.79023, 0.792107, 0.793975, 0.795837, 0.797691, 0.799537, 0.801376, 0.803208, 0.805031, 0.806848, 0.808656, 0.810457, 0.812251, 0.814036, 0.815814, 0.817585, 0.819348, 0.821103, 0.82285, 0.824589, 0.826321, 0.828045, 0.829761, 0.83147, 0.83317, 0.834863, 0.836548, 0.838225, 0.839894, 0.841555, 0.843208, 0.844854, 0.846491, 0.84812, 0.849742, 0.851355, 0.852961, 0.854558, 0.856147, 0.857729, 0.859302, 0.860867, 0.862424, 0.863973, 0.865514, 0.867046, 0.868571, 0.870087, 0.871595, 0.873095, 0.874587, 0.87607, 0.877545, 0.879012, 0.880471, 0.881921, 0.883363, 0.884797, 0.886223, 0.88764, 0.889048, 0.890449, 0.891841, 0.893224, 0.894599, 0.895966, 0.897325, 0.898674, 0.900016, 0.901349, 0.902673, 0.903989, 0.905297, 0.906596, 0.907886, 0.909168, 0.910441, 0.911706, 0.912962, 0.91421, 0.915449, 0.916679, 0.917901, 0.919114, 0.920318, 0.921514, 0.922701, 0.92388, 0.925049, 0.92621, 0.927363, 0.928506, 0.929641, 0.930767, 0.931884, 0.932993, 0.934093, 0.935184, 0.936266, 0.937339, 0.938404, 0.939459, 0.940506, 0.941544, 0.942573, 0.943593, 0.944605, 0.945607, 0.946601, 0.947586, 0.948561, 0.949528, 0.950486, 0.951435, 0.952375, 0.953306, 0.954228, 0.955141, 0.956045, 0.95694, 0.957826, 0.958703, 0.959572, 0.960431, 0.96128, 0.962121, 0.962953, 0.963776, 0.96459, 0.965394, 0.96619, 0.966976, 0.967754, 0.968522, 0.969281, 0.970031, 0.970772, 0.971504, 0.972226, 0.97294, 0.973644, 0.974339, 0.975025, 0.975702, 0.97637, 0.977028, 0.977677, 0.978317, 0.978948, 0.97957, 0.980182, 0.980785, 0.981379, 0.981964, 0.982539, 0.983105, 0.983662, 0.98421, 0.984749, 0.985278, 0.985798, 0.986308, 0.986809, 0.987301, 0.987784, 0.988258, 0.988722, 0.989177, 0.989622, 0.990058, 0.990485, 0.990903, 0.991311, 0.99171, 0.992099, 0.99248, 0.99285, 0.993212, 0.993564, 0.993907, 0.99424, 0.994565, 0.994879, 0.995185, 0.995481, 0.995767, 0.996045, 0.996313, 0.996571, 0.99682, 0.99706, 0.99729, 0.997511, 0.997723, 0.997925, 0.998118, 0.998302, 0.998476, 0.99864, 0.998795, 0.998941, 0.999078, 0.999205, 0.999322, 0.999431, 0.999529, 0.999619, 0.999699, 0.999769, 0.999831, 0.999882, 0.999925, 0.999958, 0.999981, 0.999995};
+
+double allTheC2 [680] = {1, 0.92388, 0.707107, 0.382683, 2.22045e-16, -0.382683, -0.707107, -0.92388, 1, 0.995185, 0.980785, 0.95694, 0.92388, 0.881921, 0.83147, 0.77301, 0.707107, 0.634393, 0.55557, 0.471397, 0.382683, 0.290285, 0.19509, 0.0980171, -3.33067e-16, -0.0980171, -0.19509, -0.290285, -0.382683, -0.471397, -0.55557, -0.634393, -0.707107, -0.77301, -0.83147, -0.881921, -0.92388, -0.95694, -0.980785, -0.995185, 1, 0.999699, 0.998795, 0.99729, 0.995185, 0.99248, 0.989177, 0.985278, 0.980785, 0.975702, 0.970031, 0.963776, 0.95694, 0.949528, 0.941544, 0.932993, 0.92388, 0.91421, 0.903989, 0.893224, 0.881921, 0.870087, 0.857729, 0.844854, 0.83147, 0.817585, 0.803208, 0.788346, 0.77301, 0.757209, 0.740951, 0.724247, 0.707107, 0.689541, 0.671559, 0.653173, 0.634393, 0.615232, 0.595699, 0.575808, 0.55557, 0.534998, 0.514103, 0.492898, 0.471397, 0.449611, 0.427555, 0.405241, 0.382683, 0.359895, 0.33689, 0.313682, 0.290285, 0.266713, 0.24298, 0.219101, 0.19509, 0.170962, 0.14673, 0.122411, 0.0980171, 0.0735646, 0.0490677, 0.0245412, 1.11022e-15, -0.0245412, -0.0490677, -0.0735646, -0.0980171, -0.122411, -0.14673, -0.170962, -0.19509, -0.219101, -0.24298, -0.266713, -0.290285, -0.313682, -0.33689, -0.359895, -0.382683, -0.405241, -0.427555, -0.449611, -0.471397, -0.492898, -0.514103, -0.534998, -0.55557, -0.575808, -0.595699, -0.615232, -0.634393, -0.653173, -0.671559, -0.689541, -0.707107, -0.724247, -0.740951, -0.757209, -0.77301, -0.788346, -0.803208, -0.817585, -0.83147, -0.844854, -0.857729, -0.870087, -0.881921, -0.893224, -0.903989, -0.91421, -0.92388, -0.932993, -0.941544, -0.949528, -0.95694, -0.963776, -0.970031, -0.975702, -0.980785, -0.985278, -0.989177, -0.99248, -0.995185, -0.99729, -0.998795, -0.999699, 1, 0.999981, 0.999925, 0.999831, 0.999699, 0.999529, 0.999322, 0.999078, 0.998795, 0.998476, 0.998118, 0.997723, 0.99729, 0.99682, 0.996313, 0.995767, 0.995185, 0.994565, 0.993907, 0.993212, 0.99248, 0.99171, 0.990903, 0.990058, 0.989177, 0.988258, 0.987301, 0.986308, 0.985278, 0.98421, 0.983105, 0.981964, 0.980785, 0.97957, 0.978317, 0.977028, 0.975702, 0.974339, 0.97294, 0.971504, 0.970031, 0.968522, 0.966976, 0.965394, 0.963776, 0.962121, 0.960431, 0.958703, 0.95694, 0.955141, 0.953306, 0.951435, 0.949528, 0.947586, 0.945607, 0.943593, 0.941544, 0.939459, 0.937339, 0.935184, 0.932993, 0.930767, 0.928506, 0.92621, 0.92388, 0.921514, 0.919114, 0.916679, 0.91421, 0.911706, 0.909168, 0.906596, 0.903989, 0.901349, 0.898674, 0.895966, 0.893224, 0.890449, 0.88764, 0.884797, 0.881921, 0.879012, 0.87607, 0.873095, 0.870087, 0.867046, 0.863973, 0.860867, 0.857729, 0.854558, 0.851355, 0.84812, 0.844854, 0.841555, 0.838225, 0.834863, 0.83147, 0.828045, 0.824589, 0.821103, 0.817585, 0.814036, 0.810457, 0.806848, 0.803208, 0.799537, 0.795837, 0.792107, 0.788346, 0.784557, 0.780737, 0.776888, 0.77301, 0.769103, 0.765167, 0.761202, 0.757209, 0.753187, 0.749136, 0.745058, 0.740951, 0.736817, 0.732654, 0.728464, 0.724247, 0.720003, 0.715731, 0.711432, 0.707107, 0.702755, 0.698376, 0.693971, 0.689541, 0.685084, 0.680601, 0.676093, 0.671559, 0.667, 0.662416, 0.657807, 0.653173, 0.648514, 0.643832, 0.639124, 0.634393, 0.629638, 0.624859, 0.620057, 0.615232, 0.610383, 0.605511, 0.600616, 0.595699, 0.59076, 0.585798, 0.580814, 0.575808, 0.570781, 0.565732, 0.560662, 0.55557, 0.550458, 0.545325, 0.540171, 0.534998, 0.529804, 0.52459, 0.519356, 0.514103, 0.50883, 0.503538, 0.498228, 0.492898, 0.48755, 0.482184, 0.476799, 0.471397, 0.465976, 0.460539, 0.455084, 0.449611, 0.444122, 0.438616, 0.433094, 0.427555, 0.422, 0.41643, 0.410843, 0.405241, 0.399624, 0.393992, 0.388345, 0.382683, 0.377007, 0.371317, 0.365613, 0.359895, 0.354164, 0.348419, 0.342661, 0.33689, 0.331106, 0.32531, 0.319502, 0.313682, 0.30785, 0.302006, 0.296151, 0.290285, 0.284408, 0.27852, 0.272621, 0.266713, 0.260794, 0.254866, 0.248928, 0.24298, 0.237024, 0.231058, 0.225084, 0.219101, 0.21311, 0.207111, 0.201105, 0.19509, 0.189069, 0.18304, 0.177004, 0.170962, 0.164913, 0.158858, 0.152797, 0.14673, 0.140658, 0.134581, 0.128498, 0.122411, 0.116319, 0.110222, 0.104122, 0.0980171, 0.091909, 0.0857973, 0.0796824, 0.0735646, 0.0674439, 0.0613207, 0.0551952, 0.0490677, 0.0429383, 0.0368072, 0.0306748, 0.0245412, 0.0184067, 0.0122715, 0.00613588, 4.21885e-15, -0.00613588, -0.0122715, -0.0184067, -0.0245412, -0.0306748, -0.0368072, -0.0429383, -0.0490677, -0.0551952, -0.0613207, -0.0674439, -0.0735646, -0.0796824, -0.0857973, -0.091909, -0.0980171, -0.104122, -0.110222, -0.116319, -0.122411, -0.128498, -0.134581, -0.140658, -0.14673, -0.152797, -0.158858, -0.164913, -0.170962, -0.177004, -0.18304, -0.189069, -0.19509, -0.201105, -0.207111, -0.21311, -0.219101, -0.225084, -0.231058, -0.237024, -0.24298, -0.248928, -0.254866, -0.260794, -0.266713, -0.272621, -0.27852, -0.284408, -0.290285, -0.296151, -0.302006, -0.30785, -0.313682, -0.319502, -0.32531, -0.331106, -0.33689, -0.342661, -0.348419, -0.354164, -0.359895, -0.365613, -0.371317, -0.377007, -0.382683, -0.388345, -0.393992, -0.399624, -0.405241, -0.410843, -0.41643, -0.422, -0.427555, -0.433094, -0.438616, -0.444122, -0.449611, -0.455084, -0.460539, -0.465976, -0.471397, -0.476799, -0.482184, -0.48755, -0.492898, -0.498228, -0.503538, -0.50883, -0.514103, -0.519356, -0.52459, -0.529804, -0.534998, -0.540171, -0.545325, -0.550458, -0.55557, -0.560662, -0.565732, -0.570781, -0.575808, -0.580814, -0.585798, -0.59076, -0.595699, -0.600616, -0.605511, -0.610383, -0.615232, -0.620057, -0.624859, -0.629638, -0.634393, -0.639124, -0.643832, -0.648514, -0.653173, -0.657807, -0.662416, -0.667, -0.671559, -0.676093, -0.680601, -0.685084, -0.689541, -0.693971, -0.698376, -0.702755, -0.707107, -0.711432, -0.715731, -0.720003, -0.724247, -0.728464, -0.732654, -0.736817, -0.740951, -0.745058, -0.749136, -0.753187, -0.757209, -0.761202, -0.765167, -0.769103, -0.77301, -0.776888, -0.780737, -0.784557, -0.788346, -0.792107, -0.795837, -0.799537, -0.803208, -0.806848, -0.810457, -0.814036, -0.817585, -0.821103, -0.824589, -0.828045, -0.83147, -0.834863, -0.838225, -0.841555, -0.844854, -0.84812, -0.851355, -0.854558, -0.857729, -0.860867, -0.863973, -0.867046, -0.870087, -0.873095, -0.87607, -0.879012, -0.881921, -0.884797, -0.88764, -0.890449, -0.893224, -0.895966, -0.898674, -0.901349, -0.903989, -0.906596, -0.909168, -0.911706, -0.91421, -0.916679, -0.919114, -0.921514, -0.92388, -0.92621, -0.928506, -0.930767, -0.932993, -0.935184, -0.937339, -0.939459, -0.941544, -0.943593, -0.945607, -0.947586, -0.949528, -0.951435, -0.953306, -0.955141, -0.95694, -0.958703, -0.960431, -0.962121, -0.963776, -0.965394, -0.966976, -0.968522, -0.970031, -0.971504, -0.97294, -0.974339, -0.975702, -0.977028, -0.978317, -0.97957, -0.980785, -0.981964, -0.983105, -0.98421, -0.985278, -0.986308, -0.987301, -0.988258, -0.989177, -0.990058, -0.990903, -0.99171, -0.99248, -0.993212, -0.993907, -0.994565, -0.995185, -0.995767, -0.996313, -0.99682, -0.99729, -0.997723, -0.998118, -0.998476, -0.998795, -0.999078, -0.999322, -0.999529, -0.999699, -0.999831, -0.999925, -0.999981};
+
+double allTheS2 [680] = {0, 0.382683, 0.707107, 0.92388, 1, 0.92388, 0.707107, 0.382683, 0, 0.0980171, 0.19509, 0.290285, 0.382683, 0.471397, 0.55557, 0.634393, 0.707107, 0.77301, 0.83147, 0.881921, 0.92388, 0.95694, 0.980785, 0.995185, 1, 0.995185, 0.980785, 0.95694, 0.92388, 0.881921, 0.83147, 0.77301, 0.707107, 0.634393, 0.55557, 0.471397, 0.382683, 0.290285, 0.19509, 0.0980171, 0, 0.0245412, 0.0490677, 0.0735646, 0.0980171, 0.122411, 0.14673, 0.170962, 0.19509, 0.219101, 0.24298, 0.266713, 0.290285, 0.313682, 0.33689, 0.359895, 0.382683, 0.405241, 0.427555, 0.449611, 0.471397, 0.492898, 0.514103, 0.534998, 0.55557, 0.575808, 0.595699, 0.615232, 0.634393, 0.653173, 0.671559, 0.689541, 0.707107, 0.724247, 0.740951, 0.757209, 0.77301, 0.788346, 0.803208, 0.817585, 0.83147, 0.844854, 0.857729, 0.870087, 0.881921, 0.893224, 0.903989, 0.91421, 0.92388, 0.932993, 0.941544, 0.949528, 0.95694, 0.963776, 0.970031, 0.975702, 0.980785, 0.985278, 0.989177, 0.99248, 0.995185, 0.99729, 0.998795, 0.999699, 1, 0.999699, 0.998795, 0.99729, 0.995185, 0.99248, 0.989177, 0.985278, 0.980785, 0.975702, 0.970031, 0.963776, 0.95694, 0.949528, 0.941544, 0.932993, 0.92388, 0.91421, 0.903989, 0.893224, 0.881921, 0.870087, 0.857729, 0.844854, 0.83147, 0.817585, 0.803208, 0.788346, 0.77301, 0.757209, 0.740951, 0.724247, 0.707107, 0.689541, 0.671559, 0.653173, 0.634393, 0.615232, 0.595699, 0.575808, 0.55557, 0.534998, 0.514103, 0.492898, 0.471397, 0.449611, 0.427555, 0.405241, 0.382683, 0.359895, 0.33689, 0.313682, 0.290285, 0.266713, 0.24298, 0.219101, 0.19509, 0.170962, 0.14673, 0.122411, 0.0980171, 0.0735646, 0.0490677, 0.0245412, 0, 0.00613588, 0.0122715, 0.0184067, 0.0245412, 0.0306748, 0.0368072, 0.0429383, 0.0490677, 0.0551952, 0.0613207, 0.0674439, 0.0735646, 0.0796824, 0.0857973, 0.091909, 0.0980171, 0.104122, 0.110222, 0.116319, 0.122411, 0.128498, 0.134581, 0.140658, 0.14673, 0.152797, 0.158858, 0.164913, 0.170962, 0.177004, 0.18304, 0.189069, 0.19509, 0.201105, 0.207111, 0.21311, 0.219101, 0.225084, 0.231058, 0.237024, 0.24298, 0.248928, 0.254866, 0.260794, 0.266713, 0.272621, 0.27852, 0.284408, 0.290285, 0.296151, 0.302006, 0.30785, 0.313682, 0.319502, 0.32531, 0.331106, 0.33689, 0.342661, 0.348419, 0.354164, 0.359895, 0.365613, 0.371317, 0.377007, 0.382683, 0.388345, 0.393992, 0.399624, 0.405241, 0.410843, 0.41643, 0.422, 0.427555, 0.433094, 0.438616, 0.444122, 0.449611, 0.455084, 0.460539, 0.465976, 0.471397, 0.476799, 0.482184, 0.48755, 0.492898, 0.498228, 0.503538, 0.50883, 0.514103, 0.519356, 0.52459, 0.529804, 0.534998, 0.540171, 0.545325, 0.550458, 0.55557, 0.560662, 0.565732, 0.570781, 0.575808, 0.580814, 0.585798, 0.59076, 0.595699, 0.600616, 0.605511, 0.610383, 0.615232, 0.620057, 0.624859, 0.629638, 0.634393, 0.639124, 0.643832, 0.648514, 0.653173, 0.657807, 0.662416, 0.667, 0.671559, 0.676093, 0.680601, 0.685084, 0.689541, 0.693971, 0.698376, 0.702755, 0.707107, 0.711432, 0.715731, 0.720003, 0.724247, 0.728464, 0.732654, 0.736817, 0.740951, 0.745058, 0.749136, 0.753187, 0.757209, 0.761202, 0.765167, 0.769103, 0.77301, 0.776888, 0.780737, 0.784557, 0.788346, 0.792107, 0.795837, 0.799537, 0.803208, 0.806848, 0.810457, 0.814036, 0.817585, 0.821103, 0.824589, 0.828045, 0.83147, 0.834863, 0.838225, 0.841555, 0.844854, 0.84812, 0.851355, 0.854558, 0.857729, 0.860867, 0.863973, 0.867046, 0.870087, 0.873095, 0.87607, 0.879012, 0.881921, 0.884797, 0.88764, 0.890449, 0.893224, 0.895966, 0.898674, 0.901349, 0.903989, 0.906596, 0.909168, 0.911706, 0.91421, 0.916679, 0.919114, 0.921514, 0.92388, 0.92621, 0.928506, 0.930767, 0.932993, 0.935184, 0.937339, 0.939459, 0.941544, 0.943593, 0.945607, 0.947586, 0.949528, 0.951435, 0.953306, 0.955141, 0.95694, 0.958703, 0.960431, 0.962121, 0.963776, 0.965394, 0.966976, 0.968522, 0.970031, 0.971504, 0.97294, 0.974339, 0.975702, 0.977028, 0.978317, 0.97957, 0.980785, 0.981964, 0.983105, 0.98421, 0.985278, 0.986308, 0.987301, 0.988258, 0.989177, 0.990058, 0.990903, 0.99171, 0.99248, 0.993212, 0.993907, 0.994565, 0.995185, 0.995767, 0.996313, 0.99682, 0.99729, 0.997723, 0.998118, 0.998476, 0.998795, 0.999078, 0.999322, 0.999529, 0.999699, 0.999831, 0.999925, 0.999981, 1, 0.999981, 0.999925, 0.999831, 0.999699, 0.999529, 0.999322, 0.999078, 0.998795, 0.998476, 0.998118, 0.997723, 0.99729, 0.99682, 0.996313, 0.995767, 0.995185, 0.994565, 0.993907, 0.993212, 0.99248, 0.99171, 0.990903, 0.990058, 0.989177, 0.988258, 0.987301, 0.986308, 0.985278, 0.98421, 0.983105, 0.981964, 0.980785, 0.97957, 0.978317, 0.977028, 0.975702, 0.974339, 0.97294, 0.971504, 0.970031, 0.968522, 0.966976, 0.965394, 0.963776, 0.962121, 0.960431, 0.958703, 0.95694, 0.955141, 0.953306, 0.951435, 0.949528, 0.947586, 0.945607, 0.943593, 0.941544, 0.939459, 0.937339, 0.935184, 0.932993, 0.930767, 0.928506, 0.92621, 0.92388, 0.921514, 0.919114, 0.916679, 0.91421, 0.911706, 0.909168, 0.906596, 0.903989, 0.901349, 0.898674, 0.895966, 0.893224, 0.890449, 0.88764, 0.884797, 0.881921, 0.879012, 0.87607, 0.873095, 0.870087, 0.867046, 0.863973, 0.860867, 0.857729, 0.854558, 0.851355, 0.84812, 0.844854, 0.841555, 0.838225, 0.834863, 0.83147, 0.828045, 0.824589, 0.821103, 0.817585, 0.814036, 0.810457, 0.806848, 0.803208, 0.799537, 0.795837, 0.792107, 0.788346, 0.784557, 0.780737, 0.776888, 0.77301, 0.769103, 0.765167, 0.761202, 0.757209, 0.753187, 0.749136, 0.745058, 0.740951, 0.736817, 0.732654, 0.728464, 0.724247, 0.720003, 0.715731, 0.711432, 0.707107, 0.702755, 0.698376, 0.693971, 0.689541, 0.685084, 0.680601, 0.676093, 0.671559, 0.667, 0.662416, 0.657807, 0.653173, 0.648514, 0.643832, 0.639124, 0.634393, 0.629638, 0.624859, 0.620057, 0.615232, 0.610383, 0.605511, 0.600616, 0.595699, 0.59076, 0.585798, 0.580814, 0.575808, 0.570781, 0.565732, 0.560662, 0.55557, 0.550458, 0.545325, 0.540171, 0.534998, 0.529804, 0.52459, 0.519356, 0.514103, 0.50883, 0.503538, 0.498228, 0.492898, 0.48755, 0.482184, 0.476799, 0.471397, 0.465976, 0.460539, 0.455084, 0.449611, 0.444122, 0.438616, 0.433094, 0.427555, 0.422, 0.41643, 0.410843, 0.405241, 0.399624, 0.393992, 0.388345, 0.382683, 0.377007, 0.371317, 0.365613, 0.359895, 0.354164, 0.348419, 0.342661, 0.33689, 0.331106, 0.32531, 0.319502, 0.313682, 0.30785, 0.302006, 0.296151, 0.290285, 0.284408, 0.27852, 0.272621, 0.266713, 0.260794, 0.254866, 0.248928, 0.24298, 0.237024, 0.231058, 0.225084, 0.219101, 0.21311, 0.207111, 0.201105, 0.19509, 0.189069, 0.18304, 0.177004, 0.170962, 0.164913, 0.158858, 0.152797, 0.14673, 0.140658, 0.134581, 0.128498, 0.122411, 0.116319, 0.110222, 0.104122, 0.0980171, 0.091909, 0.0857973, 0.0796824, 0.0735646, 0.0674439, 0.0613207, 0.0551952, 0.0490677, 0.0429383, 0.0368072, 0.0306748, 0.0245412, 0.0184067, 0.0122715, 0.00613588};
+
+double allTheC3 [680] = {1, 0.83147, 0.382683, -0.19509, -0.707107, -0.980785, -0.92388, -0.55557, 1, 0.989177, 0.95694, 0.903989, 0.83147, 0.740951, 0.634393, 0.514103, 0.382683, 0.24298, 0.0980171, -0.0490677, -0.19509, -0.33689, -0.471397, -0.595699, -0.707107, -0.803208, -0.881921, -0.941544, -0.980785, -0.998795, -0.995185, -0.970031, -0.92388, -0.857729, -0.77301, -0.671559, -0.55557, -0.427555, -0.290285, -0.14673, 1, 0.999322, 0.99729, 0.993907, 0.989177, 0.983105, 0.975702, 0.966976, 0.95694, 0.945607, 0.932993, 0.919114, 0.903989, 0.88764, 0.870087, 0.851355, 0.83147, 0.810457, 0.788346, 0.765167, 0.740951, 0.715731, 0.689541, 0.662416, 0.634393, 0.605511, 0.575808, 0.545325, 0.514103, 0.482184, 0.449611, 0.41643, 0.382683, 0.348419, 0.313682, 0.27852, 0.24298, 0.207111, 0.170962, 0.134581, 0.0980171, 0.0613207, 0.0245412, -0.0122715, -0.0490677, -0.0857973, -0.122411, -0.158858, -0.19509, -0.231058, -0.266713, -0.302006, -0.33689, -0.371317, -0.405241, -0.438616, -0.471397, -0.503538, -0.534998, -0.565732, -0.595699, -0.624859, -0.653173, -0.680601, -0.707107, -0.732654, -0.757209, -0.780737, -0.803208, -0.824589, -0.844854, -0.863973, -0.881921, -0.898674, -0.91421, -0.928506, -0.941544, -0.953306, -0.963776, -0.97294, -0.980785, -0.987301, -0.99248, -0.996313, -0.998795, -0.999925, -0.999699, -0.998118, -0.995185, -0.990903, -0.985278, -0.978317, -0.970031, -0.960431, -0.949528, -0.937339, -0.92388, -0.909168, -0.893224, -0.87607, -0.857729, -0.838225, -0.817585, -0.795837, -0.77301, -0.749136, -0.724247, -0.698376, -0.671559, -0.643832, -0.615232, -0.585798, -0.55557, -0.52459, -0.492898, -0.460539, -0.427555, -0.393992, -0.359895, -0.32531, -0.290285, -0.254866, -0.219101, -0.18304, -0.14673, -0.110222, -0.0735646, -0.0368072, 1, 0.999958, 0.999831, 0.999619, 0.999322, 0.998941, 0.998476, 0.997925, 0.99729, 0.996571, 0.995767, 0.994879, 0.993907, 0.99285, 0.99171, 0.990485, 0.989177, 0.987784, 0.986308, 0.984749, 0.983105, 0.981379, 0.97957, 0.977677, 0.975702, 0.973644, 0.971504, 0.969281, 0.966976, 0.96459, 0.962121, 0.959572, 0.95694, 0.954228, 0.951435, 0.948561, 0.945607, 0.942573, 0.939459, 0.936266, 0.932993, 0.929641, 0.92621, 0.922701, 0.919114, 0.915449, 0.911706, 0.907886, 0.903989, 0.900016, 0.895966, 0.891841, 0.88764, 0.883363, 0.879012, 0.874587, 0.870087, 0.865514, 0.860867, 0.856147, 0.851355, 0.846491, 0.841555, 0.836548, 0.83147, 0.826321, 0.821103, 0.815814, 0.810457, 0.805031, 0.799537, 0.793975, 0.788346, 0.782651, 0.776888, 0.771061, 0.765167, 0.759209, 0.753187, 0.747101, 0.740951, 0.734739, 0.728464, 0.722128, 0.715731, 0.709273, 0.702755, 0.696177, 0.689541, 0.682846, 0.676093, 0.669283, 0.662416, 0.655493, 0.648514, 0.641481, 0.634393, 0.627252, 0.620057, 0.61281, 0.605511, 0.598161, 0.59076, 0.583309, 0.575808, 0.568259, 0.560662, 0.553017, 0.545325, 0.537587, 0.529804, 0.521975, 0.514103, 0.506187, 0.498228, 0.490226, 0.482184, 0.4741, 0.465976, 0.457813, 0.449611, 0.441371, 0.433094, 0.42478, 0.41643, 0.408044, 0.399624, 0.39117, 0.382683, 0.374164, 0.365613, 0.357031, 0.348419, 0.339777, 0.331106, 0.322408, 0.313682, 0.304929, 0.296151, 0.287347, 0.27852, 0.269668, 0.260794, 0.251898, 0.24298, 0.234042, 0.225084, 0.216107, 0.207111, 0.198098, 0.189069, 0.180023, 0.170962, 0.161886, 0.152797, 0.143695, 0.134581, 0.125455, 0.116319, 0.107172, 0.0980171, 0.0888536, 0.0796824, 0.0705046, 0.0613207, 0.0521317, 0.0429383, 0.0337412, 0.0245412, 0.0153392, 0.00613588, -0.00306796, -0.0122715, -0.0214741, -0.0306748, -0.0398729, -0.0490677, -0.0582583, -0.0674439, -0.0766239, -0.0857973, -0.0949635, -0.104122, -0.113271, -0.122411, -0.13154, -0.140658, -0.149765, -0.158858, -0.167938, -0.177004, -0.186055, -0.19509, -0.204109, -0.21311, -0.222094, -0.231058, -0.240003, -0.248928, -0.257831, -0.266713, -0.275572, -0.284408, -0.293219, -0.302006, -0.310767, -0.319502, -0.32821, -0.33689, -0.345541, -0.354164, -0.362756, -0.371317, -0.379847, -0.388345, -0.39681, -0.405241, -0.413638, -0.422, -0.430326, -0.438616, -0.446869, -0.455084, -0.46326, -0.471397, -0.479494, -0.48755, -0.495565, -0.503538, -0.511469, -0.519356, -0.527199, -0.534998, -0.542751, -0.550458, -0.558119, -0.565732, -0.573297, -0.580814, -0.588282, -0.595699, -0.603067, -0.610383, -0.617647, -0.624859, -0.632019, -0.639124, -0.646176, -0.653173, -0.660114, -0.667, -0.673829, -0.680601, -0.687315, -0.693971, -0.700569, -0.707107, -0.713585, -0.720003, -0.726359, -0.732654, -0.738887, -0.745058, -0.751165, -0.757209, -0.763188, -0.769103, -0.774953, -0.780737, -0.786455, -0.792107, -0.797691, -0.803208, -0.808656, -0.814036, -0.819348, -0.824589, -0.829761, -0.834863, -0.839894, -0.844854, -0.849742, -0.854558, -0.859302, -0.863973, -0.868571, -0.873095, -0.877545, -0.881921, -0.886223, -0.890449, -0.894599, -0.898674, -0.902673, -0.906596, -0.910441, -0.91421, -0.917901, -0.921514, -0.925049, -0.928506, -0.931884, -0.935184, -0.938404, -0.941544, -0.944605, -0.947586, -0.950486, -0.953306, -0.956045, -0.958703, -0.96128, -0.963776, -0.96619, -0.968522, -0.970772, -0.97294, -0.975025, -0.977028, -0.978948, -0.980785, -0.982539, -0.98421, -0.985798, -0.987301, -0.988722, -0.990058, -0.991311, -0.99248, -0.993564, -0.994565, -0.995481, -0.996313, -0.99706, -0.997723, -0.998302, -0.998795, -0.999205, -0.999529, -0.999769, -0.999925, -0.999995, -0.999981, -0.999882, -0.999699, -0.999431, -0.999078, -0.99864, -0.998118, -0.997511, -0.99682, -0.996045, -0.995185, -0.99424, -0.993212, -0.992099, -0.990903, -0.989622, -0.988258, -0.986809, -0.985278, -0.983662, -0.981964, -0.980182, -0.978317, -0.97637, -0.974339, -0.972226, -0.970031, -0.967754, -0.965394, -0.962953, -0.960431, -0.957826, -0.955141, -0.952375, -0.949528, -0.946601, -0.943593, -0.940506, -0.937339, -0.934093, -0.930767, -0.927363, -0.92388, -0.920318, -0.916679, -0.912962, -0.909168, -0.905297, -0.901349, -0.897325, -0.893224, -0.889048, -0.884797, -0.880471, -0.87607, -0.871595, -0.867046, -0.862424, -0.857729, -0.852961, -0.84812, -0.843208, -0.838225, -0.83317, -0.828045, -0.82285, -0.817585, -0.812251, -0.806848, -0.801376, -0.795837, -0.79023, -0.784557, -0.778817, -0.77301, -0.767139, -0.761202, -0.755201, -0.749136, -0.743008, -0.736817, -0.730563, -0.724247, -0.71787, -0.711432, -0.704934, -0.698376, -0.691759, -0.685084, -0.67835, -0.671559, -0.664711, -0.657807, -0.650847, -0.643832, -0.636762, -0.629638, -0.622461, -0.615232, -0.60795, -0.600616, -0.593232, -0.585798, -0.578314, -0.570781, -0.563199, -0.55557, -0.547894, -0.540171, -0.532403, -0.52459, -0.516732, -0.50883, -0.500885, -0.492898, -0.484869, -0.476799, -0.468689, -0.460539, -0.45235, -0.444122, -0.435857, -0.427555, -0.419217, -0.410843, -0.402435, -0.393992, -0.385516, -0.377007, -0.368467, -0.359895, -0.351293, -0.342661, -0.334, -0.32531, -0.316593, -0.30785, -0.29908, -0.290285, -0.281465, -0.272621, -0.263755, -0.254866, -0.245955, -0.237024, -0.228072, -0.219101, -0.210112, -0.201105, -0.19208, -0.18304, -0.173984, -0.164913, -0.155828, -0.14673, -0.13762, -0.128498, -0.119365, -0.110222, -0.10107, -0.091909, -0.0827403, -0.0735646, -0.0643826, -0.0551952, -0.0460032, -0.0368072, -0.0276081, -0.0184067, -0.00920375};
+
+double allTheS3 [680] = {0, 0.55557, 0.92388, 0.980785, 0.707107, 0.19509, -0.382683, -0.83147, 0, 0.14673, 0.290285, 0.427555, 0.55557, 0.671559, 0.77301, 0.857729, 0.92388, 0.970031, 0.995185, 0.998795, 0.980785, 0.941544, 0.881921, 0.803208, 0.707107, 0.595699, 0.471397, 0.33689, 0.19509, 0.0490677, -0.0980171, -0.24298, -0.382683, -0.514103, -0.634393, -0.740951, -0.83147, -0.903989, -0.95694, -0.989177, 0, 0.0368072, 0.0735646, 0.110222, 0.14673, 0.18304, 0.219101, 0.254866, 0.290285, 0.32531, 0.359895, 0.393992, 0.427555, 0.460539, 0.492898, 0.52459, 0.55557, 0.585798, 0.615232, 0.643832, 0.671559, 0.698376, 0.724247, 0.749136, 0.77301, 0.795837, 0.817585, 0.838225, 0.857729, 0.87607, 0.893224, 0.909168, 0.92388, 0.937339, 0.949528, 0.960431, 0.970031, 0.978317, 0.985278, 0.990903, 0.995185, 0.998118, 0.999699, 0.999925, 0.998795, 0.996313, 0.99248, 0.987301, 0.980785, 0.97294, 0.963776, 0.953306, 0.941544, 0.928506, 0.91421, 0.898674, 0.881921, 0.863973, 0.844854, 0.824589, 0.803208, 0.780737, 0.757209, 0.732654, 0.707107, 0.680601, 0.653173, 0.624859, 0.595699, 0.565732, 0.534998, 0.503538, 0.471397, 0.438616, 0.405241, 0.371317, 0.33689, 0.302006, 0.266713, 0.231058, 0.19509, 0.158858, 0.122411, 0.0857973, 0.0490677, 0.0122715, -0.0245412, -0.0613207, -0.0980171, -0.134581, -0.170962, -0.207111, -0.24298, -0.27852, -0.313682, -0.348419, -0.382683, -0.41643, -0.449611, -0.482184, -0.514103, -0.545325, -0.575808, -0.605511, -0.634393, -0.662416, -0.689541, -0.715731, -0.740951, -0.765167, -0.788346, -0.810457, -0.83147, -0.851355, -0.870087, -0.88764, -0.903989, -0.919114, -0.932993, -0.945607, -0.95694, -0.966976, -0.975702, -0.983105, -0.989177, -0.993907, -0.99729, -0.999322, 0, 0.00920375, 0.0184067, 0.0276081, 0.0368072, 0.0460032, 0.0551952, 0.0643826, 0.0735646, 0.0827403, 0.091909, 0.10107, 0.110222, 0.119365, 0.128498, 0.13762, 0.14673, 0.155828, 0.164913, 0.173984, 0.18304, 0.19208, 0.201105, 0.210112, 0.219101, 0.228072, 0.237024, 0.245955, 0.254866, 0.263755, 0.272621, 0.281465, 0.290285, 0.29908, 0.30785, 0.316593, 0.32531, 0.334, 0.342661, 0.351293, 0.359895, 0.368467, 0.377007, 0.385516, 0.393992, 0.402435, 0.410843, 0.419217, 0.427555, 0.435857, 0.444122, 0.45235, 0.460539, 0.468689, 0.476799, 0.484869, 0.492898, 0.500885, 0.50883, 0.516732, 0.52459, 0.532403, 0.540171, 0.547894, 0.55557, 0.563199, 0.570781, 0.578314, 0.585798, 0.593232, 0.600616, 0.60795, 0.615232, 0.622461, 0.629638, 0.636762, 0.643832, 0.650847, 0.657807, 0.664711, 0.671559, 0.67835, 0.685084, 0.691759, 0.698376, 0.704934, 0.711432, 0.71787, 0.724247, 0.730563, 0.736817, 0.743008, 0.749136, 0.755201, 0.761202, 0.767139, 0.77301, 0.778817, 0.784557, 0.79023, 0.795837, 0.801376, 0.806848, 0.812251, 0.817585, 0.82285, 0.828045, 0.83317, 0.838225, 0.843208, 0.84812, 0.852961, 0.857729, 0.862424, 0.867046, 0.871595, 0.87607, 0.880471, 0.884797, 0.889048, 0.893224, 0.897325, 0.901349, 0.905297, 0.909168, 0.912962, 0.916679, 0.920318, 0.92388, 0.927363, 0.930767, 0.934093, 0.937339, 0.940506, 0.943593, 0.946601, 0.949528, 0.952375, 0.955141, 0.957826, 0.960431, 0.962953, 0.965394, 0.967754, 0.970031, 0.972226, 0.974339, 0.97637, 0.978317, 0.980182, 0.981964, 0.983662, 0.985278, 0.986809, 0.988258, 0.989622, 0.990903, 0.992099, 0.993212, 0.99424, 0.995185, 0.996045, 0.99682, 0.997511, 0.998118, 0.99864, 0.999078, 0.999431, 0.999699, 0.999882, 0.999981, 0.999995, 0.999925, 0.999769, 0.999529, 0.999205, 0.998795, 0.998302, 0.997723, 0.99706, 0.996313, 0.995481, 0.994565, 0.993564, 0.99248, 0.991311, 0.990058, 0.988722, 0.987301, 0.985798, 0.98421, 0.982539, 0.980785, 0.978948, 0.977028, 0.975025, 0.97294, 0.970772, 0.968522, 0.96619, 0.963776, 0.96128, 0.958703, 0.956045, 0.953306, 0.950486, 0.947586, 0.944605, 0.941544, 0.938404, 0.935184, 0.931884, 0.928506, 0.925049, 0.921514, 0.917901, 0.91421, 0.910441, 0.906596, 0.902673, 0.898674, 0.894599, 0.890449, 0.886223, 0.881921, 0.877545, 0.873095, 0.868571, 0.863973, 0.859302, 0.854558, 0.849742, 0.844854, 0.839894, 0.834863, 0.829761, 0.824589, 0.819348, 0.814036, 0.808656, 0.803208, 0.797691, 0.792107, 0.786455, 0.780737, 0.774953, 0.769103, 0.763188, 0.757209, 0.751165, 0.745058, 0.738887, 0.732654, 0.726359, 0.720003, 0.713585, 0.707107, 0.700569, 0.693971, 0.687315, 0.680601, 0.673829, 0.667, 0.660114, 0.653173, 0.646176, 0.639124, 0.632019, 0.624859, 0.617647, 0.610383, 0.603067, 0.595699, 0.588282, 0.580814, 0.573297, 0.565732, 0.558119, 0.550458, 0.542751, 0.534998, 0.527199, 0.519356, 0.511469, 0.503538, 0.495565, 0.48755, 0.479494, 0.471397, 0.46326, 0.455084, 0.446869, 0.438616, 0.430326, 0.422, 0.413638, 0.405241, 0.39681, 0.388345, 0.379847, 0.371317, 0.362756, 0.354164, 0.345541, 0.33689, 0.32821, 0.319502, 0.310767, 0.302006, 0.293219, 0.284408, 0.275572, 0.266713, 0.257831, 0.248928, 0.240003, 0.231058, 0.222094, 0.21311, 0.204109, 0.19509, 0.186055, 0.177004, 0.167938, 0.158858, 0.149765, 0.140658, 0.13154, 0.122411, 0.113271, 0.104122, 0.0949635, 0.0857973, 0.0766239, 0.0674439, 0.0582583, 0.0490677, 0.0398729, 0.0306748, 0.0214741, 0.0122715, 0.00306796, -0.00613588, -0.0153392, -0.0245412, -0.0337412, -0.0429383, -0.0521317, -0.0613207, -0.0705046, -0.0796824, -0.0888536, -0.0980171, -0.107172, -0.116319, -0.125455, -0.134581, -0.143695, -0.152797, -0.161886, -0.170962, -0.180023, -0.189069, -0.198098, -0.207111, -0.216107, -0.225084, -0.234042, -0.24298, -0.251898, -0.260794, -0.269668, -0.27852, -0.287347, -0.296151, -0.304929, -0.313682, -0.322408, -0.331106, -0.339777, -0.348419, -0.357031, -0.365613, -0.374164, -0.382683, -0.39117, -0.399624, -0.408044, -0.41643, -0.42478, -0.433094, -0.441371, -0.449611, -0.457813, -0.465976, -0.4741, -0.482184, -0.490226, -0.498228, -0.506187, -0.514103, -0.521975, -0.529804, -0.537587, -0.545325, -0.553017, -0.560662, -0.568259, -0.575808, -0.583309, -0.59076, -0.598161, -0.605511, -0.61281, -0.620057, -0.627252, -0.634393, -0.641481, -0.648514, -0.655493, -0.662416, -0.669283, -0.676093, -0.682846, -0.689541, -0.696177, -0.702755, -0.709273, -0.715731, -0.722128, -0.728464, -0.734739, -0.740951, -0.747101, -0.753187, -0.759209, -0.765167, -0.771061, -0.776888, -0.782651, -0.788346, -0.793975, -0.799537, -0.805031, -0.810457, -0.815814, -0.821103, -0.826321, -0.83147, -0.836548, -0.841555, -0.846491, -0.851355, -0.856147, -0.860867, -0.865514, -0.870087, -0.874587, -0.879012, -0.883363, -0.88764, -0.891841, -0.895966, -0.900016, -0.903989, -0.907886, -0.911706, -0.915449, -0.919114, -0.922701, -0.92621, -0.929641, -0.932993, -0.936266, -0.939459, -0.942573, -0.945607, -0.948561, -0.951435, -0.954228, -0.95694, -0.959572, -0.962121, -0.96459, -0.966976, -0.969281, -0.971504, -0.973644, -0.975702, -0.977677, -0.97957, -0.981379, -0.983105, -0.984749, -0.986308, -0.987784, -0.989177, -0.990485, -0.99171, -0.99285, -0.993907, -0.994879, -0.995767, -0.996571, -0.99729, -0.997925, -0.998476, -0.998941, -0.999322, -0.999619, -0.999831, -0.999958};
+
+void
+fft_dit4_core_p1(double *fr, double *fi, ulong ldn)
+// Auxiliary routine for fft_dit4()
+// Decimation in time (DIT) radix-4 FFT
+// Input data must be in revbin_permuted order
+// ldn := base-2 logarithm of the array length
+// Fixed isign = +1
+{
+    int tmp = 0;// for getting the precalculated cos and sin
+    
+    const ulong n = (1UL<<ldn);
+    
+    if ( n<=2 )
+    {
+        if ( n==2 )
+        {
+            sumdiff(fr[0], fr[1]);
+            sumdiff(fi[0], fi[1]);
+        }
+        return;
+    }
+    
+    
+    ulong ldm = (ldn&1);
+    if ( ldm!=0 )  // n is not a power of 4, need a radix-8 step
+    {
+        for (ulong i0=0; i0<n; i0+=8)  fft8_dit_core_p1(fr+i0, fi+i0);
+    }
+    else
+    {
+        for (ulong i0=0; i0<n; i0+=4)
+        {
+            double xr, yr, ur, vr, xi, yi, ui, vi;
+            ulong i1 = i0 + 1;
+            ulong i2 = i1 + 1;
+            ulong i3 = i2 + 1;
+            
+            sumdiff(fr[i0], fr[i1], xr, ur);
+            sumdiff(fr[i2], fr[i3], yr, vi);
+            sumdiff(fi[i0], fi[i1], xi, ui);
+            sumdiff(fi[i3], fi[i2], yi, vr);
+            
+            sumdiff(ui, vi, fi[i1], fi[i3]);
+            sumdiff(xi, yi, fi[i0], fi[i2]);
+            sumdiff(ur, vr, fr[i1], fr[i3]);
+            sumdiff(xr, yr, fr[i0], fr[i2]);
+        }
+    }
+    ldm += 2*LX;
+    
+    
+    for ( ; ldm<=ldn; ldm+=LX)
+    {
+        ulong m = (1UL<<ldm);
+        ulong m4 = (m>>LX);
+        
+        for (ulong j=0; j<m4; j++)
+        {
+            double c, s, c2, s2, c3, s3;
+            c = allTheC[tmp];
+            s = allTheS[tmp];
+            c2 = allTheC2[tmp];
+            s2 = allTheS2[tmp];
+            c3 = allTheC3[tmp];
+            s3 = allTheS3[tmp];
+            tmp = tmp + 1;
+            
+            for (ulong r=0; r<n; r+=m)
+            {
+                ulong i0 = j + r;
+                ulong i1 = i0 + m4;
+                ulong i2 = i1 + m4;
+                ulong i3 = i2 + m4;
+                
+                double xr, yr, ur, vr, xi, yi, ui, vi;
+                cmult(c2, s2, fr[i1], fi[i1], xr, xi);
+                
+                sumdiff3_r(xr, fr[i0], ur);
+                sumdiff3_r(xi, fi[i0], ui);
+                
+                cmult(c,  s,  fr[i2], fi[i2], yr, vr);
+                cmult(c3, s3, fr[i3], fi[i3], vi, yi);
+                
+                sumdiff(yr, vi);
+                sumdiff(yi, vr);
+                
+                sumdiff(ur, vr, fr[i1], fr[i3]);
+                sumdiff(ui, vi, fi[i1], fi[i3]);
+                sumdiff(xr, yr, fr[i0], fr[i2]);
+                sumdiff(xi, yi, fi[i0], fi[i2]);
+            }
+        }
+    }
+}
+// -------------------------
+//////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+//// store all the precalculated indexes in the array
+int PermutedIndex [2048] = {0, 1024, 512, 1536, 256, 1280, 768, 1792, 128, 1152, 640, 1664, 384, 1408, 896, 1920, 64, 1088, 576, 1600, 320, 1344, 832, 1856, 192, 1216, 704, 1728, 448, 1472, 960, 1984, 32, 1056, 544, 1568, 288, 1312, 800, 1824, 160, 1184, 672, 1696, 416, 1440, 928, 1952, 96, 1120, 608, 1632, 352, 1376, 864, 1888, 224, 1248, 736, 1760, 480, 1504, 992, 2016, 16, 1040, 528, 1552, 272, 1296, 784, 1808, 144, 1168, 656, 1680, 400, 1424, 912, 1936, 80, 1104, 592, 1616, 336, 1360, 848, 1872, 208, 1232, 720, 1744, 464, 1488, 976, 2000, 48, 1072, 560, 1584, 304, 1328, 816, 1840, 176, 1200, 688, 1712, 432, 1456, 944, 1968, 112, 1136, 624, 1648, 368, 1392, 880, 1904, 240, 1264, 752, 1776, 496, 1520, 1008, 2032, 8, 1032, 520, 1544, 264, 1288, 776, 1800, 136, 1160, 648, 1672, 392, 1416, 904, 1928, 72, 1096, 584, 1608, 328, 1352, 840, 1864, 200, 1224, 712, 1736, 456, 1480, 968, 1992, 40, 1064, 552, 1576, 296, 1320, 808, 1832, 168, 1192, 680, 1704, 424, 1448, 936, 1960, 104, 1128, 616, 1640, 360, 1384, 872, 1896, 232, 1256, 744, 1768, 488, 1512, 1000, 2024, 24, 1048, 536, 1560, 280, 1304, 792, 1816, 152, 1176, 664, 1688, 408, 1432, 920, 1944, 88, 1112, 600, 1624, 344, 1368, 856, 1880, 216, 1240, 728, 1752, 472, 1496, 984, 2008, 56, 1080, 568, 1592, 312, 1336, 824, 1848, 184, 1208, 696, 1720, 440, 1464, 952, 1976, 120, 1144, 632, 1656, 376, 1400, 888, 1912, 248, 1272, 760, 1784, 504, 1528, 1016, 2040, 4, 1028, 516, 1540, 260, 1284, 772, 1796, 132, 1156, 644, 1668, 388, 1412, 900, 1924, 68, 1092, 580, 1604, 324, 1348, 836, 1860, 196, 1220, 708, 1732, 452, 1476, 964, 1988, 36, 1060, 548, 1572, 292, 1316, 804, 1828, 164, 1188, 676, 1700, 420, 1444, 932, 1956, 100, 1124, 612, 1636, 356, 1380, 868, 1892, 228, 1252, 740, 1764, 484, 1508, 996, 2020, 20, 1044, 532, 1556, 276, 1300, 788, 1812, 148, 1172, 660, 1684, 404, 1428, 916, 1940, 84, 1108, 596, 1620, 340, 1364, 852, 1876, 212, 1236, 724, 1748, 468, 1492, 980, 2004, 52, 1076, 564, 1588, 308, 1332, 820, 1844, 180, 1204, 692, 1716, 436, 1460, 948, 1972, 116, 1140, 628, 1652, 372, 1396, 884, 1908, 244, 1268, 756, 1780, 500, 1524, 1012, 2036, 12, 1036, 524, 1548, 268, 1292, 780, 1804, 140, 1164, 652, 1676, 396, 1420, 908, 1932, 76, 1100, 588, 1612, 332, 1356, 844, 1868, 204, 1228, 716, 1740, 460, 1484, 972, 1996, 44, 1068, 556, 1580, 300, 1324, 812, 1836, 172, 1196, 684, 1708, 428, 1452, 940, 1964, 108, 1132, 620, 1644, 364, 1388, 876, 1900, 236, 1260, 748, 1772, 492, 1516, 1004, 2028, 28, 1052, 540, 1564, 284, 1308, 796, 1820, 156, 1180, 668, 1692, 412, 1436, 924, 1948, 92, 1116, 604, 1628, 348, 1372, 860, 1884, 220, 1244, 732, 1756, 476, 1500, 988, 2012, 60, 1084, 572, 1596, 316, 1340, 828, 1852, 188, 1212, 700, 1724, 444, 1468, 956, 1980, 124, 1148, 636, 1660, 380, 1404, 892, 1916, 252, 1276, 764, 1788, 508, 1532, 1020, 2044, 2, 1026, 514, 1538, 258, 1282, 770, 1794, 130, 1154, 642, 1666, 386, 1410, 898, 1922, 66, 1090, 578, 1602, 322, 1346, 834, 1858, 194, 1218, 706, 1730, 450, 1474, 962, 1986, 34, 1058, 546, 1570, 290, 1314, 802, 1826, 162, 1186, 674, 1698, 418, 1442, 930, 1954, 98, 1122, 610, 1634, 354, 1378, 866, 1890, 226, 1250, 738, 1762, 482, 1506, 994, 2018, 18, 1042, 530, 1554, 274, 1298, 786, 1810, 146, 1170, 658, 1682, 402, 1426, 914, 1938, 82, 1106, 594, 1618, 338, 1362, 850, 1874, 210, 1234, 722, 1746, 466, 1490, 978, 2002, 50, 1074, 562, 1586, 306, 1330, 818, 1842, 178, 1202, 690, 1714, 434, 1458, 946, 1970, 114, 1138, 626, 1650, 370, 1394, 882, 1906, 242, 1266, 754, 1778, 498, 1522, 1010, 2034, 10, 1034, 522, 1546, 266, 1290, 778, 1802, 138, 1162, 650, 1674, 394, 1418, 906, 1930, 74, 1098, 586, 1610, 330, 1354, 842, 1866, 202, 1226, 714, 1738, 458, 1482, 970, 1994, 42, 1066, 554, 1578, 298, 1322, 810, 1834, 170, 1194, 682, 1706, 426, 1450, 938, 1962, 106, 1130, 618, 1642, 362, 1386, 874, 1898, 234, 1258, 746, 1770, 490, 1514, 1002, 2026, 26, 1050, 538, 1562, 282, 1306, 794, 1818, 154, 1178, 666, 1690, 410, 1434, 922, 1946, 90, 1114, 602, 1626, 346, 1370, 858, 1882, 218, 1242, 730, 1754, 474, 1498, 986, 2010, 58, 1082, 570, 1594, 314, 1338, 826, 1850, 186, 1210, 698, 1722, 442, 1466, 954, 1978, 122, 1146, 634, 1658, 378, 1402, 890, 1914, 250, 1274, 762, 1786, 506, 1530, 1018, 2042, 6, 1030, 518, 1542, 262, 1286, 774, 1798, 134, 1158, 646, 1670, 390, 1414, 902, 1926, 70, 1094, 582, 1606, 326, 1350, 838, 1862, 198, 1222, 710, 1734, 454, 1478, 966, 1990, 38, 1062, 550, 1574, 294, 1318, 806, 1830, 166, 1190, 678, 1702, 422, 1446, 934, 1958, 102, 1126, 614, 1638, 358, 1382, 870, 1894, 230, 1254, 742, 1766, 486, 1510, 998, 2022, 22, 1046, 534, 1558, 278, 1302, 790, 1814, 150, 1174, 662, 1686, 406, 1430, 918, 1942, 86, 1110, 598, 1622, 342, 1366, 854, 1878, 214, 1238, 726, 1750, 470, 1494, 982, 2006, 54, 1078, 566, 1590, 310, 1334, 822, 1846, 182, 1206, 694, 1718, 438, 1462, 950, 1974, 118, 1142, 630, 1654, 374, 1398, 886, 1910, 246, 1270, 758, 1782, 502, 1526, 1014, 2038, 14, 1038, 526, 1550, 270, 1294, 782, 1806, 142, 1166, 654, 1678, 398, 1422, 910, 1934, 78, 1102, 590, 1614, 334, 1358, 846, 1870, 206, 1230, 718, 1742, 462, 1486, 974, 1998, 46, 1070, 558, 1582, 302, 1326, 814, 1838, 174, 1198, 686, 1710, 430, 1454, 942, 1966, 110, 1134, 622, 1646, 366, 1390, 878, 1902, 238, 1262, 750, 1774, 494, 1518, 1006, 2030, 30, 1054, 542, 1566, 286, 1310, 798, 1822, 158, 1182, 670, 1694, 414, 1438, 926, 1950, 94, 1118, 606, 1630, 350, 1374, 862, 1886, 222, 1246, 734, 1758, 478, 1502, 990, 2014, 62, 1086, 574, 1598, 318, 1342, 830, 1854, 190, 1214, 702, 1726, 446, 1470, 958, 1982, 126, 1150, 638, 1662, 382, 1406, 894, 1918, 254, 1278, 766, 1790, 510, 1534, 1022, 2046, 1, 1025, 513, 1537, 257, 1281, 769, 1793, 129, 1153, 641, 1665, 385, 1409, 897, 1921, 65, 1089, 577, 1601, 321, 1345, 833, 1857, 193, 1217, 705, 1729, 449, 1473, 961, 1985, 33, 1057, 545, 1569, 289, 1313, 801, 1825, 161, 1185, 673, 1697, 417, 1441, 929, 1953, 97, 1121, 609, 1633, 353, 1377, 865, 1889, 225, 1249, 737, 1761, 481, 1505, 993, 2017, 17, 1041, 529, 1553, 273, 1297, 785, 1809, 145, 1169, 657, 1681, 401, 1425, 913, 1937, 81, 1105, 593, 1617, 337, 1361, 849, 1873, 209, 1233, 721, 1745, 465, 1489, 977, 2001, 49, 1073, 561, 1585, 305, 1329, 817, 1841, 177, 1201, 689, 1713, 433, 1457, 945, 1969, 113, 1137, 625, 1649, 369, 1393, 881, 1905, 241, 1265, 753, 1777, 497, 1521, 1009, 2033, 9, 1033, 521, 1545, 265, 1289, 777, 1801, 137, 1161, 649, 1673, 393, 1417, 905, 1929, 73, 1097, 585, 1609, 329, 1353, 841, 1865, 201, 1225, 713, 1737, 457, 1481, 969, 1993, 41, 1065, 553, 1577, 297, 1321, 809, 1833, 169, 1193, 681, 1705, 425, 1449, 937, 1961, 105, 1129, 617, 1641, 361, 1385, 873, 1897, 233, 1257, 745, 1769, 489, 1513, 1001, 2025, 25, 1049, 537, 1561, 281, 1305, 793, 1817, 153, 1177, 665, 1689, 409, 1433, 921, 1945, 89, 1113, 601, 1625, 345, 1369, 857, 1881, 217, 1241, 729, 1753, 473, 1497, 985, 2009, 57, 1081, 569, 1593, 313, 1337, 825, 1849, 185, 1209, 697, 1721, 441, 1465, 953, 1977, 121, 1145, 633, 1657, 377, 1401, 889, 1913, 249, 1273, 761, 1785, 505, 1529, 1017, 2041, 5, 1029, 517, 1541, 261, 1285, 773, 1797, 133, 1157, 645, 1669, 389, 1413, 901, 1925, 69, 1093, 581, 1605, 325, 1349, 837, 1861, 197, 1221, 709, 1733, 453, 1477, 965, 1989, 37, 1061, 549, 1573, 293, 1317, 805, 1829, 165, 1189, 677, 1701, 421, 1445, 933, 1957, 101, 1125, 613, 1637, 357, 1381, 869, 1893, 229, 1253, 741, 1765, 485, 1509, 997, 2021, 21, 1045, 533, 1557, 277, 1301, 789, 1813, 149, 1173, 661, 1685, 405, 1429, 917, 1941, 85, 1109, 597, 1621, 341, 1365, 853, 1877, 213, 1237, 725, 1749, 469, 1493, 981, 2005, 53, 1077, 565, 1589, 309, 1333, 821, 1845, 181, 1205, 693, 1717, 437, 1461, 949, 1973, 117, 1141, 629, 1653, 373, 1397, 885, 1909, 245, 1269, 757, 1781, 501, 1525, 1013, 2037, 13, 1037, 525, 1549, 269, 1293, 781, 1805, 141, 1165, 653, 1677, 397, 1421, 909, 1933, 77, 1101, 589, 1613, 333, 1357, 845, 1869, 205, 1229, 717, 1741, 461, 1485, 973, 1997, 45, 1069, 557, 1581, 301, 1325, 813, 1837, 173, 1197, 685, 1709, 429, 1453, 941, 1965, 109, 1133, 621, 1645, 365, 1389, 877, 1901, 237, 1261, 749, 1773, 493, 1517, 1005, 2029, 29, 1053, 541, 1565, 285, 1309, 797, 1821, 157, 1181, 669, 1693, 413, 1437, 925, 1949, 93, 1117, 605, 1629, 349, 1373, 861, 1885, 221, 1245, 733, 1757, 477, 1501, 989, 2013, 61, 1085, 573, 1597, 317, 1341, 829, 1853, 189, 1213, 701, 1725, 445, 1469, 957, 1981, 125, 1149, 637, 1661, 381, 1405, 893, 1917, 253, 1277, 765, 1789, 509, 1533, 1021, 2045, 3, 1027, 515, 1539, 259, 1283, 771, 1795, 131, 1155, 643, 1667, 387, 1411, 899, 1923, 67, 1091, 579, 1603, 323, 1347, 835, 1859, 195, 1219, 707, 1731, 451, 1475, 963, 1987, 35, 1059, 547, 1571, 291, 1315, 803, 1827, 163, 1187, 675, 1699, 419, 1443, 931, 1955, 99, 1123, 611, 1635, 355, 1379, 867, 1891, 227, 1251, 739, 1763, 483, 1507, 995, 2019, 19, 1043, 531, 1555, 275, 1299, 787, 1811, 147, 1171, 659, 1683, 403, 1427, 915, 1939, 83, 1107, 595, 1619, 339, 1363, 851, 1875, 211, 1235, 723, 1747, 467, 1491, 979, 2003, 51, 1075, 563, 1587, 307, 1331, 819, 1843, 179, 1203, 691, 1715, 435, 1459, 947, 1971, 115, 1139, 627, 1651, 371, 1395, 883, 1907, 243, 1267, 755, 1779, 499, 1523, 1011, 2035, 11, 1035, 523, 1547, 267, 1291, 779, 1803, 139, 1163, 651, 1675, 395, 1419, 907, 1931, 75, 1099, 587, 1611, 331, 1355, 843, 1867, 203, 1227, 715, 1739, 459, 1483, 971, 1995, 43, 1067, 555, 1579, 299, 1323, 811, 1835, 171, 1195, 683, 1707, 427, 1451, 939, 1963, 107, 1131, 619, 1643, 363, 1387, 875, 1899, 235, 1259, 747, 1771, 491, 1515, 1003, 2027, 27, 1051, 539, 1563, 283, 1307, 795, 1819, 155, 1179, 667, 1691, 411, 1435, 923, 1947, 91, 1115, 603, 1627, 347, 1371, 859, 1883, 219, 1243, 731, 1755, 475, 1499, 987, 2011, 59, 1083, 571, 1595, 315, 1339, 827, 1851, 187, 1211, 699, 1723, 443, 1467, 955, 1979, 123, 1147, 635, 1659, 379, 1403, 891, 1915, 251, 1275, 763, 1787, 507, 1531, 1019, 2043, 7, 1031, 519, 1543, 263, 1287, 775, 1799, 135, 1159, 647, 1671, 391, 1415, 903, 1927, 71, 1095, 583, 1607, 327, 1351, 839, 1863, 199, 1223, 711, 1735, 455, 1479, 967, 1991, 39, 1063, 551, 1575, 295, 1319, 807, 1831, 167, 1191, 679, 1703, 423, 1447, 935, 1959, 103, 1127, 615, 1639, 359, 1383, 871, 1895, 231, 1255, 743, 1767, 487, 1511, 999, 2023, 23, 1047, 535, 1559, 279, 1303, 791, 1815, 151, 1175, 663, 1687, 407, 1431, 919, 1943, 87, 1111, 599, 1623, 343, 1367, 855, 1879, 215, 1239, 727, 1751, 471, 1495, 983, 2007, 55, 1079, 567, 1591, 311, 1335, 823, 1847, 183, 1207, 695, 1719, 439, 1463, 951, 1975, 119, 1143, 631, 1655, 375, 1399, 887, 1911, 247, 1271, 759, 1783, 503, 1527, 1015, 2039, 15, 1039, 527, 1551, 271, 1295, 783, 1807, 143, 1167, 655, 1679, 399, 1423, 911, 1935, 79, 1103, 591, 1615, 335, 1359, 847, 1871, 207, 1231, 719, 1743, 463, 1487, 975, 1999, 47, 1071, 559, 1583, 303, 1327, 815, 1839, 175, 1199, 687, 1711, 431, 1455, 943, 1967, 111, 1135, 623, 1647, 367, 1391, 879, 1903, 239, 1263, 751, 1775, 495, 1519, 1007, 2031, 31, 1055, 543, 1567, 287, 1311, 799, 1823, 159, 1183, 671, 1695, 415, 1439, 927, 1951, 95, 1119, 607, 1631, 351, 1375, 863, 1887, 223, 1247, 735, 1759, 479, 1503, 991, 2015, 63, 1087, 575, 1599, 319, 1343, 831, 1855, 191, 1215, 703, 1727, 447, 1471, 959, 1983, 127, 1151, 639, 1663, 383, 1407, 895, 1919, 255, 1279, 767, 1791, 511, 1535, 1023, 2047};
+
+ulong ldn = 11;
+
+
+Complex* CFFT::convolutionF(const Complex *input, const Complex *filter, long nSIG, long NFIL, long &NFFT)
+{
+    //Check for invalid inputs.
+    if(input == NULL || filter == NULL)
+    {
+        cout << "Could not perform convolution on empty arrays!" << endl;
+        return NULL;
+    }
+    
+    bool NFFTChanged = false;
+    //If NFFT not a power of 2, or it is smaller than signal or filter, prompt for new.
+    while (log(NFFT) / log(2) != (int)(log(NFFT) / log(2)) || NFFT < nSIG || NFFT < NFIL)
+    {
+        cout << "Please input a valid NFFT, which is >= nSIG(" << nSIG << ") and >= NFIL(" << NFIL <<") : ";
+        cin >> NFFT;
+        NFFTChanged = true;
+    }
+    
+    //Perform zero padding.
+    Complex *fInput, *fFilter;
+    
+    fInput = new Complex[NFFT];
+    for(int i = 0; i < nSIG; i++)
+        fInput[i] = input[i];
+    
+    fFilter = new Complex[NFFT];
+    for(int i = 0; i < NFIL; i++)
+        fFilter[i] = filter[i];
+    
+    //Store the output data.
+    Complex *output = new Complex[NFFT];
+    
+    //Perform FFT on both input and filter.
+    CFFT::Forward(fInput, (unsigned int)NFFT);
+    CFFT::Forward(fFilter, (unsigned int)NFFT);
+    
+    for(int i = 0; i < NFFT; i++)
+        output[i] = fInput[i] * fFilter[i];// SIMD?
+    
+    return output;
+}
+
+//Convolution function which returns the time representation of the result.
+//NFFT is the FFT size, nSIG is the size for the input, NFIL is the size of the filter.
+Complex* CFFT::convolutionT(const Complex *input,const Complex *filter, long nSIG, long NFIL, long &NFFT)
+{
+    //Store the output data.
+    Complex *output = convolutionF(input, filter, nSIG, NFIL, NFFT);
+    
+    //Perform IFFT on the ouput.
+    CFFT::Inverse(output, (unsigned int)NFFT);
+    
+    return output;
+}
+
+Complex* CFFT::stereoConvMonoInputF(const Complex *input,const Complex *filterLeft,const Complex *filterRight, long nSIG, long NFILL, long NFILR, long &NFFT)
+{
+    Complex *result = stereoConvMonoInputT(input, filterLeft, filterRight, nSIG, NFILL, NFILR, NFFT);
+    CFFT::Forward(result, (unsigned int)NFFT);
+    return result;
+}
+
+//The size of the ouput will be 2 times of the size of FFT for the input signal and the NFFT value
+//will be doubled after running the function
+Complex* CFFT::stereoConvMonoInputT(const Complex *input,const Complex *filterLeft,const Complex *filterRight, long nSIG, long NFILL, long NFilR, long &NFFT)
+{
+    Complex *tempLeft = new Complex[NFFT];
+    Complex *tempRight = new Complex[NFFT];
+    Complex *result = new Complex[2*NFFT];
+    
+    tempLeft = CFFT::convolutionT(input, filterLeft, nSIG, NFILL, NFFT);
+    tempRight = CFFT::convolutionT(input, filterRight, nSIG, NFILL, NFFT);
+    NFFT = NFFT * 2;
+    for (int i = 0; i < NFFT / 2; i++)
+    {
+        result[2 * i] = tempLeft[i];
+        result[2 * i + 1] = tempRight[i];
+    }
+    
+    delete tempLeft;
+    delete tempRight;
+    
+    return result;
+}
+
+/*
+ Complex* CFFT::stereoConvStereoInputT(const Complex *input, const Complex *filterLeft, const Complex *filterRight, long nSIG, long NFILL, long NFILR, long &NFFT)
+ {
+	NFFT = NFFT / 2;
+	Complex *leftTemp = new Complex[NFFT];
+	Complex *rightTemp = new Complex[NFFT];
+ 
+	
+ 
+	
+ 
+	return 0;
+ 
+ }
+ 
+ */
+
+
+//storing the an array into a text file
+//filename is the file name you want to store the data into
+//datatype represents the data you wanna store: real/real+imag/amplitude
+void CFFT::storingData(Complex *data, int NFFT,string temp ,char datatype)
+{
+    //string temp= filename;
+    ofstream outputFile(temp.c_str());
+    if (outputFile.is_open())
+    {
+        switch (datatype)
+        {
+            case 'r':
+                //ofstream outputFile("real.txt");
+                for (int i = 0; i < NFFT; i++)
+                    outputFile << data[i].real() << endl;
+                break;
+                
+            case 'c':
+                for (int i = 0; i < NFFT; i++)
+                    outputFile << data[i].real() << "    " << data[i].imag() << endl;
+                break;
+            case 'a':
+                for (int i = 0; i < NFFT; i++)
+                    outputFile << sqrt(data[i].real()*data[i].real()+data[i].imag()*data[i].imag()) << endl;
+                break;
+        }
+        outputFile.close();
+    }
+}
+
+
+
+//   FORWARD FOURIER TRANSFORM
+//     Input  - input data
+//     Output - transform result
+//     N      - length of both input data and result
+bool CFFT::Forward(const Complex *const Input, Complex *const Output, const unsigned int N)
+{
+    //   Check input parameters
+    if (!Input || !Output || N < 1 || N & (N - 1))
+        return false;
+    //   Initialize data
+    Rearrange(Input, Output, N);
+    //   Call FFT implementation
+    Perform(Output, N);
+    //   Succeeded
+    return true;
+}
+
+//   FORWARD FOURIER TRANSFORM, INPLACE VERSION
+//     Data - both input data and output
+//     N    - length of input data
+bool CFFT::Forward(Complex *const Data, const unsigned int N)
+{
+    //   Check input parameters
+    if (!Data || N < 1 || N & (N - 1))
+        return false;
+    //   Rearrange
+    Rearrange(Data, N);
+    //   Call FFT implementation
+    Perform(Data, N);
+    //   Succeeded
+    return true;
+}
+
+//   INVERSE FOURIER TRANSFORM
+//     Input  - input data
+//     Output - transform result
+//     N      - length of both input data and result
+//     Scale  - if to scale result
+bool CFFT::Inverse(const Complex *const Input, Complex *const Output, const unsigned int N, const bool Scale /* = true */)
+{
+    //   Check input parameters
+    if (!Input || !Output || N < 1 || N & (N - 1))
+        return false;
+    //   Initialize data
+    Rearrange(Input, Output, N);
+    //   Call FFT implementation
+    Perform(Output, N, true);
+    //   Scale if necessary
+    if (Scale)
+        CFFT::Scale(Output, N);
+    //   Succeeded
+    return true;
+}
+
+//   INVERSE FOURIER TRANSFORM, INPLACE VERSION
+//     Data  - both input data and output
+//     N     - length of both input data and result
+//     Scale - if to scale result
+bool CFFT::Inverse(Complex *const Data, const unsigned int N, const bool Scale /* = true */)
+{
+    //   Check input parameters
+    if (!Data || N < 1 || N & (N - 1))
+        return false;
+    //   Rearrange
+    Rearrange(Data, N);
+    //   Call FFT implementation
+    Perform(Data, N, true);
+    //   Scale if necessary
+    if (Scale)
+        CFFT::Scale(Data, N);
+    //   Succeeded
+    return true;
+}
+
+//   Rearrange function
+void CFFT::Rearrange(const Complex *const Input, Complex *const Output, const unsigned int N)
+{
+    //	//   Data entry position
+    //	unsigned int Target = 0;
+    //	//   Process all positions of input signal
+    for (unsigned int Position = 0; Position < N; ++Position)
+    {
+        Output[Position] = Input[Position];
+    }// SIMD?
+}
+
+//   Inplace version of rearrange function
+void CFFT::Rearrange(Complex *const Data, const unsigned int N){}
+
+
+//   FFT implementation
+void CFFT::Perform(Complex *const Data, const unsigned int N, const bool Inverse /* = false */)
+{
+
+    if (Inverse) {
+        // Inverse transform from Complex number to real number(Complex number with real part defined and imaginary part set to zero)  radix-4 decimation in Time method
+        double dataReal [N];
+        double dataImag [N];
+        
+        for (int i = 0; i < N; i++) {
+            //cout << "i: " << i << endl;
+            //cout << "getPermutedIndex(i): " << getPermutedIndex(i) << endl;
+            dataReal[i] = Data[PermutedIndex[i]].real();
+            dataImag[i] = Data[PermutedIndex[i]].imag();
+        }
+        double *Kr = dataReal;
+        double *Ki = dataImag;
+        
+        
+        //revbin_permute(Kr, 1UL<<ldn);
+        //revbin_permute(Ki, 1UL<<ldn);
+        
+        fft_dit4_core_p1(Kr, Ki, ldn);
+        
+        for (unsigned int Position = 0; Position < N; ++Position){
+            Data[Position] = std::complex <double> (Kr[Position], 0);
+        }
+        
+        
+    }
+    else{
+        // forward transform from real number to Complex number using FHT method
+        // first permute Data
+        double *f = new double [N];
+        for (int i = 0; i < N; i++) {
+            f[i] = Data[i].real();
+        }
+        revbin_permute(f, 1UL<<ldn);
+        // now tmpdata is the permuted version of Data
+        fht_loc_dit2_core(f, ldn);
+        /////////////  FHT to FFT conversion start ///////////////
+        const ulong n = (1UL<<ldn);
+        //int is = 1;
+        //if ( is>0 )  for (ulong i=1, j=n-1;  i<j;  i++, j--)  sumdiff05(f[i], f[j]);
+        //else         for (ulong i=1, j=n-1;  i<j;  i++, j--)  sumdiff05_r(f[i], f[j]);
+        for (ulong i=1, j=n-1; i<j;  i++, j--)  sumdiff05(f[i], f[j]);
+        
+        for (int i = 1; i<512; i++) {
+            Complex tmp(f[i],-f[N-i]);
+            Data[i] = tmp;
+        }
+        Data[512] = Complex(f[512],0);
+        for (int i = 513; i<N; i++) {
+            Complex tmp(f[N-i],f[i]);
+            Data[i] = tmp;
+        }
+        Data[0] = Complex(f[0],0);
+        delete f;
+        /////////////  FHT to FFT conversion finished ///////////////
+    }
+    
+}
+
+//   Scaling of inverse FFT result
+void CFFT::Scale(Complex *const Data, const unsigned int N)
+{
+    const double Factor = 1. / double(N);
+    //   Scale all data entries
+    for (unsigned int Position = 0; Position < N; ++Position)
+        Data[Position] *= Factor;
+}// SIMD?
diff -r 000000000000 -r 2076b4d80327 fhtdit.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/fhtdit.cc	Tue Apr 07 21:09:22 2015 +0000
@@ -0,0 +1,262 @@
+// This file is part of the FXT library.
+// Copyright (C) 2010, 2012 Joerg Arndt
+// License: GNU General Public License version 3 or later,
+// see the file COPYING.txt in the main directory.
+
+#include "fxttypes.h"
+#include "cmult.h"
+#include "sumdiff.h"
+#include "constants.h"  // COS_1_PI_8, SIN_1_PI_8
+#include "sincos.h"
+#include "revbinpermute.h"
+
+#include <cmath>  // M_PI, M_SQRT1_2
+
+// tuning parameter:
+// define to use trig recurrence:
+// (and possibly lose some precision, see below)
+//#define USE_TRIG_REC
+// with type 'long double' slight speed loss on my machine,
+// with type 'double' little speed gain.
+
+#if defined USE_TRIG_REC
+//#warning 'FYI: fht(double *, ulong) uses trig recursion'
+#endif
+
+// tuning parameter:
+#define  INITIAL_RADIX_16  1  // 0 or 1 (default)
+//
+#if  ( INITIAL_RADIX_16==1 )
+//#warning 'FYI: INITIAL_RADIX_16 set in fht_dit(double *, ulong)'
+#else
+#warning 'INITIAL_RADIX_16 is NOT SET in fht_dit(double *, ulong)'
+#endif
+
+
+void
+fht_dit_core(double *f, ulong ldn)
+// Decimation in time (DIT) FHT.
+// Input data must be in revbin_permuted order.
+// ldn := base-2 logarithm of the array length.
+{
+    if ( ldn<=2 )
+    {
+        if ( ldn==1 )  // two point fht
+        {
+            sumdiff(f[0], f[1]);
+        }
+        else if ( ldn==2 )  // four point fht
+        {
+            double f0, f1, f2, f3;
+            sumdiff(f[0], f[1], f0, f1);
+            sumdiff(f[2], f[3], f2, f3);
+            sumdiff(f0, f2, f[0], f[2]);
+            sumdiff(f1, f3, f[1], f[3]);
+        }
+        return;
+    }
+
+    const ulong n = (1UL<<ldn);
+    const double *fn = f + n;
+    ulong ldk = ldn & 1;
+
+    if ( ldk==0 )    // ldn is multiple of 2  => n is a power of 4
+    {
+#if  ( INITIAL_RADIX_16==1 )
+        for (double *fi=f; fi<fn; fi+=16)  // radix-16 step
+        {
+            double f0, f1, f2, f3;
+            sumdiff(fi[0], fi[1], f0, f1);
+            sumdiff(fi[2], fi[3], f2, f3);
+            sumdiff(f0, f2, fi[0], fi[2]);
+            sumdiff(f1, f3, fi[1], fi[3]);
+
+            sumdiff(fi[4], fi[5], f0, f1);
+            sumdiff(fi[6], fi[7], f2, f3);
+            sumdiff(f0, f2, fi[4], fi[6]);
+            sumdiff(f1, f3, fi[5], fi[7]);
+
+            sumdiff(fi[8], fi[9], f0, f1);
+            sumdiff(fi[10], fi[11], f2, f3);
+            sumdiff(f0, f2, fi[8], fi[10]);
+            sumdiff(f1, f3, fi[9], fi[11]);
+
+            sumdiff(fi[12], fi[13], f0, f1);
+            sumdiff(fi[14], fi[15], f2, f3);
+            sumdiff(f0, f2, fi[12], fi[14]);
+            sumdiff(f1, f3, fi[13], fi[15]);
+
+            sumdiff(fi[0], fi[4], f0, f1);
+            sumdiff(fi[8], fi[12], f2, f3);
+            sumdiff(f0, f2, fi[0], fi[8]);
+            sumdiff(f1, f3, fi[4], fi[12]);
+
+            sumdiff(fi[2], fi[6], f0, f1);
+            f3 = M_SQRT2 * fi[14];
+            f2 = M_SQRT2 * fi[10];
+            sumdiff(f0, f2, fi[2], fi[10]);
+            sumdiff(f1, f3, fi[6], fi[14]);
+
+            double a, b, g0, g1, g2, g3;
+            sumdiff(fi[5], fi[7], a, b);
+            a *= M_SQRT1_2;
+            b *= M_SQRT1_2;
+            sumdiff(fi[1], a, f0, f1);
+            sumdiff(fi[3], b, g0, g1);
+            sumdiff(fi[13], fi[15], a, b);
+            a *= M_SQRT1_2;
+            b *= M_SQRT1_2;
+            sumdiff(fi[9], a, f2, f3);
+            sumdiff(fi[11], b, g2, g3);
+            double c1 = COS_1_PI_8;  // jjkeep
+            double s1 = SIN_1_PI_8;  // jjkeep
+            cmult(s1, c1, f2, g3, b, a);
+            sumdiff(f0, a, fi[1], fi[9]);
+            sumdiff(g1, b, fi[7], fi[15]);
+            cmult(c1, s1, g2, f3, b, a);
+            sumdiff(g0, a, fi[3], fi[11]);
+            sumdiff(f1, b, fi[5], fi[13]);
+        }
+        ldk = 4;
+#else  // INITIAL_RADIX_16
+        for (double *fi=f; fi<fn; fi+=4)  // radix-4 step
+        {
+            double f0, f1, f2, f3;
+            sumdiff(fi[0], fi[1], f0, f1);
+            sumdiff(fi[2], fi[3], f2, f3);
+            sumdiff(f0, f2, fi[0], fi[2]);
+            sumdiff(f1, f3, fi[1], fi[3]);
+        }
+        ldk = 2;
+#endif  // INITIAL_RADIX_16
+    }
+    else          // ldk==1,  n is no power of 4
+    {
+        for (double *fi=f; fi<fn; fi+=8)  // radix-8 step
+        {
+            double g0, f0, f1, g1;
+            sumdiff(fi[0], fi[1], f0, g0);
+            sumdiff(fi[2], fi[3], f1, g1);
+
+            sumdiff(f0, f1);
+            sumdiff(g0, g1);
+
+            double s1, c1, s2, c2;
+            sumdiff(fi[4], fi[5], s1, c1);
+            sumdiff(fi[6], fi[7], s2, c2);
+
+            sumdiff(s1, s2);
+
+            sumdiff(f0, s1, fi[0], fi[4]);
+            sumdiff(f1, s2, fi[2], fi[6]);
+
+            c1 *= M_SQRT2;
+            c2 *= M_SQRT2;
+
+            sumdiff(g0, c1, fi[1], fi[5]);
+            sumdiff(g1, c2, fi[3], fi[7]);
+        }
+        ldk = 3;
+    }
+
+
+    for (  ; ldk<ldn;  ldk+=2)
+    {
+        ulong k   = 1UL<<ldk;
+        ulong kh  = k >> 1;
+        ulong k2  = k << 1;
+        ulong k3  = k2 + k;
+        ulong k4  = k2 << 1;
+
+        for (double *fi=f, *gi=f+kh;  fi<fn;  fi+=k4, gi+=k4)
+        {
+            double f0, f1, f2, f3;
+            sumdiff(fi[0], fi[k], f0, f1);
+            sumdiff(fi[k2], fi[k3], f2, f3);
+            sumdiff(f0, f2, fi[0], fi[k2]);
+            sumdiff(f1, f3, fi[k], fi[k3]);
+
+            sumdiff(gi[0], gi[k], f0, f1);
+            f3 = M_SQRT2 * gi[k3];
+            f2 = M_SQRT2 * gi[k2];
+            sumdiff(f0, f2, gi[0], gi[k2]);
+            sumdiff(f1, f3, gi[k], gi[k3]);
+        }
+
+        double tt = M_PI_4/(double)kh;  // jjkeep
+#if defined USE_TRIG_REC
+        double s1 = 0.0,  c1 = 1.0;  // jjkeep
+        double al = sin(0.5*tt);  // jjkeep
+        al *= (2.0*al);
+        double be = sin(tt);  // jjkeep
+#endif  // USE_TRIG_REC
+
+        for (ulong i=1; i<kh; i++)
+        {
+#if defined USE_TRIG_REC
+            {
+                double t1 = c1;  // jjkeep
+                c1 -= (al*t1+be*s1);
+                s1 -= (al*s1-be*t1);
+            }
+#else
+            double s1, c1;  // jjkeep
+            SinCos(tt*(double)i, &s1, &c1);  // jjkeep
+#endif  // USE_TRIG_REC
+
+            double c2, s2;  // jjkeep
+            csqr(c1, s1, c2, s2);
+
+            for (double *fi=f+i, *gi=f+k-i;  fi<fn;  fi+=k4, gi+=k4)
+            {
+                double a, b, g0, f0, f1, g1, f2, g2, f3, g3;
+                cmult(s2, c2, fi[k], gi[k], b, a);
+                sumdiff(fi[0], a, f0, f1);
+                sumdiff(gi[0], b, g0, g1);
+
+                cmult(s2, c2, fi[k3], gi[k3], b, a);
+                sumdiff(fi[k2], a, f2, f3);
+                sumdiff(gi[k2], b, g2, g3);
+
+                cmult(s1, c1, f2, g3, b, a);
+                sumdiff(f0, a, fi[0], fi[k2]);
+                sumdiff(g1, b, gi[k], gi[k3]);
+
+                cmult(c1, s1, g2, f3, b, a);
+                sumdiff(g0, a, gi[0], gi[k2]);
+                sumdiff(f1, b, fi[k], fi[k3]);
+            }
+        }
+    }
+}
+// -------------------------
+
+
+void
+fht_dit(double *f, ulong ldn)
+// Fast Hartley Transform.
+// Split-radix decimation in time (DIT) algorithm.
+// ldn := base-2 logarithm of the array length.
+{
+    if ( ldn<=2 )
+    {
+        if ( ldn==1 )  // two point fht
+        {
+            sumdiff(f[0], f[1]);
+        }
+        else if ( ldn==2 )  // four point fht
+        {
+            double f0, f1, f2, f3;
+            sumdiff(f[0], f[2], f0, f1);
+            sumdiff(f[1], f[3], f2, f3);
+            sumdiff(f0, f2, f[0], f[2]);
+            sumdiff(f1, f3, f[1], f[3]);
+        }
+
+        return;
+    }
+
+    revbin_permute(f, 1UL<<ldn);
+    fht_dit_core(f, ldn);
+}
+// -------------------------
diff -r 000000000000 -r 2076b4d80327 fhtdit2.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/fhtdit2.cc	Tue Apr 07 21:09:22 2015 +0000
@@ -0,0 +1,130 @@
+// This file is part of the FXT library.
+// Copyright (C) 2010, 2012 Joerg Arndt
+// License: GNU General Public License version 3 or later,
+// see the file COPYING.txt in the main directory.
+
+#include "fxttypes.h"
+#include "complextype.h"
+#include "sincos.h"
+#include "revbinpermute.h"
+#include "sumdiff.h"
+
+#include <cmath>  // M_PI
+
+
+void
+fht_depth_first_dit2(double *f, ulong ldn)
+// Radix-2 decimation in time (DIT) FHT.
+// Depth-first version.
+// Compared to usual fht this one
+// - does more trig computations
+// - is (far) better memory local
+{
+    const ulong n = 1UL<<ldn;
+    revbin_permute(f, n);
+
+    for (ulong ldm=1; ldm<=ldn; ++ldm)
+    {
+        const ulong m = (1UL<<ldm);
+        const ulong mh = (m>>1);
+        const ulong m4 = (mh>>1);
+        const double phi0 = M_PI/(double)mh;
+
+        for (ulong r=0; r<n; r+=m)
+        {
+            {  // j == 0:
+                ulong t1 = r;
+                ulong t2 = t1 + mh;
+                sumdiff(f[t1], f[t2]);
+            }
+
+            if ( m4 )
+            {
+                ulong t1 = r + m4;
+                ulong t2 = t1 + mh;
+                sumdiff(f[t1], f[t2]);
+            }
+
+            for (ulong j=1, k=mh-1;  j<k;  ++j, --k)
+            {
+                double s, c;
+                SinCos(phi0*(double)j, &s, &c);
+
+                ulong tj = r + mh + j;
+                ulong tk = r + mh + k;
+                double fj = f[tj];
+                double fk = f[tk];
+                f[tj] = fj * c + fk * s;
+                f[tk] = fj * s - fk * c;
+
+                ulong t1 = r + j;
+                ulong t2 = tj;  // == t1 + mh;
+                sumdiff(f[t1], f[t2]);
+
+                t1 = r + k;
+                t2 = tk;  // == t1 + mh;
+                sumdiff(f[t1], f[t2]);
+            }
+        }
+    }
+}
+// -------------------------
+
+
+void
+fht_dit2(double *f, ulong ldn)
+// Radix-2 decimation in time (DIT) FHT.
+{
+    const ulong n = 1UL<<ldn;
+    revbin_permute(f, n);
+
+    for (ulong ldm=1; ldm<=ldn; ++ldm)
+    {
+        const ulong m = (1UL<<ldm);
+        const ulong mh = (m>>1);
+        const ulong m4 = (mh>>1);
+        const double phi0 = M_PI/(double)mh;
+
+        for (ulong r=0; r<n; r+=m)
+        {
+            {  // j == 0:
+                ulong t1 = r;
+                ulong t2 = t1 + mh;
+                sumdiff(f[t1], f[t2]);
+            }
+
+            if ( m4 )
+            {
+                ulong t1 = r + m4;
+                ulong t2 = t1 + mh;
+                sumdiff(f[t1], f[t2]);
+            }
+        }
+
+        for (ulong j=1, k=mh-1;  j<k;  ++j, --k)
+        {
+            double s, c;
+            SinCos(phi0*(double)j, &s, &c);
+
+            for (ulong r=0; r<n; r+=m)
+            {
+                ulong tj = r + mh + j;
+                ulong tk = r + mh + k;
+                double fj = f[tj];
+                double fk = f[tk];
+                f[tj] = fj * c + fk * s;
+                f[tk] = fj * s - fk * c;
+
+                ulong t1 = r + j;
+                ulong t2 = tj;  // == t1 + mh;
+                sumdiff(f[t1], f[t2]);
+
+                t1 = r + k;
+                t2 = tk;  // == t1 + mh;
+                sumdiff(f[t1], f[t2]);
+            }
+        }
+    }
+}
+// -------------------------
+
diff -r 000000000000 -r 2076b4d80327 location.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/location.cpp	Tue Apr 07 21:09:22 2015 +0000
@@ -0,0 +1,54 @@
+#include "../include/location.h"
+#include <iostream>
+
+Location& Location::operator= (const Location& rhs)
+{
+    x = (rhs.x);
+    y = (rhs.y);
+    z = (rhs.z);
+    return *this;
+}
+
+Location& Location::operator+= (const Location& rhs)
+{
+    x = x + rhs.x;
+    y = y + rhs.y;
+    z = z + rhs.z;
+    return *this;
+}
+
+Location operator+ (Location lhs, const Location& rhs)
+{
+    lhs += rhs;
+    return lhs;
+}
+
+bool Location::operator< (const Location& rhs) const
+{
+    if (x == rhs.x){
+        if (y == rhs.y){
+            return z < rhs.z;
+        }return y < rhs.y;
+    }return x < rhs.x;
+}
+
+ostream& operator<< (ostream& os, Location obj)
+{
+    os << obj.getX() <<","<< obj.getY() <<","<< obj.getZ();
+    return os;
+}
+
+float Location::getX (void) const
+{
+    return this->x;
+}
+
+float Location::getY (void) const
+{
+    return this->y;
+}
+
+float Location::getZ (void) const
+{
+    return this->z;
+}
diff -r 000000000000 -r 2076b4d80327 mit_hrtf_lib.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mit_hrtf_lib.c	Tue Apr 07 21:09:22 2015 +0000
@@ -0,0 +1,467 @@
+/*############################################################################*/
+/*#                                                                          #*/
+/*#  MIT HRTF C Library                                                      #*/
+/*#  Copyright � 2007 Aristotel Digenis                                      #*/
+/*#                                                                          #*/
+/*#  Filename:  mit_hrtf_lib.c                                               #*/
+/*#  Version:   0.1                                                          #*/
+/*#  Date:      04/05/2007                                                   #*/
+/*#  Author(s): Aristotel Digenis (adigenis@users.sourceforge.net)           #*/
+/*#  Credit:    Bill Gardner and Keith Martin                                #*/
+/*#  Licence:   GNU Library or Lesser General Public License (LGPL)          #*/
+/*#                                                                          #*/
+/*############################################################################*/
+
+
+#include "../include/mit_hrtf_lib.h"
+#include "../MIT_HRTF_Library/source/normal/mit_hrtf_normal_44100.h"
+#include "../MIT_HRTF_Library/source/diffuse/mit_hrtf_diffuse_44100.h"
+#include "../MIT_HRTF_Library/source/normal/mit_hrtf_normal_48000.h"
+#include "../MIT_HRTF_Library/source/diffuse/mit_hrtf_diffuse_48000.h"
+#include "../MIT_HRTF_Library/source/normal/mit_hrtf_normal_88200.h"
+#include "../MIT_HRTF_Library/source/diffuse/mit_hrtf_diffuse_88200.h"
+#include "../MIT_HRTF_Library/source/normal/mit_hrtf_normal_96000.h"
+#include "../MIT_HRTF_Library/source/diffuse/mit_hrtf_diffuse_96000.h"
+
+
+/*	Internal functions for handling the indexing of the -/+40 degree elevation 
+	data which has irregular azimuth positions. */
+int mit_hrtf_findAzimuthFor40Elev(int azimuth);
+int mit_hrtf_findIndexFor40Elev(int azimuth);
+
+
+
+unsigned int mit_hrtf_availability(int azimuth, int elevation, unsigned int samplerate, unsigned int diffused)
+{
+	if(elevation > 90 || elevation < -40)
+		return 0;
+
+	if(azimuth > 180 || azimuth < -180)
+		return 0;
+    
+    if(diffused > 1)
+        return 0;
+
+	if(samplerate == 44100)
+		return MIT_HRTF_44_TAPS;
+	else if(samplerate == 48000)
+		return MIT_HRTF_48_TAPS;
+	else if(samplerate == 88200)
+		return MIT_HRTF_88_TAPS;
+	else if(samplerate == 96000)
+		return MIT_HRTF_96_TAPS;
+
+	return 0;
+}
+
+
+
+unsigned int mit_hrtf_get(int* pAzimuth, int* pElevation, unsigned int samplerate, unsigned int diffused, short* psLeft, short* psRight)
+{
+	int nInternalElevation = 0;
+	float fAzimuthIncrement = 0;
+	int nInternalAzimuth = 0;
+	int nSwitchLeftRight = 0;
+	int nAzimuthIndex = 0;
+	const mit_hrtf_filter_set_44* pFilter44 = 0;
+	const mit_hrtf_filter_set_48* pFilter48 = 0;
+	const mit_hrtf_filter_set_88* pFilter88 = 0;
+	const mit_hrtf_filter_set_96* pFilter96 = 0;
+	const short* psLeftTaps = 0;
+	const short* psRightTaps = 0;
+	const short* psTempTaps = 0;
+	unsigned int nTotalTaps = 0;
+	unsigned int niTap = 0;
+
+	//Check if the requested HRTF exists
+	if(!mit_hrtf_availability(*pAzimuth, *pElevation, samplerate, diffused))
+		return 0;
+
+	//Snap elevation to the nearest available elevation in the filter set
+	if(*pElevation < 0)
+		nInternalElevation = ((*pElevation - 5) / 10) * 10;
+	else
+		nInternalElevation = ((*pElevation + 5) / 10) * 10;
+
+	// Elevation of 50 has a maximum 176 in the azimuth plane so we need to handle that.
+	if(nInternalElevation == 50)
+	{
+		if(*pAzimuth < 0)
+			*pAzimuth = *pAzimuth < -176 ? -176 : *pAzimuth;
+		else
+			*pAzimuth = *pAzimuth > 176 ? 176 : *pAzimuth;
+	}
+
+	//Snap azimuth to the nearest available azimuth in the filter set.
+	switch(nInternalElevation)
+	{
+	case 0:		fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_00 - 1);		break;	// 180 5
+	case 10:	
+	case -10:	fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_10 - 1);		break;	// 180 5
+	case 20:	
+	case -20:	fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_20 - 1);		break;	// 180 5
+	case 30:	
+	case -30:	fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_30 - 1);		break;	// 180 6
+	case 40:	
+	case -40:	fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_40 - 1);		break;	// 180 6.43
+	case 50:	fAzimuthIncrement = 176.f / (MIT_HRTF_AZI_POSITIONS_50 - 1);		break;	// 176 8
+	case 60:	fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_60 - 1);		break;	// 180 10
+	case 70:	fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_70 - 1);		break;	// 180 15
+	case 80:	fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_80 - 1);		break;	// 180 30
+	case 90:	fAzimuthIncrement = 0;												break;	// 0   1
+	};
+
+	if(*pAzimuth < 0)
+	{
+		nInternalAzimuth = (int)((int)((-*pAzimuth + fAzimuthIncrement / 2.f) / fAzimuthIncrement) * fAzimuthIncrement + 0.5f);
+		nSwitchLeftRight = 1;
+	}
+	else
+	{
+		nInternalAzimuth = (int)((int)((*pAzimuth + fAzimuthIncrement / 2.f) / fAzimuthIncrement) * fAzimuthIncrement + 0.5f);
+	}
+
+	//Determine array index for azimuth based on elevation
+	switch(nInternalElevation)
+	{
+	case 0:		nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_00 - 1));	break;
+	case 10: 
+	case -10:	nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_10 - 1));	break;
+	case 20:	
+	case -20:	nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_20 - 1));	break;
+	case 30:	
+	case -30:	nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_30 - 1));	break;
+	case 40:	
+	case -40:	nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_40 - 1));	break;
+	case 50:	nAzimuthIndex = (int)((nInternalAzimuth / 176.f) * (MIT_HRTF_AZI_POSITIONS_50 - 1));	break;
+	case 60:	nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_60 - 1));	break;
+	case 70:	nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_70 - 1));	break;
+	case 80:	nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_80 - 1));	break;
+	case 90:	nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_90 - 1));	break;
+	};
+
+	//The azimuths for +/- elevations need special handling
+	if(nInternalElevation == 40 || nInternalElevation == -40)
+	{
+		nInternalAzimuth = mit_hrtf_findAzimuthFor40Elev(nInternalAzimuth);
+		nAzimuthIndex = mit_hrtf_findIndexFor40Elev(nInternalAzimuth);
+	}
+
+	//Assign pointer to appropriate array depending on saple rate, normal or diffused filters, elevation, and azimuth index.
+	switch(samplerate)
+	{
+	case 44100:
+		if(diffused)
+			pFilter44 = &diffuse_44;
+		else
+			pFilter44 = &normal_44;
+
+		switch(nInternalElevation)
+		{	
+		case -10:	psLeftTaps = pFilter44->e_10[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e_10[nAzimuthIndex].right;		break;
+		case -20:	psLeftTaps = pFilter44->e_20[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e_20[nAzimuthIndex].right;		break;
+		case -30:	psLeftTaps = pFilter44->e_30[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e_30[nAzimuthIndex].right;		break;
+		case -40:	psLeftTaps = pFilter44->e_40[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e_40[nAzimuthIndex].right;		break;
+		case 0:		psLeftTaps = pFilter44->e00[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e00[nAzimuthIndex].right;		break;
+		case 10:	psLeftTaps = pFilter44->e10[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e10[nAzimuthIndex].right;		break;
+		case 20:	psLeftTaps = pFilter44->e20[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e20[nAzimuthIndex].right;		break;
+		case 30:	psLeftTaps = pFilter44->e30[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e30[nAzimuthIndex].right;		break;
+		case 40:	psLeftTaps = pFilter44->e40[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e40[nAzimuthIndex].right;		break;
+		case 50:	psLeftTaps = pFilter44->e50[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e50[nAzimuthIndex].right;		break;
+		case 60:	psLeftTaps = pFilter44->e60[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e60[nAzimuthIndex].right;		break;
+		case 70:	psLeftTaps = pFilter44->e70[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e70[nAzimuthIndex].right;		break;
+		case 80:	psLeftTaps = pFilter44->e80[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e80[nAzimuthIndex].right;		break;
+		case 90:	psLeftTaps = pFilter44->e90[nAzimuthIndex].left;
+					psRightTaps = pFilter44->e90[nAzimuthIndex].right;		break;
+		};
+		
+		//How many taps to copy later to user's buffers
+		nTotalTaps = MIT_HRTF_44_TAPS;
+		break;
+	case 48000:
+		if(diffused)
+			pFilter48 = &diffuse_48;
+		else
+			pFilter48 = &normal_48;
+
+		switch(nInternalElevation)
+		{	
+		case -10:	psLeftTaps = pFilter48->e_10[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e_10[nAzimuthIndex].right;		break;
+		case -20:	psLeftTaps = pFilter48->e_20[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e_20[nAzimuthIndex].right;		break;
+		case -30:	psLeftTaps = pFilter48->e_30[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e_30[nAzimuthIndex].right;		break;
+		case -40:	psLeftTaps = pFilter48->e_40[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e_40[nAzimuthIndex].right;		break;
+		case 0:		psLeftTaps = pFilter48->e00[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e00[nAzimuthIndex].right;		break;
+		case 10:	psLeftTaps = pFilter48->e10[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e10[nAzimuthIndex].right;		break;
+		case 20:	psLeftTaps = pFilter48->e20[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e20[nAzimuthIndex].right;		break;
+		case 30:	psLeftTaps = pFilter48->e30[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e30[nAzimuthIndex].right;		break;
+		case 40:	psLeftTaps = pFilter48->e40[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e40[nAzimuthIndex].right;		break;
+		case 50:	psLeftTaps = pFilter48->e50[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e50[nAzimuthIndex].right;		break;
+		case 60:	psLeftTaps = pFilter48->e60[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e60[nAzimuthIndex].right;		break;
+		case 70:	psLeftTaps = pFilter48->e70[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e70[nAzimuthIndex].right;		break;
+		case 80:	psLeftTaps = pFilter48->e80[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e80[nAzimuthIndex].right;		break;
+		case 90:	psLeftTaps = pFilter48->e90[nAzimuthIndex].left;
+					psRightTaps = pFilter48->e90[nAzimuthIndex].right;		break;
+		};
+		
+		//How many taps to copy later to user's buffers
+		nTotalTaps = MIT_HRTF_48_TAPS;
+		break;
+	case 88200:
+		if(diffused)
+			pFilter88 = &diffuse_88;
+		else
+			pFilter88 = &normal_88;
+
+		switch(nInternalElevation)
+		{	
+		case -10:	psLeftTaps = pFilter88->e_10[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e_10[nAzimuthIndex].right;		break;
+		case -20:	psLeftTaps = pFilter88->e_20[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e_20[nAzimuthIndex].right;		break;
+		case -30:	psLeftTaps = pFilter88->e_30[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e_30[nAzimuthIndex].right;		break;
+		case -40:	psLeftTaps = pFilter88->e_40[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e_40[nAzimuthIndex].right;		break;
+		case 0:		psLeftTaps = pFilter88->e00[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e00[nAzimuthIndex].right;		break;
+		case 10:	psLeftTaps = pFilter88->e10[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e10[nAzimuthIndex].right;		break;
+		case 20:	psLeftTaps = pFilter88->e20[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e20[nAzimuthIndex].right;		break;
+		case 30:	psLeftTaps = pFilter88->e30[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e30[nAzimuthIndex].right;		break;
+		case 40:	psLeftTaps = pFilter88->e40[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e40[nAzimuthIndex].right;		break;
+		case 50:	psLeftTaps = pFilter88->e50[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e50[nAzimuthIndex].right;		break;
+		case 60:	psLeftTaps = pFilter88->e60[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e60[nAzimuthIndex].right;		break;
+		case 70:	psLeftTaps = pFilter88->e70[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e70[nAzimuthIndex].right;		break;
+		case 80:	psLeftTaps = pFilter88->e80[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e80[nAzimuthIndex].right;		break;
+		case 90:	psLeftTaps = pFilter88->e90[nAzimuthIndex].left;
+					psRightTaps = pFilter88->e90[nAzimuthIndex].right;		break;
+		};
+		
+		//How many taps to copy later to user's buffers
+		nTotalTaps = MIT_HRTF_88_TAPS;
+		break;
+	case 96000:
+		if(diffused)
+			pFilter96 = &diffuse_96;
+		else
+			pFilter96 = &normal_96;
+
+		switch(nInternalElevation)
+		{	
+		case -10:	psLeftTaps = pFilter96->e_10[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e_10[nAzimuthIndex].right;		break;
+		case -20:	psLeftTaps = pFilter96->e_20[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e_20[nAzimuthIndex].right;		break;
+		case -30:	psLeftTaps = pFilter96->e_30[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e_30[nAzimuthIndex].right;		break;
+		case -40:	psLeftTaps = pFilter96->e_40[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e_40[nAzimuthIndex].right;		break;
+		case 0:		psLeftTaps = pFilter96->e00[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e00[nAzimuthIndex].right;		break;
+		case 10:	psLeftTaps = pFilter96->e10[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e10[nAzimuthIndex].right;		break;
+		case 20:	psLeftTaps = pFilter96->e20[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e20[nAzimuthIndex].right;		break;
+		case 30:	psLeftTaps = pFilter96->e30[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e30[nAzimuthIndex].right;		break;
+		case 40:	psLeftTaps = pFilter96->e40[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e40[nAzimuthIndex].right;		break;
+		case 50:	psLeftTaps = pFilter96->e50[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e50[nAzimuthIndex].right;		break;
+		case 60:	psLeftTaps = pFilter96->e60[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e60[nAzimuthIndex].right;		break;
+		case 70:	psLeftTaps = pFilter96->e70[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e70[nAzimuthIndex].right;		break;
+		case 80:	psLeftTaps = pFilter96->e80[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e80[nAzimuthIndex].right;		break;
+		case 90:	psLeftTaps = pFilter96->e90[nAzimuthIndex].left;
+					psRightTaps = pFilter96->e90[nAzimuthIndex].right;		break;
+		};
+		
+		//How many taps to copy later to user's buffers
+		nTotalTaps = MIT_HRTF_96_TAPS;
+		break;
+	};
+
+	//Switch left and right ear if the azimuth is to the left of front centre (azimuth < 0)
+	if(nSwitchLeftRight)
+	{
+		psTempTaps = psRightTaps;
+		psRightTaps = psLeftTaps;
+		psLeftTaps = psTempTaps;
+	}
+	
+	//Copy taps to user's arrays
+	for(niTap = 0; niTap < nTotalTaps; niTap++)
+	{
+		psLeft[niTap] = psLeftTaps[niTap];
+		psRight[niTap] = psRightTaps[niTap];
+	}
+
+	//Assign the real azimuth and elevation used
+	*pAzimuth = nInternalAzimuth;
+	*pElevation = nInternalElevation;
+
+	return nTotalTaps;	
+}
+
+
+
+int mit_hrtf_findAzimuthFor40Elev(int azimuth)
+{
+	if(azimuth >= 0 && azimuth < 4)
+		return 0;
+	else if(azimuth >= 4 && azimuth < 10)
+		return 6;
+	else if(azimuth >= 10 && azimuth < 17)
+		return 13;
+	else if(azimuth >= 17 && azimuth < 23)
+		return 19;
+	else if(azimuth >= 23 && azimuth < 30)
+		return 26;
+	else if(azimuth >= 30 && azimuth < 36)
+		return 32;
+	else if(azimuth >= 36 && azimuth < 43)
+		return 39;
+	else if(azimuth >= 43 && azimuth < 49)
+		return 45;
+	else if(azimuth >= 49 && azimuth < 55)
+		return 51;
+	else if(azimuth >= 55 && azimuth < 62)
+		return 58;
+	else if(azimuth >= 62 && azimuth < 68)
+		return 64;
+	else if(azimuth >= 68 && azimuth < 75)
+		return 71;
+	else if(azimuth >= 75 && azimuth < 81)
+		return 77;
+	else if(azimuth >= 81 && azimuth < 88)
+		return 84;
+	else if(azimuth >= 88 && azimuth < 94)
+		return 90;
+	else if(azimuth >= 94 && azimuth < 100)
+		return 96;
+	else if(azimuth >= 100 && azimuth < 107)
+		return 103;
+	else if(azimuth >= 107 && azimuth < 113)
+		return 109;
+	else if(azimuth >= 113 && azimuth < 120)
+		return 116;
+	else if(azimuth >= 120 && azimuth < 126)
+		return 122;
+	else if(azimuth >= 126 && azimuth < 133)
+		return 129;
+	else if(azimuth >= 133 && azimuth < 139)
+		return 135;
+	else if(azimuth >= 139 && azimuth < 145)
+		return 141;
+	else if(azimuth >= 145 && azimuth < 152)
+		return 148;
+	else if(azimuth >= 152 && azimuth < 158)
+		return 154;
+	else if(azimuth >= 158 && azimuth < 165)
+		return 161;
+	else if(azimuth >= 165 && azimuth < 171)
+		return 167;
+	else if(azimuth >= 171 && azimuth < 178)
+		return 174;
+	else
+		return 180;
+};
+
+
+
+int mit_hrtf_findIndexFor40Elev(int azimuth)
+{
+	if(azimuth >= 0 && azimuth < 4)
+		return 0;
+	else if(azimuth >= 4 && azimuth < 10)
+		return 1;
+	else if(azimuth >= 10 && azimuth < 17)
+		return 2;
+	else if(azimuth >= 17 && azimuth < 23)
+		return 3;
+	else if(azimuth >= 23 && azimuth < 30)
+		return 4;
+	else if(azimuth >= 30 && azimuth < 36)
+		return 5;
+	else if(azimuth >= 36 && azimuth < 43)
+		return 6;
+	else if(azimuth >= 43 && azimuth < 49)
+		return 7;
+	else if(azimuth >= 49 && azimuth < 55)
+		return 8;
+	else if(azimuth >= 55 && azimuth < 62)
+		return 9;
+	else if(azimuth >= 62 && azimuth < 68)
+		return 10;
+	else if(azimuth >= 68 && azimuth < 75)
+		return 11;
+	else if(azimuth >= 75 && azimuth < 81)
+		return 12;
+	else if(azimuth >= 81 && azimuth < 88)
+		return 13;
+	else if(azimuth >= 88 && azimuth < 94)
+		return 14;
+	else if(azimuth >= 94 && azimuth < 100)
+		return 15;
+	else if(azimuth >= 100 && azimuth < 107)
+		return 16;
+	else if(azimuth >= 107 && azimuth < 113)
+		return 17;
+	else if(azimuth >= 113 && azimuth < 120)
+		return 18;
+	else if(azimuth >= 120 && azimuth < 126)
+		return 19;
+	else if(azimuth >= 126 && azimuth < 133)
+		return 20;
+	else if(azimuth >= 133 && azimuth < 139)
+		return 21;
+	else if(azimuth >= 139 && azimuth < 145)
+		return 22;
+	else if(azimuth >= 145 && azimuth < 152)
+		return 23;
+	else if(azimuth >= 152 && azimuth < 158)
+		return 24;
+	else if(azimuth >= 158 && azimuth < 165)
+		return 25;
+	else if(azimuth >= 165 && azimuth < 171)
+		return 26;
+	else if(azimuth >= 171 && azimuth < 178)
+		return 27;
+	else
+		return 28;
+}
\ No newline at end of file
diff -r 000000000000 -r 2076b4d80327 velocity.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/velocity.cpp	Tue Apr 07 21:09:22 2015 +0000
@@ -0,0 +1,54 @@
+#include "../include/velocity.h"
+#include <iostream>
+
+Velocity& Velocity::operator= (const Velocity& rhs)
+{
+    dx = (rhs.dx);
+    dy = (rhs.dy);
+    dz = (rhs.dz);
+    return *this;
+}
+
+Velocity& Velocity::operator+= (const Velocity& rhs)
+{
+    dx = dx + rhs.dx;
+    dy = dy + rhs.dy;
+    dz = dz + rhs.dz;
+    return *this;
+}
+
+Velocity operator+ (Velocity lhs, const Velocity& rhs)
+{
+    lhs += rhs;
+    return lhs;
+}
+
+bool Velocity::operator< (const Velocity& rhs) const
+{
+    if (dx == rhs.dx){
+        if (dy == rhs.dy){
+            return dz < rhs.dz;
+        }return dy < rhs.dy;
+    }return dx < rhs.dx;
+}
+
+std::ostream& operator<< (std::ostream& os, const Velocity& obj)
+{
+    os << obj.getdX() <<","<< obj.getdY() <<","<< obj.getdZ();
+    return os;
+}
+
+float Velocity::getdX (void) const
+{
+    return this->dx;
+}
+
+float Velocity::getdY (void) const
+{
+    return this->dy;
+}
+
+float Velocity::getdZ (void) const
+{
+    return this->dz;
+}