%fname = 'h-bond2.cube/h-bond2.cube';

if (exist('fname', 'var') == 0)
    fname = input('Enter file name: ', "s")
end

fid = fopen(fname);
hdr1 = fgetl(fid);
hdr2 = fgetl(fid);

c = clock;
['Start time: ', num2str(c(4)), ':', num2str(c(5))]

%First line with actual useful content
%Number of atoms and the start coordinate of the info
L1 = fgetl(fid);
L1spl = split(L1);
natoms = str2num(L1spl{2});
StartPosX = str2num(L1spl{3});
StartPosY = str2num(L1spl{4});
StartPosZ = str2num(L1spl{5});

%Next 3 lines define the box with information
%NVoxAx1, Ax1X, Ax1Y and Ax1Z
L2 = fgetl(fid);
L2spl = split(L2);
NVox1 = str2num(L2spl{2});
Ax1X = str2num(L2spl{3});
Ax1Y = str2num(L2spl{4});
Ax1Z = str2num(L2spl{5});

%NVoxAx2, Ax2X, Ax2Y and Ax2Z
L3 = fgetl(fid);
L3spl = split(L3);
NVox2 = str2num(L3spl{2});
Ax2X = str2num(L3spl{3});
Ax2Y = str2num(L3spl{4});
Ax2Z = str2num(L3spl{5});

%NVoxAx1, Ax1X, Ax1Y and Ax1Z
L4 = fgetl(fid);
L4spl = split(L4);
NVox3 = str2num(L4spl{2});
Ax3X = str2num(L4spl{3});
Ax3Y = str2num(L4spl{4});
Ax3Z = str2num(L4spl{5});

%Now for a for loop - Atom number, charge, and x, y, z
AtNum = zeros(1, natoms);
AtQ = zeros(1, natoms);
AtPos = zeros(3, natoms);
for cnt = 1:natoms
    %read in line
    Lat = fgetl(fid);
    Latsp = split(Lat);
    AtNum(cnt) = str2num(Latsp{2});
    AtQ(cnt) = str2num(Latsp{3});
    AtPos(1, cnt) = str2num(Latsp{4});
    AtPos(2, cnt) = str2num(Latsp{5});
    AtPos(3, cnt) = str2num(Latsp{6});
end

cntx = 1;
cnty = 1;
cntz = 1;

f = waitbar(0, 'Loading data');

Vox = zeros(NVox1, NVox2, NVox3);
Lvox = fgetl(fid);
while(Lvox ~= -1)
    Lvoxsp = split(Lvox);
    nvoxels = length(Lvoxsp)-2;
    for cnt = 2:(nvoxels+1)
        Vox(cntx, cnty, cntz) = str2num(Lvoxsp{cnt});
        cntz = cntz + 1;
        if (cntz > NVox3)
            cntz = 1;
            cnty = cnty + 1;
            if (cnty > NVox2)
                cnty = 1;
                cntx = cntx + 1;
                waitbar(cntx/NVox1, f);
            end
        end
    end
    Lvox = fgetl(fid);
end
close(f)
fclose(fid);

c = clock;
['End time: ', num2str(c(4)), ':', num2str(c(5))]