//
/*
 * io_trr.cpp
 *
 *  Created on: Feb 7, 2019
 *      Author: Arghya 
 */

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

#define DIM 3
#define GROMACS_MAGIC 1993
#define NM2ANG 10.

using namespace std;

int CIO::readTRR(ifstream& ifTrrFileStream, vector<CFrame>& vctOfFrames, vector<delphi_real>& vctAtomRadii, const int& iFrameFirst, const int& iFrameLast, const int& iFrameStride) 
{
   
    int natoms = -1;
    int iFrame = 0;
    
    // gather size of the trr file
    ifTrrFileStream.seekg(0,ios_base::end);
    long int totalSize = ifTrrFileStream.tellg();
    
    // reset to the beggining
    ifTrrFileStream.seekg(0,ios_base::beg);
    
    // read the file and generate new frames 
    while (ifTrrFileStream.tellg() != totalSize) {
        
        // endianness check using MAGIC NUMBER
        int varInt = this->decode<int>(ifTrrFileStream);
        while ( varInt != GROMACS_MAGIC ) {
            this->bSwapBytes = true;
            varInt = this->swapbytes<int>(varInt);

            if ( varInt != GROMACS_MAGIC ) throw CBadTRR(); 
        }
        
        // version
        int version = this->decode<int>(ifTrrFileStream);
        
        // title
        string title(this->decode_gro_string(ifTrrFileStream));
        
        // sizes and other information
        int ir_size = this->decode<int>(ifTrrFileStream);
        int e_size = this->decode<int>(ifTrrFileStream);
        int box_size = this->decode<int>(ifTrrFileStream);
        int vir_size = this->decode<int>(ifTrrFileStream);
        int pres_size = this->decode<int>(ifTrrFileStream);
        int top_size = this->decode<int>(ifTrrFileStream);
        int sym_size = this->decode<int>(ifTrrFileStream);
        
        int sizes[7] = {ir_size, e_size, box_size, vir_size, pres_size, top_size, sym_size };
        
        // read `x_size, v_size, f_size`
        int x_size = this->decode<int>(ifTrrFileStream);
        int v_size = this->decode<int>(ifTrrFileStream);
        int f_size = this->decode<int>(ifTrrFileStream);
        
        // natoms
        natoms = this->decode<int>(ifTrrFileStream);
        
        // skip step and nre
        skip_n<int>(ifTrrFileStream, 2);
        
        // time tag of this snapshot
        float t = this->decode<float>(ifTrrFileStream);
        // cout << " Time : " << t << endl;
        
        // skip lambda (float)
        skip_n<float>(ifTrrFileStream, 1);
        
        
        // skip the number of bytes contained in sizes
        for ( int k = 0; k < 7; k++) skipBytes(ifTrrFileStream, sizes[k]); 
        
        // create an instance of Cframe if the index of this frame fits the criteria set by first, last and stride parameters
        bool bRecordFrame = shouldRecordFrame(iFrame, iFrameFirst, iFrameLast, iFrameStride);
        
        CFrame newFrame(natoms);
        
		#ifdef VERBOSE
        if (bRecordFrame)
            cout << "io_trr> Reading frame : " << iFrame << endl;
        #endif

        // coordinates
        int precision = (int) (x_size/(3 * natoms));
        int idx = 0;
        
        if ( precision == sizeof(float))
        {
            if (bRecordFrame)
            {
                while ( idx != natoms )
                {
                    newFrame.coordinates[idx].nX = NM2ANG * this->decode<float>(ifTrrFileStream);
                    newFrame.coordinates[idx].nY = NM2ANG * this->decode<float>(ifTrrFileStream);
                    newFrame.coordinates[idx].nZ = NM2ANG * this->decode<float>(ifTrrFileStream);
                    
                    idx++;
                }
            }
            else
                skip_n<float>(ifTrrFileStream, DIM * natoms);

            if ( v_size > 0 ) skip_n<float>(ifTrrFileStream, DIM * natoms);
            if ( f_size > 0 ) skip_n<float>(ifTrrFileStream, DIM * natoms);
            
        }// if precision == float
        else if  ( precision == sizeof(double))
        {
            if (bRecordFrame)
            {
                while ( idx != natoms )
                {
                    newFrame.coordinates[idx].nX = NM2ANG * this->decode<double>(ifTrrFileStream);
                    newFrame.coordinates[idx].nY = NM2ANG * this->decode<double>(ifTrrFileStream);
                    newFrame.coordinates[idx].nZ = NM2ANG * this->decode<double>(ifTrrFileStream);
                    
                    idx++;
                    
                } // while idx != natoms
            }
            else
                skip_n<double>(ifTrrFileStream, DIM * natoms);
            
            if ( v_size > 0 ) skip_n<double>(ifTrrFileStream, DIM * natoms);
            if ( f_size > 0 ) skip_n<double>(ifTrrFileStream, DIM * natoms);

        } // if precision == double
        
        if (bRecordFrame)
        {
            newFrame.computeFrameGeometry(vctAtomRadii);
            //newFrame.dumpFrame();
            vctOfFrames.push_back(newFrame);
        }
        
        // break if current frame is the last frame
        if (iFrame == iFrameLast) break;
        
        // increment the index of the frame read/skipped
        iFrame++;

        // cout << "^^^^^^^^^^^^^^^^^^ frame : " << iFrame << " ^^^^^^^^^^^^^^^^^^^^^^^" << endl;
        // cout << " Now at " << ifTrrFileStream.tellg() << "/" << totalSize << endl;
        
    } // while (ifTrrFileStream.tellg() != totalSize)
    
    ifTrrFileStream.close();
    
    
    return natoms;
} // END OF FUNCTION


