#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <math.h>
using namespace std;

// this program will parse two extinction profiles generated by UltraScan-II's extinction fitter
// and return the angle between the two spectra. The calculation made is cos x = (u*v)/(||u||*||v||)
// where u and v are the two spectrum vectors.

int maxsize;
struct spectrum
{
   const char * filename;
	vector <double> lambda, extinction;
	int arraysize;
};

void parse(struct spectrum *);
void check_overlap(vector <struct spectrum> *);
double calc_angle(vector <struct spectrum> *, int, int);

int main () 
{
   vector <struct spectrum> spec;
	spectrum tmp_spec;
	spec.clear();
   unsigned int i,j;
	tmp_spec.lambda.clear();
	tmp_spec.extinction.clear();
   tmp_spec.filename = "CytC-oxidized.extinction.dat";
	spec.push_back(tmp_spec);
   tmp_spec.filename = "CytC-reduced.extinction.dat";
	spec.push_back(tmp_spec);
   for (j=0; j<spec.size(); j++)
   {
      parse (&spec[j]);
   }
	check_overlap(&spec);
	cout << calc_angle(&spec, 0, 1) << " degrees)" << endl;
   return 0;
}

void parse(struct spectrum *tmp_spec)
{
   ifstream myfile;
   string line;
   char seps[]   = " ,\t\n";
   char *token;
   int i=0, j=0;
   myfile.open (tmp_spec->filename);
   if (myfile.is_open())
   {
		getline (myfile, line); // get rid of first line, which is text
      while ( getline (myfile, line) )
      {
         //cout << line << '\n';
         char *string = new char[line.length() + 1];
         strcpy(string, line.c_str());
         token = strtok( string, seps );
			tmp_spec->lambda.push_back(atof(token));
			token = strtok( NULL, seps );
			tmp_spec->extinction.push_back(atof(token));
         delete [] string;
      }
   }
   myfile.close();
}

void check_overlap(vector <struct spectrum> *spec)
{
	maxsize=2147483647;
	int j;
	double lowest=-1.0;
	for (j=0; j<spec->size(); j++)
	{
		lowest = max(lowest, (*spec)[j].lambda[0]); // find the lowest value common to all spectra
	}
	for (j=0; j<spec->size(); j++)
	{
		while (lowest > (*spec)[j].lambda[0])
		{
			(*spec)[j].lambda.erase((*spec)[j].lambda.begin());
			(*spec)[j].extinction.erase((*spec)[j].extinction.begin());
		}
		(*spec)[j].arraysize = (*spec)[j].lambda.size();
	}
	for (j=0; j<spec->size(); j++)
	{
	//	cout << (*spec)[j].filename << ": " << (*spec)[j].lambda[0] << ", ";
	//	cout << (*spec)[j].arraysize << endl;
		maxsize = min(maxsize, (*spec)[j].arraysize);
	}
}

double calc_angle(vector <struct spectrum> *spec, int vec1, int vec2)
{
	int i, j;
	double val1=0.0, val2=0.0, val3=0.0, tmp;
	for (i=0; i<maxsize; i++)
	{
		val1 += (*spec)[vec1].extinction[i] * (*spec)[vec2].extinction[i];
		val2 += (*spec)[vec1].extinction[i] * (*spec)[vec1].extinction[i];
		val3 += (*spec)[vec2].extinction[i] * (*spec)[vec2].extinction[i];
	}
	tmp = val1/(pow(val2, 0.5) * pow(val3, 0.5));
	cout << "cos x between " << (*spec)[vec1].filename << " and " << (*spec)[vec2].filename << ": "<< tmp << " (";
	return(180 * (acos(tmp))/M_PI);
}