//
//  io_dcd.cpp
//
//  Created by Arghya Chakravorty on 1/23/19.
//
#include "io.h"

#define SUCCESS 0
#define BAD_DCD -1

void CIO::readBlockSize(ifstream& fileStream, int& blockSize)
{
    //fileStream.read( (char *)&blockSize, sizeof(int) );
    blockSize = this->decode<int>(fileStream);    
}

void CIO::skipBytes(ifstream& fileStream, int& nbytes)
{
    fileStream.seekg(nbytes, ios_base::cur);
}

int CIO::readDCD(ifstream& ifDcdFileStream, vector<CFrame>& vctOfFrames, vector<delphi_real>& vctAtomRadii, const int& iFrameFirst, const int& iFrameLast, const int& iFrameStride)
{
    /*
     * In the DCD file, every block has an integer that starts it and the same integer that marks
     * its end. this integer is the bytesize or size of that block, i.e. the number of bytes of information
     * to be expected till before the end of the block is reached.
     * */
    
    // local variables
    int natoms;
    long int totalSize;
    int idx, jdx, kdx, iFrame = 0;
    int startBlockSize, endBlockSize;
    float x, y, z;
    
    // gather size of the dcd file
    ifDcdFileStream.seekg(0,ios_base::end);
    totalSize = ifDcdFileStream.tellg();
    
    // reset to the beggining
    ifDcdFileStream.seekg(0,ios_base::beg);
    
    /**********************************
     *      HEADER BLOCK
     * read and skip dcd header
     ***********************************/
    readBlockSize(ifDcdFileStream, startBlockSize);
    // cout << " Starting size is : " << startBlockSize << " bytes " << endl;
    
    if ( startBlockSize != 84 ) return BAD_DCD;
    
    // skip the next `startBlockSize` bytes
    skipBytes(ifDcdFileStream, startBlockSize);
    
    readBlockSize(ifDcdFileStream, endBlockSize);
    // cout << " Ending size is : " << endBlockSize << " bytes " << endl;
    
    if ( startBlockSize != endBlockSize ) cerr << " Bad reading operation of the header block " << endl;
    
    
    /**********************************
     *      TITLE BLOCK
     * read and skip TITLE block
     ***********************************/
    readBlockSize(ifDcdFileStream, startBlockSize);
    // cout << " Starting size is : " << startBlockSize << " bytes " << endl;
    
    // skip the next `startBlockSize` bytes
    skipBytes(ifDcdFileStream, startBlockSize);
    
    readBlockSize(ifDcdFileStream, endBlockSize);
    // cout << " Ending size is : " << endBlockSize << " bytes " << endl;
    
    if ( startBlockSize != endBlockSize ) cerr << " Bad reading operation of the title block " << endl;
    
    
    /**********************************
     *   NUMBER OF ATOMS BLOCK
     ***********************************/
    readBlockSize(ifDcdFileStream, startBlockSize);
    // cout << " Starting size is : " << startBlockSize << " bytes " << endl;
    
    // fetch number of atoms
    natoms = this->decode<int>(ifDcdFileStream);
    // ifDcdFileStream.read( (char *) &natoms, sizeof(int) );
    cout << " Each DCD frame is expected to have " << natoms << " atoms" << endl;
    
    readBlockSize(ifDcdFileStream, endBlockSize);
    // cout << " Ending size is : " << endBlockSize << " bytes " << endl;
    
    if ( startBlockSize != endBlockSize ) cerr << " Bad reading operation of the #atoms block " << endl;
    
    /**********************************
     *   COORDINATES PER FRAME BLOCK
     ***********************************/
    while ( ifDcdFileStream.tellg() != totalSize)
    {
        // could be the start of the unit cell information or X-block
        readBlockSize(ifDcdFileStream, startBlockSize);
        // cout << " Starting size is : " << startBlockSize << " bytes " << endl;
        bool bRecordFrame = shouldRecordFrame(iFrame, iFrameFirst, iFrameLast, iFrameStride);
        // create an instance of Class CFrame
        CFrame newFrame(natoms);
        
        
        // has unit cell information
        if ( startBlockSize == 48 )
        {
            // there is unit cell dimension information
            // skip 48 bytes
            skipBytes(ifDcdFileStream, startBlockSize);
            
            readBlockSize(ifDcdFileStream, endBlockSize);
            // cout << " Ending size is : " << endBlockSize << " bytes " << endl;
            
            // start of the X-block if unit cell information is present
            readBlockSize(ifDcdFileStream, startBlockSize);
            // cout << " Starting size for x block : " << startBlockSize << " bytes " << endl;
            
        }
        
        if (startBlockSize == sizeof(float)*natoms)
        {
            // read X which are consecutively placed for each atoms
            for (idx = 0; idx < natoms; idx++)
            {
                // ifDcdFileStream.read( (char *) &x, sizeof(float) );
                // newFrame.coordinates[idx].nX = x;
                newFrame.coordinates[idx].nX = this->decode<float>(ifDcdFileStream);
                //cout << newFrame.coordinates[idx].nX << endl;
            }
            // end of X-Block
            readBlockSize(ifDcdFileStream, endBlockSize);
            //--------------------------------------------------------------------------//


            // read Y which are consecutively placed for each atoms.
            // Here, the startBlockSize needs to be read again
            readBlockSize(ifDcdFileStream, startBlockSize);
            
            for (jdx = 0; jdx < natoms; jdx++)
            {
                // ifDcdFileStream.read( (char *) &y, sizeof(float) );
                // newFrame.coordinates[jdx].nY = y;
                newFrame.coordinates[jdx].nY = this->decode<float>(ifDcdFileStream);
            }
            // end of Y-Block
            readBlockSize(ifDcdFileStream, endBlockSize);
            //--------------------------------------------------------------------------//

            
            // read Z which are consecutively placed for each atoms.
            // Here, the startBlockSize needs to be read again
            readBlockSize(ifDcdFileStream, startBlockSize);
            
            for (kdx = 0; kdx < natoms; kdx++)
            {
                // ifDcdFileStream.read( (char *) &z, sizeof(float) );
                // newFrame.coordinates[kdx].nZ = z;
                newFrame.coordinates[kdx].nZ = this->decode<float>(ifDcdFileStream);
            }
            // end of Z-Block
            readBlockSize(ifDcdFileStream, endBlockSize);
            //--------------------------------------------------------------------------//
           
            
        }
        if (bRecordFrame) {
            newFrame.computeFrameGeometry(vctAtomRadii);
            //newFrame.dumpFrame();
            vctOfFrames.push_back(newFrame);
        }
        iFrame++;
    }
    return natoms;
}
