#include "PDB.h"
#include "GDB.h"   

void printSyntax()
{
	cerr << "Calculate coordinates rmsd for PDB coordinates (12/2007)." << endl;
	cerr << "Options: "<< endl;

	cerr << " -pdb0      Reference PDB coordinate file " << endl;
	cerr << " -pdbs      Target PDB coordinate file(s)" << endl;
	cerr << " -atoms     set of atoms for calculating RMSD (default: CA,C',N)" << endl;
	cerr << " -seg       ranges of residues (a1 a2 b1 b2 ...) for calculating RMSD (default: All residues)" << endl;

	cerr << "Usage: " << endl;
	cerr << "  pdbrms -pdb0 a.pdb -pdbs b.pdb c.pdb -atoms CA -seg 6 16 30 66" << endl;
	cerr << "            Caluclate CA rmsd between three PDB coordinates for residues 6-16 and 30-66" << endl << endl;

	exit(0);
}

int main( int argc, char ** argv )
{
	string argList, pdbDir, slash_char;

	if(argc == 1) {
		printSyntax();
	}

	for(int i = 1; i < argc; i++){
		argList+= (argv[i]);
		argList+= (" ");
	}

	PDB pdb1;
	
	string atoms = "N CA C";
	
	vector<string> atomList, pdbList;
	atomList.push_back("N");
	atomList.push_back("CA");
	atomList.push_back("C");

	vector<int> segList;
	bool segFlag = false;

	string temp = getenv( "PATH" );
	if(temp.find("/") != string::npos ) slash_char = "/"; // unix
	else if(temp.find("\\") != string::npos ) slash_char = "\\"; // Windows 

	vector<string> fields = GDB::split("-", argList);
	for(int i = 0; i < fields.size(); i++){
		vector<string> temp;
		int pos = fields[i].find_first_of(' ');
		temp.push_back( fields[i].substr(0,pos) );
		temp.push_back( fields[i].substr(pos+1,fields[i].length()-pos-1 ));
		string arg = PDB::simplifyWhiteSpace(temp[1]);

		if(temp[0] == "pdb0" )
			pdb1.loadPDB(arg);
		else if(temp[0] == "pdbs" ) //pdb2.loadPDB(arg);
		{
			pdbList.clear();
                	vector<string> temp = GDB::split(" ", arg);
                	for(int i = 0; i < temp.size(); i++) 
                        	pdbList.push_back(temp[i]);
		}
		else if(temp[0] == "atoms" )
		{
			atomList.clear();
                	vector<string> temp = GDB::split(" ", arg);
                	for(int i = 0; i < temp.size(); i++) 
                        	atomList.push_back(temp[i]);
				
			atoms = arg;
		}
		else if(temp[0] == "seg" )
		{
			segList.clear();
			segFlag = true;
			
                	vector<string> temp = GDB::split(" ", arg);
			if( temp.size() % 2 == 1) 
			{
				cerr << "\t***Error: segment numbers must be paired -" << arg << endl;
                                exit(0);

			}
                	int cnt = 1;
                	for(int i = 0; i < temp.size(); i++) 
                        	segList.push_back( atoi(temp[i].c_str() ) );
		}
		else 
			printSyntax();
	}


	DoubleD2 vec_pdb1, vec_pdb2;
	int s_pdb1 = pdb1.getVectors(vec_pdb1, atomList, segList), s_pdb2;

	for(int i =0; i < pdbList.size(); i++)
	{
		vec_pdb2.clear();
		
		PDB pdb2(pdbList[i]);
		s_pdb2 = pdb2.getVectors(vec_pdb2, atomList, segList);
	
		if( s_pdb2 == s_pdb1 )
			printf("%6.3f\t\t%6s\t\t%6s\n", pdb2.rmsd(vec_pdb1,vec_pdb2,s_pdb1), pdb2.PDBfileName.c_str(), pdb1.PDBfileName.c_str() );
		else 
			cout << " ***Error: atom sets " << atoms << " not indentical for "<< pdb1.PDBfileName << " and " << pdb2.PDBfileName << endl;
	}		
}