nivedita velagaleti
ese 519
Revision 0:2076b4d80327, committed 2015-04-07
- Comitter:
- niv17
- Date:
- Tue Apr 07 21:09:22 2015 +0000
- Commit message:
- sonic initial april 7
Changed in this revision
--- /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;
+}
--- /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;
+}
--- /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?
--- /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);
+}
+// -------------------------
--- /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]);
+ }
+ }
+ }
+}
+// -------------------------
+
--- /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;
+}
--- /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
--- /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;
+}