From b2943be7c1cf3654f04104eb0d1336a30ecfb2bd Mon Sep 17 00:00:00 2001
From: Yohai Meiron <sahmes@gmail.com>
Date: Sat, 28 Mar 2020 01:34:37 -0400
Subject: [PATCH] Improved HDF5 reader and added ASCII writer with custom
 precision

---
 io.cpp | 48 +++++++++++++++++++++++++++++++++---------------
 1 file changed, 33 insertions(+), 15 deletions(-)

diff --git a/io.cpp b/io.cpp
index b1286af..dce1d55 100644
--- a/io.cpp
+++ b/io.cpp
@@ -1,4 +1,5 @@
 #include "hdf5.h"
+#include <cmath>
 #include <string>
 #include <vector>
 #include <algorithm>
@@ -6,6 +7,18 @@
 struct double3 {double x,y,z;
 };
 
+void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const double *m, const double3 *x, const double3 *v, int precision=10)
+{
+    int id_width = (int)log10(N-1) + 1;
+    char string_template[256];
+    sprintf(string_template, "%%0%dd%%%d.%dE%%%d.%dE%%%d.%dE%%%d.%dE%%%d.%dE%%%d.%dE%%%d.%dE\n", id_width, precision+7, precision, precision+8, precision, precision+8, precision, precision+8, precision, precision+8, precision, precision+8, precision, precision+8, precision);
+    printf("%d\n", step_num);
+    printf("%d\n", N);
+    printf("%.16E\n", t);
+    for (int i=0; i<N; i++) {
+        printf(string_template, i, m[i], x[i].x, x[i].y, x[i].z, v[i].x, v[i].y, v[i].z);
+    }
+}
 
 void h5_write(const std::string file_name, const int step_num, const int N, const double t, const double *m, const double3 *x, const double3 *v, const double *pot, const double3 *acc, const double3 *jrk, const int write_mode=0, const bool use_double_precision=true)
 {
@@ -70,13 +83,17 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub
     *step_num = *max_element(step_num_arr.begin(), step_num_arr.end());
 
     // Prepare to read the data sets dimensionality data.
-    char dataset_path[256];
-    hid_t dataset_id, dataspace_id;
+    char path[256];
+    hid_t attr_id, dataset_id, dataspace_id;
     int ndims;
     hsize_t dims[2];
 
-    sprintf(dataset_path, "/Step#%d/MASS", *step_num);
-    dataset_id = H5Dopen2(file_id, dataset_path, H5P_DEFAULT);
+    sprintf(path, "/Step#%d", *step_num);
+    attr_id = H5Aopen_by_name(file_id, path, "Time", H5P_DEFAULT, H5P_DEFAULT);
+    H5Aread(attr_id, H5T_NATIVE_DOUBLE, t);
+
+    sprintf(path, "/Step#%d/MASS", *step_num);
+    dataset_id = H5Dopen2(file_id, path, H5P_DEFAULT);
     dataspace_id = H5Dget_space(dataset_id);
     ndims = H5Sget_simple_extent_ndims(dataspace_id);
     if (ndims != 1)
@@ -85,24 +102,22 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub
     H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, m);
     *N = dims[0];
 
-    sprintf(dataset_path, "/Step#%d/X", *step_num);
-    dataset_id = H5Dopen2(file_id, dataset_path, H5P_DEFAULT);
+    sprintf(path, "/Step#%d/X", *step_num);
+    dataset_id = H5Dopen2(file_id, path, H5P_DEFAULT);
     dataspace_id = H5Dget_space(dataset_id);
     ndims = H5Sget_simple_extent_ndims(dataspace_id);
-    if (ndims != 2) throw std::runtime_error("Dataset X in Step#" + std::to_string(*step_num) + " of file " + file_name + " is " + std::to_string(ndims) + "-dimensional (expected 2-dimensional)");
+    if (ndims != 2) throw std::runtime_error("Bad dimensionality");
     H5Sget_simple_extent_dims(dataspace_id, dims, NULL);
-    if (dims[0] != *N) throw std::runtime_error("Expected " + std::to_string(*N) + "rows in Dataset X in Step#" + std::to_string(*step_num) + " of file " + file_name + " (based on MASS)");
-    if (dims[1] != 3) throw std::runtime_error("Dataset X in Step#" + std::to_string(*step_num) + " of file " + file_name + " has " + std::to_string(dims[1]) + " columns (expected 3 columns)");
+    if ((dims[0] != *N) || (dims[1] != 3))  throw std::runtime_error("Bad dimensionality");
     H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (double*)x);
 
-    sprintf(dataset_path, "/Step#%d/V", *step_num);
-    dataset_id = H5Dopen2(file_id, dataset_path, H5P_DEFAULT);
+    sprintf(path, "/Step#%d/V", *step_num);
+    dataset_id = H5Dopen2(file_id, path, H5P_DEFAULT);
     dataspace_id = H5Dget_space(dataset_id);
     ndims = H5Sget_simple_extent_ndims(dataspace_id);
-    if (ndims != 2) throw std::runtime_error("Dataset V in Step#" + std::to_string(*step_num) + " of file " + file_name + " is " + std::to_string(ndims) + "-dimensional (expected 2-dimensional)");
+    if (ndims != 2) throw std::runtime_error("Bad dimensionality");
     H5Sget_simple_extent_dims(dataspace_id, dims, NULL);
-    if (dims[0] != *N) throw std::runtime_error("Expected " + std::to_string(*N) + "rows in Dataset V in Step#" + std::to_string(*step_num) + " of file " + file_name + " (based on MASS)");
-    if (dims[1] != 3) throw std::runtime_error("Dataset V in Step#" + std::to_string(*step_num) + " of file " + file_name + " has " + std::to_string(dims[1]) + " columns (expected 3 columns)");
+    if ((dims[0] != *N) || (dims[1] != 3)) throw std::runtime_error("Bad dimensionality");
     H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (double*)v);
 }
 
@@ -127,7 +142,10 @@ int main()
     h5_write("123.h5", 17, N, 123.987, (double*)x, x, x, (double*)x, x, x, 7, true);
     h5_read("123.h5", &step_num, &N, &t, m, x, v);
     printf("%d %d %f\n", step_num, N, t);
-#error still havent implemented reading the time attribute
+    
+    
+    
+    ascii_write("123.dat", 0, N, t, m, x, v, 10);
     
 //     h5_write(file_name, 34, 987.654, N, (double*)x, x, x, (double*)x, x, x, 7);