//
/*
 * io_ambertraj.cpp
 *
 *  Created on: Apr 18, 2019
 *      Author: Shailesh
 */

#include "io.h"
#include <vector>

#define DIM 3    // number of coordinate dimensions

using namespace std;

delphi_integer CIO::readMDCRD(ifstream& ifMdcrdFileStream,
		SPrmtopHeadPointers prmtopHeadPointer, vector<CFrame>& vctOfFrames,
		vector<delphi_real>& vctAtomRadii, const int& iFrameFirst,
		const int& iFrameLast, const int& iFrameStride) 
{
	
	int valuesPerFrame = (prmtopHeadPointer.natoms * DIM);
	#ifdef VERBOSE
	cout << "io_mdcrd> prmtopHeadPointer.natoms = " << prmtopHeadPointer.natoms << endl;
	#endif

	int linesPerFrame = ceil(valuesPerFrame / 10.0);

	int firstFrameOffset = 1; // trajectory title spans first line of mdcrd file
	int currentLineIndex = 0;

	string strLine;
	int currentFrameIndex = 0;
	
	// read title line of the trajectory file
	getline(ifMdcrdFileStream, strLine);
	
	while (!ifMdcrdFileStream.eof()) 
	{
		bool bRecordFrame = shouldRecordFrame(currentFrameIndex, iFrameFirst, iFrameLast, iFrameStride);
		
		if (bRecordFrame) 
		{
			CFrame newFrame(prmtopHeadPointer.natoms);
			int valueIndex = 0;
			int valuesToRead = valuesPerFrame;
			
			for (int i = 0; i < linesPerFrame; ++i) 
			{
				getline(ifMdcrdFileStream, strLine);
				int valuesInThisLine = (valuesToRead > 10) ? 10 : valuesToRead;
				
				// This section is exclusively for trajectories containg box but box missing in prmtop
				if (prmtopHeadPointer.ifbox == 0 &&  i == 0 && ( strlen(strLine.c_str()) != valuesInThisLine * 8 ) )
					throw CBadMDCRDFile("INCOMPATIBILITY DETECTED: PRMTOP has no box but MDCRD possibly has box information.");
			

				for (int j = 0; j < valuesInThisLine; ++j, ++valueIndex) 
				{
					int atomIdx = valueIndex / DIM;
					int axisIdx = valueIndex % DIM;
					double value = atof(
							removeSpace(strLine.substr(j * 8, 8)).c_str());
					switch (axisIdx) 
					{
						case 0:
							newFrame.coordinates[atomIdx].nX = value;
							break;
						case 1:
							newFrame.coordinates[atomIdx].nY = value;
							break;
						case 2:
							newFrame.coordinates[atomIdx].nZ = value;
							break;
					}
					valuesToRead--;
				}// for j 0 to valuesInThisLine

			} // for i 0 to linesPerFrame

			newFrame.computeFrameGeometry(vctAtomRadii);
			//newFrame.dumpFrame();
			vctOfFrames.push_back(newFrame);
		} 
		else 
		{
			// skip this frame
			for (int i = 0; i < linesPerFrame; ++i) 
			{
				getline(ifMdcrdFileStream, strLine);
			}
		}
		
		// if box information is present skip it
		if (prmtopHeadPointer.ifbox) 
			getline(ifMdcrdFileStream, strLine);
		
		currentFrameIndex++;
		if (currentFrameIndex == iFrameLast)
			break;
	}

	return prmtopHeadPointer.natoms;
}
