DTGrid3D.hpp

Go to the documentation of this file.

 
#ifndef DTGRID3D_HPP_
#define DTGRID3D_HPP_
//
#include <stdlib.h>
#include <string.h>
#include <cstring>
#include <iostream>
#include <fstream>
#include <sstream>
#include <typeinfo>
//
#include "ArrayReg.hpp"
#include "ArrayPtr.hpp"
#include "ArrayPrs.hpp"
#include "DTGrid2D.hpp"
 
template<typename T,typename K>
class DTGrid3D {
public:
    DTGrid3D(char *, T, T);
    ~DTGrid3D();
    
private:
 
    ArrayReg <T> * Value_List_3D;
    ArrayPrs <K> * Zcoord_List;
    ArrayPtr <T,K> * Ptr_List_3D;
 
    DTGrid2D <T,K> * proj2D;
 
    int Length_of_Value_List_3D;
    int Length_of_Zcoord_List;
    int Length_of_Ptr_List_3D;
    
    T gamma_i;
    T gamma_o;
 
    int iform;
 
    int grid_type;
    
    T scfct; 
    float scale;
    float spacing;
    
    
    /*
    typename ArrayPrs <K>::PairArray * coord3;
    T ** ptt3;
    typename ArrayPrs <K>::PairArray * coo3;
    bool key3;
    short skey1;
    short skey2;
    short skey3;
    K checkpoint_31;
    K checkpoint_32;
    K countdown_3;
    */
 
    struct {
        int i;
        int j;
        int k;
    }dimensions;
    
    struct {
        float x;
        float y;
        float z;
    }origin;
    
    //
    void set_value_list_3D(T*,int);
    void set_zcoord_list(K*,int);
    void set_ptr_list_3D(K*,int);
    void set_gamma_values(T, T);
    void set_projection2D(K*, K*, K*, K*, K*, K*, int, int, int, int, int, int);
    void set_uhbd_parameters(float, float, float, float, float, int, int, int);
    void Reset_pointers_3D(void);
 
    //
    void test_format_dtgrid(char *);
    void size_ds(unsigned int);
    void format_ascii(char *);
    void format_binary(char *);
 
public:
    struct {
        char * Title;
        int grdflag;
        int one;
        int idummy1;
        int idummy2;
        int idummy3;
        int idummy4;
        int idummy5;
        float scale;
        float spacing;
        float dummy1;
        float dummy2;
        float dummy3;
        float dummy4;
        float dummy5;
        float dummy6;
        float dummy7;
        struct dimensions
        {
            int i;
            int j;
            int k;
        }dim;
        struct origin
        {
            float x;
            float y;
            float z;
        }o;
    }uhbd;
    T accessxyz(int, int, int);
    bool excl_accessxyz(int, int, int);
    T access_fast(int, int, int);
    T * getCell(T *, int, int, int);
    T * getCell2(T *, int, int, int);
 
    int type_of_grid (void) {return grid_type;}
protected:
 
    //template <class M, class H>
    //friend class DTGrid2D;
    //template <class M, class H>
    //friend class DTGrid1D;
 
    void getCellz(T *, int, int, int, typename ArrayPrs <K>::PairArray *);
    void getCellyz(T *, int, int, int, typename ArrayPrs <K>::PairArray *);
 
    T access_fast_z(int, typename ArrayPrs <K>::PairArray *);
    T access_fast_yz(int, int, typename ArrayPrs <K>::PairArray *);
 
};
template <typename T, typename K>
DTGrid3D<T,K>::~DTGrid3D(){
 
    if (grid_type == 0){
        delete proj2D;
        delete Zcoord_List;
    }
    else if (grid_type == 1){
        delete proj2D;
        delete Value_List_3D;
        delete Zcoord_List;
        delete Ptr_List_3D;
    }
    else{
        // Undefined
    }
 
}
 
 
template <typename T, typename K>
DTGrid3D<T,K>::DTGrid3D(char * file_name, T scfct_i, T gamma){
    scfct = scfct_i;
    set_gamma_values(gamma, gamma);
    test_format_dtgrid(file_name);
    if (iform == 0)
        format_binary(file_name);
    else if (iform == 1)
        format_ascii(file_name);
    else{
        std::cout << "Error: Unknown file format" << std::endl;
        exit(0);
    }
}
 
 
template <typename T, typename K>
void DTGrid3D<T,K>::set_value_list_3D(T* values, int length){
    if (grid_type == 0){
        // No values. All True if grid_type = 0
        Length_of_Value_List_3D = 0;
    }
    else if (grid_type == 1){
        Value_List_3D = new (std::nothrow) ArrayReg <T>;
        Value_List_3D->Set_Regular_Array_Values(values,length);
        Length_of_Value_List_3D = length;
    }
    else{
        // Undefined
    }
}
 
 
template <typename T, typename K>
void DTGrid3D<T,K>::set_zcoord_list(K* values, int length){
    Zcoord_List = new (std::nothrow) ArrayPrs <K>;
    Zcoord_List->Set_Pair_Array_Values(values,length);
    Length_of_Zcoord_List = length;
}
 
 
template <typename T, typename K>
void DTGrid3D<T,K>::set_ptr_list_3D(K* values, int length){
    if (grid_type == 0){
        // No values. All True if grid_type = 0
        Length_of_Ptr_List_3D = 0;
    }
    else if (grid_type == 1){
        Ptr_List_3D = new (std::nothrow) ArrayPtr <T,K>;
        Ptr_List_3D->Set_Pointer_Array_Values(Value_List_3D->Array,values,Length_of_Value_List_3D,length);
        Length_of_Ptr_List_3D = length;
    }
    else{
        // Undefined
    }
}
 
 
template <typename T, typename K>
void DTGrid3D<T,K>::set_projection2D(K* values1, K* values2, K* values3, \
        K* values4, K* values5, K* values6, int length1, int length2, \
        int length3, int length4, int length5, int length6){
 
    proj2D = new (std::nothrow) DTGrid2D <T,K>;
    proj2D->set_value_list_2D(values1,length1);
    proj2D->set_ycoord_list(values2,length2);
    proj2D->set_ptr_list_2D(values3,length3);
    proj2D->set_projection1D(values4, values5, values6, length4, length5, length6);
}
 
 
/*template <typename T, typename K>
void DTGrid3D<T,K>::Reset_pointers_3D(void){
    Value_List_3D->Reset_Regular_Array_Values();
    Zcoord_List->Reset_Set_Pair_Array_Values();
    Ptr_List_3D->Reset_Set_Pointer_Array_Values();
    proj2D->Reset_pointers_2D();
}*/
 
 
template <typename T, typename K>
void DTGrid3D<T,K>::set_uhbd_parameters(float sc, float sp, float ox, float oy, float oz, int dimi, int dimj, int dimk){
    scale = sc;
    spacing = sp;
    
    origin.x = ox;
    origin.y = oy;
    origin.z = oz;
    
    dimensions.i = dimi;
    dimensions.j = dimk;
    dimensions.k = dimk;
}
 
 
template <typename T, typename K>
void DTGrid3D<T,K>::set_gamma_values(T g_i, T g_o){
    gamma_i = g_i;
    gamma_o = g_o;
}
 
 
template <typename T, typename K>
T DTGrid3D<T,K>::access_fast(int a, int b, int c){
    
    bool key1 = false; // a local logical operation variable.
 
    int countdown_1 = proj2D->proj1D->set_countdown_1; // a local counter to iterate over the connected components in a tubular grid.
    typename ArrayPrs <K>::PairArray ** ptt1 = proj2D->proj1D->Ptr_List_1D->PtrArray; // a temporary pointer to store an acc array member.
    typename ArrayPrs <K>::PairArray * coo1 = proj2D->proj1D->Xcoord_List->PArray; // a temporary pointer to store an xCoord array pair.
 
    // iterate over the connected components in DTGrid1D
    for (;countdown_1>0;countdown_1--){
        // check whether index a is in one of the connected components using the xCoord index array.
        key1 = a>=coo1->v1 ? (a<=coo1->v2 ? true : false) : false;
        if (key1 == true){
            return access_fast_yz(b, c, *(ptt1) + a - coo1->v1);
        }
        else {
            ++ptt1;
            ++coo1;
            continue;
        }
    }
    return gamma_o;
}
 
 
template <typename T, typename K>
T DTGrid3D<T,K>::access_fast_yz(int b, int c, typename ArrayPrs <K>::PairArray * coord2){
 
    bool key2 = false; // a local logical operation variable.
 
    K checkpoint_21 = coord2->v2; // the index of the first connected component in a tubular grid in the yCoord array.
    K checkpoint_22 = (coord2+1)->v2; // the index of the last connected component in a tubular grid in the yCoord array.
    K countdown_2 = checkpoint_22 - checkpoint_21; // Determine how many connected components are there in the corresponding tubular grid.
 
    // Return when no connected component found
    if (countdown_2 == 0){ return gamma_o;}
 
    typename ArrayPrs <K>::PairArray ** ptt2 = proj2D->Ptr_List_2D->PtrArray+checkpoint_21; // a temporary pointer to store an acc array member.
    typename ArrayPrs <K>::PairArray * coo2 = proj2D->Ycoord_List->PArray+checkpoint_21; // a temporary pointer to store an yCoord array pair.
 
    // iterate over the connected components in DTGrid2D
    for (; countdown_2 >= 1; countdown_2--){
        // check whether index b is inside one of the connected components using the yCoord index array.
        key2 = b>=coo2->v1 ? (b<=coo2->v2 ? true : false) : false;
        if (key2 == true){
            return access_fast_z(c, *(ptt2) + b - coo2->v1);
        }
        else {
            ++ptt2;
            ++coo2;
            continue;
        }
    }
    // Return if nothing found
    return gamma_o;
}
 
 
template <typename T, typename K>
T DTGrid3D<T,K>::access_fast_z(int c, typename ArrayPrs <K>::PairArray * coord3){
 
    bool key3 = false; // a local logical operation variable.
 
    K checkpoint_31 = coord3->v2; // the index of the first connected component in a tubular grid in the zCoord array.
    K checkpoint_32 = (coord3+1)->v2; // the index of the last connected component in a tubular grid in the zCoord array.
    K countdown_3 = checkpoint_32 - checkpoint_31; // Determine how many connected components are there in the corresponding tubular grid (3D).
 
    if (countdown_3 == 0) return gamma_o;
 
    T ** ptt3 = Ptr_List_3D->PtrArrayReg+checkpoint_31; // a temporary pointer to store an acc array member.
    typename ArrayPrs <K>::PairArray * coo3 = Zcoord_List->PArray+checkpoint_31; // a temporary pointer to store an zCoord array pair.
 
    //int count = check;
    // iterate over the connected components in DTGrid3D
    for (; countdown_3 >= 1; countdown_3--){
        // check whether index c is inside one of the connected components using the zCoord index array.
        key3 = c>=coo3->v1 ? (c<=coo3->v2 ? true : false) : false;
        if (key3 == true){
            return *(*(ptt3)+c-coo3->v1);
        }
        else {
            //if (count != 1) {return_type = c<coo3->v1 ? (c>(coo3-1)->v2 ? gamma_i : gamma_o) : gamma_o;}
            ++ptt3;
            ++coo3;
            continue;
        }
    }
    return gamma_o;
}
 
 
template <typename T, typename K>
T DTGrid3D<T,K>::accessxyz(int a, int b, int c){
    
    // Call the random access function in the DTGrid2D.
    typename ArrayPrs <K>::PairArray * coord3 = proj2D->accessxy(a,b); // the index of the first connected component in a tubular grid in the zCoord array.
 
    // Return the default value if recursive function call returns NULL.
    if (coord3 == NULL) return gamma_o; 
 
    bool key3 = false; // a local logical operation variable.
    //K return_type = gamma_o;
 
    K checkpoint_31 = coord3->v2; // the index of the first connected component in a tubular grid in the zCoord array.
    K checkpoint_32 = (coord3+1)->v2; // the index of the last connected component in a tubular grid in the zCoord array.
    K countdown_3 = checkpoint_32 - checkpoint_31; // Determine how many connected components are there in the corresponding tubular grid (3D).
    
    // Return when no connected component found
    if (countdown_3 == 0) return gamma_o; 
    
    T ** ptt3 = Ptr_List_3D->PtrArrayReg+checkpoint_31; // a temporary pointer to store an acc array member.
    typename ArrayPrs <K>::PairArray * coo3 = Zcoord_List->PArray+checkpoint_31; // a temporary pointer to store an zCoord array pair.
    
    //int count = check;
    // iterate over the connected components in DTGrid3D
    for (; countdown_3 >= 1; countdown_3--){
        // check whether index c is inside one of the connected components using the zCoord index array.
        key3 = c>=coo3->v1 ? (c<=coo3->v2 ? true : false) : false;
        if (key3 == true){
            return *(*(ptt3)+c-coo3->v1);
        }
        else {
            //if (count != 1) {return_type = c<coo3->v1 ? (c>(coo3-1)->v2 ? gamma_i : gamma_o) : gamma_o;}
            ++ptt3;
            ++coo3;
            continue;
        }
    }
    return gamma_o;
}
 
 
/*template <typename T, typename K>
bool DTGrid3D<T,K>::excl_accessxyz(int a, int b, int c){
    typename ArrayPrs <K>::PairArray * coord3 = proj2D->accessxy(a,b);
    if (coord3 == NULL) return false;
    bool key3 = false;
    K checkpoint_31 = coord3->v2;
    K checkpoint_32 = (coord3+1)->v2;
    K countdown_3 = checkpoint_32 - checkpoint_31;
    if (countdown_3 == 0) return false;
    typename ArrayPrs <K>::PairArray * coo3 = Zcoord_List->PArray+checkpoint_31;
    for (; countdown_3 >= 1; countdown_3--){
        key3 = c>=coo3->v1 ? (c<=coo3->v2 ? true : false) : false;
        if (key3 == true){
            return true;
        }
        else {
            ++coo3;
            continue;
        }
    }
    return false;
}*/
 
 
template <typename T, typename K>
void DTGrid3D<T,K>::getCellz(T * Cell, int x, int y, int c, typename ArrayPrs <K>::PairArray * coord3){
 
    short skey3; // a local integer variable to identify the grid cell type. 
    int cinc = c+1; // the index of the origin's neighbor in z direction.
 
    K checkpoint_31 = coord3->v2; // the index of the first connected component in a tubular grid in the zCoord array.
    K checkpoint_32 = (coord3+1)->v2; // the index of the last connected component in a tubular grid in the zCoord array.
    K countdown_3 = checkpoint_32 - checkpoint_31; // Determine how many connected components are there in the corresponding tubular grid (3D).
 
    T ** ptt3 = Ptr_List_3D->PtrArrayReg+checkpoint_31; // a temporary pointer to store an acc array member.
    typename ArrayPrs <K>::PairArray * coo3 = Zcoord_List->PArray+checkpoint_31; // a temporary pointer to store an zCoord array pair.
 
    // iterate over the connected components in DTGrid3D
    for (; countdown_3 >= 1; countdown_3--){
        // check whether index c (and/or c+1) is inside one of the connected components using the zCoord index array.
        skey3 = c>=coo3->v1 ? (cinc<=coo3->v2 ? 1 : (c<=coo3->v2 ? 2 : 0) ) : (cinc>=coo3->v1 ? (cinc<=coo3->v2 ? 3 : 0) : 0);
        switch (skey3)
        {
        // Neither c, nor c+1 is inside the connected component.
        case 0:
            ++ptt3;
            ++coo3;
            break;
        // Both c and c+1 are inside the connected component.
        case 1:
            Cell[4*x+2*y] = *(*(ptt3)+c-coo3->v1);
            Cell[4*x+2*y+1] = *(*(ptt3)+cinc-coo3->v1);
            return;
        // Only c is inside the connected component
        case 2:
            Cell[4*x+2*y] = *(*(ptt3)+c-coo3->v1);
            return;
        // Only c+1 is inside the connected component
        case 3:
            Cell[4*x+2*y+1] = *(*(ptt3)+cinc-coo3->v1);
            return;
        }
    }
    return;
}
 
 
template <typename T, typename K>
void DTGrid3D<T,K>::getCellyz(T * Cell, int x, int b, int c, typename ArrayPrs <K>::PairArray * coord2){
 
    short skey2; // a local integer variable to identify the grid cell type. 
    int binc = b+1; // the index of the origin's neighbor in y direction.
 
    K checkpoint_21 = coord2->v2; // the index of the first connected component in a tubular grid in the yCoord array.
    K checkpoint_22 = (coord2+1)->v2; // the index of the last connected component in a tubular grid in the yCoord array.
    K countdown_2 = checkpoint_22 - checkpoint_21; // Determine how many connected components are there in the corresponding tubular grid (2D).
 
    if (countdown_2 == 0){ return;}
 
    typename ArrayPrs <K>::PairArray ** ptt2 = proj2D->Ptr_List_2D->PtrArray+checkpoint_21; // a temporary pointer to store an acc array member.
    typename ArrayPrs <K>::PairArray * coo2 = proj2D->Ycoord_List->PArray+checkpoint_21; // a temporary pointer to store an yCoord array pair.
 
    // iterate over the connected components in DTGrid2D
    for (; countdown_2 >= 1; countdown_2--){
        // check whether index b (and/or b+1) is inside one of the connected components using the yCoord index array.
        skey2 = b>=coo2->v1 ? (binc<=coo2->v2 ? 1 : (b<=coo2->v2 ? 2 : 0) ) : (binc>=coo2->v1 ? (binc<=coo2->v2 ? 3 : 0) : 0);
        switch (skey2)
        {
        // Neither b, nor b+1 is inside the current connected component.
        case 0:
            ++ptt2;
            ++coo2;
            break;
        // Both b and b+1 are inside the connected component.           
        case 1:
            getCellz(Cell, x,0,c, *(ptt2)+b-coo2->v1);
            //++b;
            getCellz(Cell, x,1,c, *(ptt2)+binc-coo2->v1);
            return;
        case 2:
            getCellz(Cell, x,0,c, *(ptt2)+b-coo2->v1);
            return;
        case 3:
            getCellz(Cell, x,1,c, *(ptt2)+binc-coo2->v1);
            return;
        }
    }
    return;
}
 
 
template <typename T, typename K>
T * DTGrid3D<T,K>::getCell(T * Cell, int a, int b, int c){
    
    short skey1; // a local integer variable to identify the grid cell type. 
    int ainc = a+1; // the index of the origin's neighbor in x direction.
    int countdown_1 = proj2D->proj1D->set_countdown_1; // Determine how many connected components are there in the corresponding tubular grid (1D).
 
    typename ArrayPrs <K>::PairArray ** ptt1 = proj2D->proj1D->Ptr_List_1D->PtrArray; // a temporary pointer to store an acc array member.
    typename ArrayPrs <K>::PairArray * coo1 = proj2D->proj1D->Xcoord_List->PArray; // a temporary pointer to store an xCoord array pair.
 
    // iterate over the connected components in DTGrid1D
    for (;countdown_1>0;countdown_1--){
        // check whether index a (and/or a+1) is inside one of the connected components using the xCoord index array.
        skey1 = a>=coo1->v1 ? (ainc<=coo1->v2 ? 1 : (a<=coo1->v2 ? 2 : 0) ) : (ainc>=coo1->v1 ? (ainc<=coo1->v2 ? 3 : 0) : 0);
        switch (skey1)
        {
        // Neither a, nor a+1 is inside the connected component.
        case 0:
            ++ptt1;
            ++coo1;
            break;
        // Both a and a+1 are inside the connected component.       
        case 1:
            getCellyz(Cell, 0,b,c, *(ptt1) + a - coo1->v1);
            getCellyz(Cell, 1,b,c, *(ptt1) + ainc - coo1->v1);
            return Cell;
        // Only a is inside the connected component
        case 2:
            getCellyz(Cell, 0,b,c, *(ptt1) + a - coo1->v1);
            return Cell;
        // Only a+1 is inside the connected component
        case 3:
            getCellyz(Cell, 1,b,c, *(ptt1) + ainc - coo1->v1);
            return Cell;
        }
    }
    return Cell;
}
 
 
template <typename T, typename K>
T * DTGrid3D<T,K>::getCell2(T * Cell, int a, int b, int c){
 
    short skey1; // a local integer variable to identify the grid cell type. 
    int ainc = a+1; // the index of the origin's neighbor in x direction.
 
    typename ArrayPrs <K>::PairArray ** ptt1 = proj2D->proj1D->Ptr_List_1D->PtrArray; // a temporary pointer to store an acc array member.
    typename ArrayPrs <K>::PairArray * coo1 = proj2D->proj1D->Xcoord_List->PArray; // a temporary pointer to store an xCoord array pair.
    
    // check whether index a (and/or a+1) is inside one of the connected components using the xCoord index array.
    skey1 = a>=coo1->v1 ? (ainc<=coo1->v2 ? 1 : (a<=coo1->v2 ? 2 : 0) ) : (ainc>=coo1->v1 ? (ainc<=coo1->v2 ? 3 : 0) : 0);
    switch (skey1)
    {
    // Both a and a+1 are inside the connected component.
    case 1:
        getCellyz(Cell, 0,b,c, *(ptt1) + a - coo1->v1);
        getCellyz(Cell, 1,b,c, *(ptt1) + ainc - coo1->v1);
        return Cell;
    // Only a is inside the connected component
    case 2:
        getCellyz(Cell, 0,b,c, *(ptt1) + a - coo1->v1);
        return Cell;
    // Only a+1 is inside the connected component
    case 3:
        getCellyz(Cell, 1,b,c, *(ptt1) + ainc - coo1->v1);
        return Cell;
    default:
        return Cell;
    }
}
 
 
template <typename T, typename K>
void DTGrid3D<T,K>::test_format_dtgrid(char * Input)
{
    std::ifstream file_bt;
    file_bt.open(Input, std::ifstream::in | std::ifstream::binary);
 
    int flag; // flag variable, binary or ASCII?
    file_bt.read((char *) &flag, sizeof(int));
 
    char * buffer = new char [120]; // Read the first 120 characters
    file_bt.read(buffer, 120);
 
    char * binary_test; // check if buffer contains a return character
    binary_test = (char*) memchr (buffer, '\n', 120);
 
    if (flag == 160){
        iform = 0; // Binary file
        std::cout << "Test format DT-Grid: No return statement, binary zeros found" << std::endl;
        std::cout << "The file was detected to be of type binary" << std::endl;
    }
    else{
        if (binary_test == NULL) { // NO return statement found in the buffer
            binary_test = (char*) memchr (buffer, '\0', 120); // check if there is a null termination character.
            if (binary_test == NULL){
                iform = 1; // ASCII file
                std::cout << "Test format DT-Grid: No return statement, no binary zeros" << std::endl;
                std::cout << "The file was detected to be of type ASCII" << std::endl;
            }
            else{
                iform = 0; // Binary file
                std::cout << "Test format DT-Grid: No return statement, binary zeros found" << std::endl;
                std::cout << "The file was detected to be of type binary" << std::endl;
            }
        }
        else{
            iform = 1; // ASCII file
            std::cout << "Test format DT-Grid: Return statement found" << std::endl;
            std::cout << "The file was detected to be of type ASCII" << std::endl;
        }
    }
    file_bt.close();
    delete[] buffer;
}
 
 
template <typename T, typename K>
void DTGrid3D<T,K>::format_binary(char * file_binary)
{
    std::ifstream dtgrid_file;
    dtgrid_file.open(file_binary, std::ifstream::in | std::ifstream::binary);
 
    float sc; // scaling factor
    float sp; // grid spacing
    float ox, oy, oz; // origin coordinates
    float fval; // local variable to store float values temporarily
    int dimi, dimj, dimk; // grid dimensions
    int ival; // local variable to store int values temporarily
    int flag_binary; // the first integer stored in a dt-grid (or UHBD) binary file.
    int flag_grid_type; // excluded volume (0) / potential (1)
 
    // DT-Grid/UHBD header parameters
    //
    dtgrid_file.read((char *) &flag_binary, sizeof(int));
    uhbd.Title = new char [72];
    dtgrid_file.read(uhbd.Title, 72); // content description upto 72 characters.
    dtgrid_file.read((char *) &uhbd.scale, sizeof(float));
    dtgrid_file.read((char *) &uhbd.dummy1, sizeof(float)); // float dummy. value type not assigned.
    dtgrid_file.read((char *) &uhbd.grdflag, sizeof(int));
    dtgrid_file.read((char *) &uhbd.idummy1, sizeof(int)); // integer dummy. value type not assigned.
    dtgrid_file.read((char *) &uhbd.idummy2, sizeof(int));
    dtgrid_file.read((char *) &uhbd.one, sizeof(int));
    dtgrid_file.read((char *) &uhbd.idummy3, sizeof(int));
    dtgrid_file.read((char *) &uhbd.dim.i, sizeof(int));
    dtgrid_file.read((char *) &uhbd.dim.j, sizeof(int));
    dtgrid_file.read((char *) &uhbd.dim.k, sizeof(int));
    dtgrid_file.read((char *) &uhbd.spacing, sizeof(float));
    dtgrid_file.read((char *) &uhbd.o.x, sizeof(float));
    dtgrid_file.read((char *) &uhbd.o.y, sizeof(float));
    dtgrid_file.read((char *) &uhbd.o.z, sizeof(float));
    dtgrid_file.read((char *) &uhbd.dummy2, sizeof(float));
    dtgrid_file.read((char *) &uhbd.dummy3, sizeof(float));
    dtgrid_file.read((char *) &uhbd.dummy4, sizeof(float));
    dtgrid_file.read((char *) &uhbd.dummy5, sizeof(float));
    dtgrid_file.read((char *) &uhbd.dummy6, sizeof(float));
    dtgrid_file.read((char *) &uhbd.dummy7, sizeof(float));
    dtgrid_file.read((char *) &uhbd.idummy4, sizeof(int));
    dtgrid_file.read((char *) &uhbd.idummy5, sizeof(int));
 
    // DT-Grid header parameters
    // Length of each of the DT-Grid data structure arrays
    int len_Value, len_xcoord, len_ycoord, len_zcoord;
    int len_acc, len_acc2, len_acc3, len_projection, len_projection2;
 
    // Read the array size from DT-Grid file
    dtgrid_file.read((char *) &flag_grid_type, sizeof(int));
    dtgrid_file.read((char *) &len_Value, sizeof(int));
    dtgrid_file.read((char *) &len_xcoord, sizeof(int));
    dtgrid_file.read((char *) &len_ycoord, sizeof(int));
    dtgrid_file.read((char *) &len_zcoord, sizeof(int));
    dtgrid_file.read((char *) &len_acc, sizeof(int));
    dtgrid_file.read((char *) &len_acc2, sizeof(int));
    dtgrid_file.read((char *) &len_acc3, sizeof(int));
    dtgrid_file.read((char *) &len_projection, sizeof(int));
    dtgrid_file.read((char *) &len_projection2, sizeof(int));
 
    grid_type = flag_grid_type;
    if (grid_type == 0)
        std::cout << "Type of grid: Excluded volume" << std::endl;
    else if (grid_type == 1)
        std::cout << "Type of grid: Potential/Energy" << std::endl;
    else
        std::cout << "Type of grid: Unknown/Not-Set" << std::endl;
 
    // The following block verbose. Informative but necessary to print every time?
    std::cout << "DT-Grid Header information: " << std::endl;
    std::cout << "Length of Value array:\t" << len_Value << std::endl;
    std::cout << "Length of xCoord array:\t" << len_xcoord<< std::endl;
    std::cout << "Length of yCoord array:\t" << len_ycoord << std::endl;
    std::cout << "Length of zCoord array:\t" << len_zcoord << std::endl;
    std::cout << "Length of acc1D array:\t" << len_acc << std::endl;
    std::cout << "Length of acc2D array:\t" << len_acc2 << std::endl;
    std::cout << "Length of acc3D array:\t" << len_acc3 << std::endl;
    std::cout << "Length of proj1D array:\t" << len_projection << std::endl;
    std::cout << "Length of proj2D array:\t" << len_projection2 << std::endl << std::endl;
 
    // DT-Grid file size. Estimate from the array size given above.
    int dtgrid_size;
    dtgrid_size = sizeof(T) * len_Value + sizeof(K) * (len_xcoord + len_ycoord + len_zcoord + \
            len_acc + len_acc2 + len_acc3) + sizeof(K *) * (len_projection + len_projection2);
    size_ds(dtgrid_size); // convert to human-readable form
 
    //DT-Grid main body.
 
    // read the 'value' array from the file
    T * Value;
    if (grid_type == 0){
        // excluded volume grid
        Value = new T;
    }
    else if (grid_type == 1){
        // potential/energy grid
        Value = new T [len_Value];
        for (int i=0;i<len_Value;i++){
            dtgrid_file.read((char *) &fval, sizeof(float));
            if (typeid(T) == typeid(bool))
                Value[i] = fval;
            else
                Value[i] = (T) fval * scfct;
        }
    }
    else{
        // Grid types other than excluded volume and potential/energy. Perhaps in the future.
        Value = new T;
    }
 
    // Read the zcoord array from file
    K * zcoord;
    zcoord = new K [len_zcoord];
    for (int i=0;i<len_zcoord;i++){
        dtgrid_file.read((char *) &ival, sizeof(int));
        zcoord[i] = (K) ival;
    }
 
    // Read the proj1D acc array from file
    K * acc;
    acc = new K [len_acc];
    for (int i=0;i<len_acc;i++){
        dtgrid_file.read((char *) &ival, sizeof(int));
        acc[i] = (K) ival;
    }
 
    // Read the proj1D 'value' array
    K * projection;
    projection = new K [len_projection];
    for (int i=0;i<len_projection;i++){
        dtgrid_file.read((char *) &ival, sizeof(int));
        projection[i] = (K) ival;
    }
 
    // Read the ycoord array
    K * ycoord;
    ycoord = new K [len_ycoord];
    for (int i=0;i<len_ycoord;i++){
        dtgrid_file.read((char *) &ival, sizeof(int));
        ycoord[i] = (K) ival;
    }
 
    // Read the proj2D acc array
    K * acc2;
    acc2 = new K [len_acc2];
    for (int i=0;i<len_acc2;i++){
        dtgrid_file.read((char *) &ival, sizeof(int));
        acc2[i] = (K) ival;
    }
 
    // Read the proj2D 'value' array
    K * projection2;
    projection2 = new K [len_projection2];
    for (int i=0;i<len_projection2;i++){
        dtgrid_file.read((char *) &ival, sizeof(int));
        projection2[i] = (K) ival;
    }
 
    // Read the xcoord array
    K * xcoord;
    xcoord = new K [len_xcoord];
    for (int i=0;i<len_xcoord;i++){
        dtgrid_file.read((char *) &ival, sizeof(int));
        xcoord[i] = (K) ival;
    }
 
    // Read the acc array that points into DTGrid3D 'value' array
    K * acc3;
    acc3 = new K [len_acc3];
    for (int i=0;i<len_acc3;i++){
        dtgrid_file.read((char *) &ival, sizeof(int));
        acc3[i] = (K) ival;
    }
 
    // Initialize the DT-Grid data structure arrays.
    set_uhbd_parameters(uhbd.scale, uhbd.spacing, uhbd.o.x, uhbd.o.y, uhbd.o.z, uhbd.dim.i, uhbd.dim.j, uhbd.dim.k);
    set_value_list_3D(Value,len_Value);
    set_zcoord_list(zcoord,len_zcoord);
    set_ptr_list_3D(acc,len_acc);
    set_projection2D(projection,ycoord,acc2,projection2,xcoord,acc3,len_projection, \
            len_ycoord,len_acc2,len_projection2,len_xcoord,len_acc3);
 
    // deallocate the arrays generated for temporarily use.
    if (grid_type == 0){
        delete Value;
    }
    else if (grid_type == 1){
        delete[] Value;
    }
    else{
        delete Value;
    }
    delete[] zcoord;
    delete[] acc;
    delete[] projection;
    delete[] ycoord;
    delete[] acc2;
    delete[] projection2;
    delete[] xcoord;
    delete[] acc3;
 
    dtgrid_file.close();
}
 
 
template <typename T, typename K>
void DTGrid3D<T,K>::format_ascii(char * file_ascii)
{
    std::ifstream dtgrid_file;
    dtgrid_file.open(file_ascii, std::ifstream::in);
 
    std::string title_str, line;
    std::string scale_str, dum1_str, grdflag_str, idum1_str, idum2_str, one_str, idum3_str;
    std::string i_str, j_str, k_str;
    std::string spacing_str, ox_str, oy_str, oz_str;
    std::string dum2_str, dum3_str, dum4_str, dum5_str;
    std::string dum6_str, dum7_str, idum4_str, idum5_str;
 
    T fval;
    K ival;
    int flag_grid_type;
    int len_Value, len_xcoord, len_ycoord, len_zcoord;
    int len_acc, len_acc2, len_acc3, len_projection, len_projection2;
 
    // DT-Grid (UHBD) Title
    getline(dtgrid_file,line);
    title_str = line.substr(0,72);
 
    // DT-Grid/UHBD header parameters
    //
    getline(dtgrid_file,line);
    scale_str = line.substr(0,12);
    dum1_str = line.substr(12,12);
    grdflag_str = line.substr(24,7);
    idum1_str = line.substr(31,7);
    idum2_str = line.substr(38,7);
    one_str = line.substr(45,7);
    idum3_str = line.substr(52,7);
    getline(dtgrid_file,line);
    i_str = line.substr(0,7);
    j_str = line.substr(7,7);
    k_str = line.substr(14,7);
    spacing_str = line.substr(21,12);
    ox_str = line.substr(33,12);
    oy_str = line.substr(45,12);
    oz_str = line.substr(57,12);
    getline(dtgrid_file,line);
    dum2_str = line.substr(0,12);
    dum3_str = line.substr(12,12);
    dum4_str = line.substr(24,12);
    dum5_str = line.substr(36,12);
    getline(dtgrid_file,line);
    dum6_str = line.substr(0,12);
    dum7_str = line.substr(12,12);
    idum4_str = line.substr(24,7);
    idum5_str = line.substr(31,7);
 
    // assign the grid parameters from the UHBD file header
    uhbd.Title = new char [72];
    std::strncpy (uhbd.Title, title_str.c_str(), 72);
 
    std::stringstream (scale_str) >> uhbd.scale;
    std::stringstream (dum1_str) >> uhbd.dummy1;
    std::stringstream (grdflag_str) >> uhbd.grdflag;
    std::stringstream (idum1_str) >> uhbd.idummy1;
    std::stringstream (idum2_str) >> uhbd.idummy2;
    std::stringstream (one_str) >> uhbd.one;
    std::stringstream (idum3_str) >> uhbd.idummy3;
    std::stringstream (i_str) >> uhbd.dim.i;
    std::stringstream (j_str) >> uhbd.dim.j;
    std::stringstream (k_str) >> uhbd.dim.k;
    std::stringstream (spacing_str) >> uhbd.spacing;
    std::stringstream (ox_str) >> uhbd.o.x;
    std::stringstream (oy_str) >> uhbd.o.y;
    std::stringstream (oz_str) >> uhbd.o.z;
    std::stringstream (dum2_str) >> uhbd.dummy2;
    std::stringstream (dum3_str) >> uhbd.dummy3;
    std::stringstream (dum4_str) >> uhbd.dummy4;
    std::stringstream (dum5_str) >> uhbd.dummy5;
    std::stringstream (dum6_str) >> uhbd.dummy6;
    std::stringstream (dum7_str) >> uhbd.dummy7;
    std::stringstream (idum4_str) >> uhbd.idummy4;
    std::stringstream (idum5_str) >> uhbd.idummy5;
 
 
    // DT-Grid header parameters
    // Length of each of the DT-Grid data structure arrays
    dtgrid_file >> flag_grid_type;
    dtgrid_file >> len_Value;
    dtgrid_file >> len_xcoord;
    dtgrid_file >> len_ycoord;
    dtgrid_file >> len_zcoord;
    dtgrid_file >> len_acc;
    dtgrid_file >> len_acc2;
    dtgrid_file >> len_acc3;
    dtgrid_file >> len_projection;
    dtgrid_file >> len_projection2;
 
    grid_type = flag_grid_type;
    if (grid_type == 0)
        std::cout << "Type of grid: Excluded volume" << std::endl;
    else if (grid_type == 1)
        std::cout << "Type of grid: Potential/Energy" << std::endl;
    else
        std::cout << "Type of grid: Unknown/Not-Set" << std::endl;
 
    // The following block verbose. Informative but necessary to print every time?
    std::cout << "DT-Grid Header information: " << std::endl;
    std::cout << "Length of Value array:\t" << len_Value << std::endl;
    std::cout << "Length of xCoord array:\t" << len_xcoord<< std::endl;
    std::cout << "Length of yCoord array:\t" << len_ycoord << std::endl;
    std::cout << "Length of zCoord array:\t" << len_zcoord << std::endl;
    std::cout << "Length of acc1D array:\t" << len_acc << std::endl;
    std::cout << "Length of acc2D array:\t" << len_acc2 << std::endl;
    std::cout << "Length of acc3D array:\t" << len_acc3 << std::endl;
    std::cout << "Length of proj1D array:\t" << len_projection << std::endl;
    std::cout << "Length of proj2D array:\t" << len_projection2 << std::endl << std::endl;
 
    unsigned int dtgrid_size;
    dtgrid_size = sizeof(T) * len_Value + sizeof(K) * (len_xcoord + len_ycoord + len_zcoord + \
            len_acc + len_acc2 + len_acc3) + sizeof(K *) * (len_projection + len_projection2);
    size_ds(dtgrid_size); // convert to human-readable form
 
    // DT-Grid file main body
    //
    T * Value;
    if (grid_type == 0){
        // excluded volume grid
        Value = new T;
    }
    else if (grid_type == 1){
        // potential grid
        Value = new T [len_Value];
        for (int i=0;i<len_Value;i++){
            dtgrid_file >> fval;
            if (typeid(T) == typeid(bool))
                Value[i] = fval;
            else
                Value[i] = fval * scfct;
        }
    }
    else{
        // Grid types other than excluded volume and potential/energy. Perhaps in the future.
        Value = new T;
    }
 
    // Read the zcoord array
    K * zcoord;
    zcoord = new K [len_zcoord];
    for (int i=0;i<len_zcoord;i++){
        dtgrid_file >> ival;
        zcoord[i] = ival;
    }
 
    // Read the proj1D acc array
    K * acc;
    acc = new K [len_acc];
    for (int i=0;i<len_acc;i++){
        dtgrid_file >> ival;
        acc[i] = ival;
    }
 
    // Read the proj1D 'value' array
    K * projection;
    projection = new K [len_projection];
    for (int i=0;i<len_projection;i++){
        dtgrid_file >> ival;
        projection[i] = ival;
    }
 
    // Read the ycoord array
    K * ycoord;
    ycoord = new K [len_ycoord];
    for (int i=0;i<len_ycoord;i++){
        dtgrid_file >> ival;
        ycoord[i] = ival;
    }
 
    // Read the proj2D acc array
    K * acc2;
    acc2 = new K [len_acc2];
    for (int i=0;i<len_acc2;i++){
        dtgrid_file >> ival;
        acc2[i] = ival;
    }
    // Read the proj2D 'value' array
    K * projection2;
    projection2 = new K [len_projection2];
    for (int i=0;i<len_projection2;i++){
        dtgrid_file >> ival;
        projection2[i] = ival;
    }
 
    // Read the xcoord array
    K * xcoord;
    xcoord = new K [len_xcoord];
    for (int i=0;i<len_xcoord;i++){
        dtgrid_file >> ival;
        xcoord[i] = ival;
    }
 
    // Read the acc array that points into DTGrid3D 'value' array
    K * acc3;
    acc3 = new K [len_acc3];
    for (int i=0;i<len_acc3;i++){
        dtgrid_file >> ival;
        acc3[i] = ival;
    }
 
    // Initialize the DT-Grid data structure arrays.
    set_uhbd_parameters(uhbd.scale, uhbd.spacing, uhbd.o.x, uhbd.o.y, uhbd.o.z, uhbd.dim.i, uhbd.dim.j, uhbd.dim.k);
    set_value_list_3D(Value,len_Value);
    set_zcoord_list(zcoord,len_zcoord);
    set_ptr_list_3D(acc,len_acc);
    set_projection2D(projection,ycoord,acc2,projection2,xcoord,acc3,len_projection, \
            len_ycoord,len_acc2,len_projection2,len_xcoord,len_acc3);
 
    if (grid_type == 0){
        delete Value;
    }
    else if (grid_type == 1){
        delete[] Value;
    }
    else{
        delete Value;
    }
    delete[] zcoord;
    delete[] acc;
    delete[] projection;
    delete[] ycoord;
    delete[] acc2;
    delete[] projection2;
    delete[] xcoord;
    delete[] acc3;
 
    dtgrid_file.close();
}
 
 
template <typename T, typename K>
void DTGrid3D<T,K>::size_ds(unsigned int val)
{
    float dtgrid_size;
    int kb, mb, gb;
 
    kb = 1024; // kilobyte
    mb = 1024 * 1024; // megabyte
    gb = 1024 * 1024 * 1024; // gigabyte
 
    if (((float)val / (float) gb) >= 1.0) {
        dtgrid_size = (float) val / (float) gb;
        std::cout << "Size of DT-Grid data structure is " << dtgrid_size << " GB" << std::endl;
    }
    else {
        if (((float)val / (float) mb) >= 1.0) {
            dtgrid_size = (float) val / (float) mb;
            std::cout << "Size of DT-Grid data structure is " << dtgrid_size << " MB" << std::endl;
        }
        else {
            dtgrid_size = (float) val / (float) kb;
            std::cout << "Size of DT-Grid data structure is " << dtgrid_size << " KB" << std::endl;
        }
    }
}
 
#endif /* DTGRID3D_HPP_ */