Commit d058d28e authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Fix transition matrix element output to binary format

parent 086aba2a

src/libnptm/TransitionMatrix.cpp

+50 −20
Original line number Diff line number Diff line
@@ -110,22 +110,22 @@ TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) { 
    // This vector will be passed to the new object. DO NOT DELETE HERE!

    complex<double> *_elements;

    double _radius = 0.0;

    status = hdf_file->read("/IS", "INT32", &_is);

    status = hdf_file->read("/L_MAX", "INT32", &_lm);

    status = hdf_file->read("/VK", "FLOAT64", &_vk);

    status = hdf_file->read("/EXRI", "FLOAT64", &_exri);

    status = hdf_file->read("IS", "INT32", &_is);

    status = hdf_file->read("L_MAX", "INT32", &_lm);

    status = hdf_file->read("VK", "FLOAT64", &_vk);

    status = hdf_file->read("EXRI", "FLOAT64", &_exri);

    if (_is == 1111) {

      int num_elements = 2 * _lm;

      double *file_vector = new double[2 * num_elements]();

      hid_t file_id = hdf_file->get_file_id();

      hid_t dset_id = H5Dopen2(file_id, "/ELEMENTS", H5P_DEFAULT);

      hid_t dset_id = H5Dopen2(file_id, "ELEMENTS", H5P_DEFAULT);

      status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, file_vector);

      _elements = new complex<double>[num_elements]();

      for (int ei = 0; ei < num_elements; ei++) {

	_elements[ei] = complex<double>(file_vector[2 * ei], file_vector[2 * ei + 1]);

      }

      status = H5Dclose(dset_id);

      status = hdf_file->read("/RADIUS", "FLOAT64", &_radius);

      status = hdf_file->read("RADIUS", "FLOAT64", &_radius);

      tm = new TransitionMatrix(_is, _lm, _vk, _exri, _elements, _radius);

      delete[] file_vector;

    } else if (_is == 1) {
@@ -133,14 +133,14 @@ TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) { 
      int num_elements = nlemt * nlemt;

      double *file_vector = new double[2 * num_elements]();

      hid_t file_id = hdf_file->get_file_id();

      hid_t dset_id = H5Dopen2(file_id, "/ELEMENTS", H5P_DEFAULT);

      hid_t dset_id = H5Dopen2(file_id, "ELEMENTS", H5P_DEFAULT);

      status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, file_vector);

      _elements = new complex<double>[num_elements]();

      for (int ei = 0; ei < num_elements; ei++) {

	_elements[ei] = complex<double>(file_vector[2 * ei], file_vector[2 * ei + 1]);

      }

      status = H5Dclose(dset_id);

      status = hdf_file->read("/RADIUS", "FLOAT64", &_radius);

      status = hdf_file->read("RADIUS", "FLOAT64", &_radius);

      tm = new TransitionMatrix(_is, _lm, _vk, _exri, _elements, _radius);

      delete[] file_vector;

    }
@@ -215,33 +215,32 @@ void TransitionMatrix::write_hdf5(string file_name) { 
    List<string> rec_type_list(1);

    List<void *> rec_ptr_list(1);

    string str_type, str_name;

    rec_name_list.set(0, "/IS");

    rec_name_list.set(0, "IS");

    rec_type_list.set(0, "INT32_(1)");

    rec_ptr_list.set(0, &is);

    rec_name_list.append("/L_MAX");

    rec_name_list.append("L_MAX");

    rec_type_list.append("INT32_(1)");

    rec_ptr_list.append(&l_max);

    rec_name_list.append("/VK");

    rec_name_list.append("VK");

    rec_type_list.append("FLOAT64_(1)");

    rec_ptr_list.append(&vk);

    rec_name_list.append("/EXRI");

    rec_name_list.append("EXRI");

    rec_type_list.append("FLOAT64_(1)");

    rec_ptr_list.append(&exri);

    rec_name_list.append("/ELEMENTS");

    rec_name_list.append("ELEMENTS");

    str_type = "FLOAT64_(" + to_string(shape[0]) + "," + to_string(2 * shape[1]) + ")";

    rec_type_list.append(str_type);

    // The (N x M) matrix of complex is converted to a (N x 2M) matrix of double,

    // where REAL(E_i,j) -> E_i,(2 j) and IMAG(E_i,j) -> E_i,(2 j + 1)

    double *ptr_elements = new double[shape[0] * 2 * shape[1]]();

    for (int ei = 0; ei < shape[0]; ei++) {

      for (int ej = 0; ej < shape[1]; ej++) {

	ptr_elements[shape[1] * ei + 2 * ej] = elements[shape[1] * ei + ej].real();

	ptr_elements[shape[1] * ei + 2 * ej + 1] = elements[shape[1] * ei  + ej].imag();

      }

    int num_elements = 2 * shape[0] * shape[1];

    double *ptr_elements = new double[num_elements]();

    for (int ei = 0; ei < num_elements / 2; ei++) {

      ptr_elements[2 * ei] = elements[ei].real();

      ptr_elements[2 * ei + 1] = elements[ei].imag();

    }

    rec_ptr_list.append(ptr_elements);

    if (is == 1111) {

      rec_name_list.append("/RADIUS");

      rec_name_list.append("RADIUS");

      rec_type_list.append("FLOAT64_(1)");

      rec_ptr_list.append(&sphere_radius);

    }
@@ -299,3 +298,34 @@ void TransitionMatrix::write_legacy(string file_name) { 
    printf("ERROR: could not open Transition Matrix file for writing.\n");

  }

}



bool TransitionMatrix::operator ==(TransitionMatrix &other) {

  if (is != other.is) {

    return false;

  }

  if (l_max != other.l_max) {

    return false;

  }

  if (vk != other.vk) {

    return false;

  }

  if (exri != other.exri) {

    return false;

  }

  if (sphere_radius != other.sphere_radius) {

    return false;

  }

  if (shape[0] != other.shape[0]) {

    return false;

  }

  if (shape[1] != other.shape[1]) {

    return false;

  }

  int num_elements = shape[0] * shape[1];

  for (int ei = 0; ei < num_elements; ei++) {

    if (elements[ei] != other.elements[ei]) {

      return false;

    }

  }

  return true;

}