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

#include "io_ncamber.h"
#include "io_ncbase.h"

#ifdef ENABLE_TRAJECTORY_NC

int CNCAmberTraj::setupRead(CIO::SPrmtopHeadPointers &prmtopHeadPointer) {
    int result = -1;
    if (!this->trajfilename.empty()) {
        CNCAmberTraj::NC_openRead(this->trajfilename);
        if (this->ncid_ != -1) {
            if (CNetcdfBase::NCTYPE::NC_AMBERTRAJ == this->getNCConventions()) {
                this->setupFrameDim();
                this->setupCoords();
                if (this->Ncatom() != prmtopHeadPointer.natoms) {
                    cerr << "Number of atoms in topology and trajectory mismatch!" << endl;
                    exit(0);
                } else {
                    result = 0;
                }
            } else {
                cout << "Error unrecognized trajectory convention, it should be AMBER netcdf!" << endl;
                exit(0);
            }
        } else {
            cout << "Error opening netcdf file:: '" << this->trajfilename << "'" << endl;
            exit(0);
        }
    } else {
        cout << "The filename of the NetCDF is missing!" << endl;
        exit(0);
    }
    return result;
}

delphi_integer CNCAmberTraj::readAmberNC(CIO::SPrmtopHeadPointers prmtopHeadPointer, vector<CFrame> &vctOfFrames, vector<delphi_real> &vctAtomRadii,
                const int &iFrameFirst, const int &iFrameLast, const int &iFrameStride) {
    size_t start[3], count[3];
    delphi_integer natoms = 0;
    bool bReadFrames = false;
    if (this->ncid_ == -1) {
        if (setupRead(prmtopHeadPointer) == 0) {
            bReadFrames = true;
        }
    } else {
        bReadFrames = true;
    }
    if (bReadFrames) {
        int currentFrameIndex = iFrameFirst;
        while (currentFrameIndex < ncframe_) {
            bool bRecordFrame = this->shouldRecordFrame(currentFrameIndex, iFrameFirst, iFrameLast, iFrameStride);

            if (bRecordFrame) {
                start[0] = currentFrameIndex;
                count[0] = 1;
                start[1] = 0;
                count[1] = prmtopHeadPointer.natoms;
                start[2] = 0;
                count[2] = 3;

                CFrame newFrame(prmtopHeadPointer.natoms);
                vector<float> coord_tmp(prmtopHeadPointer.natoms * 3, 0.0);
                if (NCCheckErr(nc_get_vara_float(ncid_, coordVID_, start, count, coord_tmp.data()))) {
                    cerr << "Error: Reading coordinates data." << endl;
                    return 1;
                } else {
                    int vid = 0;
                    for (int atomIdx = 0; atomIdx < prmtopHeadPointer.natoms; atomIdx++) {
                        for (int axisIdx = 0; axisIdx < DIM; axisIdx++) {
                            switch (axisIdx) {
                            case 0:
                                newFrame.coordinates[atomIdx].nX = coord_tmp[vid];
                                break;
                            case 1:
                                newFrame.coordinates[atomIdx].nY = coord_tmp[vid];
                                ;
                                break;
                            case 2:
                                newFrame.coordinates[atomIdx].nZ = coord_tmp[vid];
                                ;
                                break;
                            }
                            vid++;
                        }
                    }
                    natoms = prmtopHeadPointer.natoms;
                    newFrame.computeFrameGeometry(vctAtomRadii);
                    //newFrame.dumpFrame();
                    vctOfFrames.push_back(newFrame);
                }

            }
            currentFrameIndex += iFrameStride;
        }
    }
    return natoms;
}
#endif // end block if nc is enabled
