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

Use HDF5 format to read the T-matrix

parent eb47cd0e

src/trapping/clffft.cpp

+347 −327
Original line number Original line Diff line number Diff line
@@ -23,6 +23,7 @@ 
#include <cstdio>

#include <cstdio>

#include <exception>

#include <exception>

#include <fstream>

#include <fstream>

#include <hdf5.h>

#include <regex>

#include <regex>

#include <string>

#include <string>
@@ -58,6 +59,22 @@ 
#include "../include/tra_subs.h"

#include "../include/tra_subs.h"

#endif

#endif





#ifndef INCLUDE_ERRORS_H_

#include "../include/errors.h"

#endif



#ifndef INCLUDE_LIST_H_

#include "../include/List.h"

#endif



#ifndef INCLUDE_TRANSITIONMATRIX_H_

#include "../include/TransitionMatrix.h"

#endif



#ifndef INCLUDE_FILE_IO_H_

#include "../include/file_io.h"

#endif



#ifdef USE_NVTX

#ifdef USE_NVTX

#include <nvtx3/nvToolsExt.h>

#include <nvtx3/nvToolsExt.h>

#endif

#endif
@@ -85,6 +102,7 @@ void lffft(string data_file, string output_path) { 
  fstream tlfff, tlfft;

  fstream tlfff, tlfft;

  double ****zpv = NULL;

  double ****zpv = NULL;

  dcomplex *ac = NULL, *ws = NULL, *wsl = NULL;

  dcomplex *ac = NULL, *ws = NULL, *wsl = NULL;

  dcomplex *vec_am0m = NULL;

  dcomplex **am0m = NULL;

  dcomplex **am0m = NULL;

  dcomplex **amd = NULL;

  dcomplex **amd = NULL;

  int **indam = NULL;

  int **indam = NULL;
@@ -107,83 +125,87 @@ void lffft(string data_file, string output_path) { 
    else if (mi == 2) jtw = stoi(m.str());

    else if (mi == 2) jtw = stoi(m.str());

    str_target = m.suffix().str();

    str_target = m.suffix().str();

  } // mi loop

  } // mi loop

  string ttms_name = output_path + "/c_TTMS";

  string ttms_name = output_path + "/c_TTMS.hd5";

  fstream ttms;

  // fstream ttms;

  ttms.open(ttms_name, ios::in | ios::binary);

  // ttms.open(ttms_name, ios::in | ios::binary);

  if (ttms.is_open()) {

  // if (ttms.is_open()) {

    ttms.read(reinterpret_cast<char *>(&is), sizeof(int));

  TransitionMatrix *ttms = TransitionMatrix::from_binary(ttms_name, "HDF5");

    ttms.read(reinterpret_cast<char *>(&le), sizeof(int));

  if (ttms == NULL)

    ttms.read(reinterpret_cast<char *>(&vks), sizeof(double));

    throw(ObjectAllocationException("ERROR: could not load the T-matrix file!\n"));

    ttms.read(reinterpret_cast<char *>(&exris), sizeof(double));

  // ttms.read(reinterpret_cast<char *>(&is), sizeof(int));

  // ttms.read(reinterpret_cast<char *>(&le), sizeof(int));

  // ttms.read(reinterpret_cast<char *>(&vks), sizeof(double));

  // ttms.read(reinterpret_cast<char *>(&exris), sizeof(double));

  is = ttms->is;

  le = ttms->l_max;

  vks = ttms->vk;

  exris = ttms->exri;

  cil->le = le;

  cil->le = le;

  cil->nlem = le * (le + 2);

  cil->nlem = le * (le + 2);

  cil->nlemt = cil->nlem + cil->nlem;

  cil->nlemt = cil->nlem + cil->nlem;

  if (is >= 2222) { // label 120

  if (is >= 2222) { // label 120

      tms = new dcomplex*[le];

    // READING OF T-MATRIX WITH IS = 2222 IS NOT SUPPORTED!

      for (int ti = 0; ti < le; ti++) tms[ti] = new dcomplex[3]();

    // tms = new dcomplex*[le];

      int lm = le;

    // for (int ti = 0; ti < le; ti++) tms[ti] = new dcomplex[3]();

      for (int i = 0; i < lm; i++) {

    // int lm = le;

	double vreal, vimag;

    // for (int i = 0; i < lm; i++) {

	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));

    //   double vreal, vimag;

	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));

    //   ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));

	tms[i][0] = vreal + vimag * I;

    //   ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));

	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));

    //   tms[i][0] = vreal + vimag * I;

	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));

    //   ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));

	tms[i][1] = vreal + vimag * I;

    //   ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));

	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));

    //   tms[i][1] = vreal + vimag * I;

	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));

    //   ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));

	tms[i][2] = vreal + vimag * I;

    //   ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));

      } // i loop

    //   tms[i][2] = vreal + vimag * I;

    //} // i loop

    throw(UnrecognizedConfigurationException("ERROR: T-matrix with IS>=2222 not supported!\n"));

  } else if (is >= 1111) { // label 125

  } else if (is >= 1111) { // label 125

    tmsm = new dcomplex[le]();

    tmsm = new dcomplex[le]();

    tmse = new dcomplex[le]();

    tmse = new dcomplex[le]();

    for (int i = 0; i < le; i++) {

    for (int i = 0; i < le; i++) {

	double vreal, vimag;

      tmsm[i] = ttms->elements[2 * i];

	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));

      tmse[i] = ttms->elements[2 * i + 1];

	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));

	tmsm[i] = vreal + vimag * I;

	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));

	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));

	tmse[i] = vreal + vimag * I;

    } // i loop

    } // i loop

  } else if (is >= 0) { // label 135

  } else if (is >= 0) { // label 135

    vec_am0m = new dcomplex[cil->nlemt * cil->nlemt];

    am0m = new dcomplex*[cil->nlemt];

    am0m = new dcomplex*[cil->nlemt];

      for (int ai = 0; ai < cil->nlemt; ai++) am0m[ai] = new dcomplex[cil->nlemt]();

    for (int ai = 0; ai < cil->nlemt; ai++) am0m[ai] = vec_am0m + cil->nlemt * ai;

    for (int i = 0; i < cil->nlemt; i++) {

    for (int i = 0; i < cil->nlemt; i++) {

      for (int j = 0; j < cil->nlemt; j++) {

      for (int j = 0; j < cil->nlemt; j++) {

	  double vreal, vimag;

	am0m[i][j] = ttms->elements[cil->nlemt * i + j];

	  ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));

	  ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));

	  am0m[i][j] = vreal + vimag * I;

      } // j loop

      } // j loop

    } // i loop

    } // i loop

  } else if (is < 0) {

  } else if (is < 0) {

      nvam = le * le + (le * (le + 1) * (le * 2 + 1)) / 3;

    // nvam = le * le + (le * (le + 1) * (le * 2 + 1)) / 3;

      amd = new dcomplex*[nvam];

    // amd = new dcomplex*[nvam];

      for (int ai = 0; ai < nvam; ai++) amd[ai] = new dcomplex[4]();

    // for (int ai = 0; ai < nvam; ai++) amd[ai] = new dcomplex[4]();

      for (int i = 0; i < nvam; i++) {

    // for (int i = 0; i < nvam; i++) {

	for (int j = 0; j < 4; j++) {

    //   for (int j = 0; j < 4; j++) {

	  double vreal, vimag;

    // 	double vreal, vimag;

	  ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));

    // 	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));

	  ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));

    // 	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));

	  amd[i][j] = vreal + vimag * I;

    // 	amd[i][j] = vreal + vimag * I;

	} // j loop

    //   } // j loop

      } // i loop

    // } // i loop

      indam = new int*[le];

    // indam = new int*[le];

      int vint;

    // int vint;

      for (int ii = 0; ii < le; ii++) indam[ii] = new int[le]();

    // for (int ii = 0; ii < le; ii++) indam[ii] = new int[le]();

      for (int i = 0; i < le; i++) {

    // for (int i = 0; i < le; i++) {

	for (int j = 0; j < le; j++) {

    //   for (int j = 0; j < le; j++) {

	  ttms.read(reinterpret_cast<char *>(&vint), sizeof(int));

    // 	ttms.read(reinterpret_cast<char *>(&vint), sizeof(int));

	  indam[i][j] = vint;

    // 	indam[i][j] = vint;

	} // j loop

    //   } // j loop

      } // i loop

    // } // i loop

      ttms.read(reinterpret_cast<char *>(&vint), sizeof(int));

    // ttms.read(reinterpret_cast<char *>(&vint), sizeof(int));

      cil->mxmpo = vint;

    // cil->mxmpo = vint;

      cil->mxim = vint * 2 - 1;

    // cil->mxim = vint * 2 - 1;

    throw(UnrecognizedConfigurationException("ERROR: T-matrix with IS<0 not supported!\n"));

  }

  }

  // label 150

  // label 150

    ttms.close();

  // ttms.close();

  delete ttms;

  TFRFME *tfrfme = NULL;

  TFRFME *tfrfme = NULL;

  string binary_name;

  string binary_name;

  if (jss != 1) binary_name = output_path + "/c_TFRFME.hd5";

  if (jss != 1) binary_name = output_path + "/c_TFRFME.hd5";
@@ -330,9 +352,6 @@ void lffft(string data_file, string output_path) { 
	  for (int iy475 = 0; iy475 < nyv; iy475++) {

	  for (int iy475 = 0; iy475 < nyv; iy475++) {

	    for (int ix475 = 0; ix475 < nxv; ix475++) {

	    for (int ix475 = 0; ix475 < nxv; ix475++) {

	      for (int i = 0; i < nlmmt; i++) {

	      for (int i = 0; i < nlmmt; i++) {

		  //double vreal, vimag;

		  //binary_input.read(reinterpret_cast<char *>(&vreal), sizeof(double));

		  //binary_input.read(reinterpret_cast<char *>(&vimag), sizeof(double));

		int row = i;

		int row = i;

		int col = (nyv * nxv * iz475) + (nxv * iy475) + ix475;

		int col = (nyv * nxv * iz475) + (nxv * iy475) + ix475;

		dcomplex value = tfrfme->vec_wsum[nrvc * row + col];

		dcomplex value = tfrfme->vec_wsum[nrvc * row + col];
@@ -393,7 +412,7 @@ void lffft(string data_file, string output_path) { 
			  output_frc, " %18.16lE%18.16lE%18.16lE\n",

			  output_frc, " %18.16lE%18.16lE%18.16lE\n",

			  fffe[2], fffs[2], fffe[2] - fffs[2]

			  fffe[2], fffs[2], fffe[2] - fffs[2]

			  );

			  );

		  } else { // label 450

		} else { // label 450, jss != 1

		  for (int i = 0; i < 3; i++) {

		  for (int i = 0; i < 3; i++) {

		    double value = fffe[i] - fffs[i];

		    double value = fffe[i] - fffs[i];

		    tlfff.write(reinterpret_cast<char *>(&value), sizeof(double));

		    tlfff.write(reinterpret_cast<char *>(&value), sizeof(double));
@@ -430,7 +449,7 @@ void lffft(string data_file, string output_path) { 
			  output_trq, " %18.16lE%18.16lE%18.16lE\n",

			  output_trq, " %18.16lE%18.16lE%18.16lE\n",

			  ffte[2], ffts[2], ffte[2] - ffts[2]

			  ffte[2], ffts[2], ffte[2] - ffts[2]

			  );

			  );

		  } else { // label 470

		} else { // label 470, jss != 1

		  for (int i = 0; i < 3; i++) {

		  for (int i = 0; i < 3; i++) {

		    double value = ffte[i] - ffts[i];

		    double value = ffte[i] - ffts[i];

		    tlfft.write(reinterpret_cast<char *>(&value), sizeof(double));

		    tlfft.write(reinterpret_cast<char *>(&value), sizeof(double));
@@ -462,9 +481,9 @@ void lffft(string data_file, string output_path) { 
  } else {

  } else {

    printf("ERROR: could not open binary input file %s.\n", binary_name.c_str());

    printf("ERROR: could not open binary input file %s.\n", binary_name.c_str());

  }

  }

  } else {

  // } else {

    printf("ERROR: could not open TTMS file.\n");

  //   printf("ERROR: could not open TTMS file.\n");

  }

  // }

  

  

  // Clean up memory

  // Clean up memory

  if (ac != NULL) delete[] ac;

  if (ac != NULL) delete[] ac;
@@ -477,7 +496,8 @@ void lffft(string data_file, string output_path) { 
    delete[] tms;

    delete[] tms;

  }

  }

  if (am0m != NULL) {

  if (am0m != NULL) {

    for (int ai = cil->nlemt - 1; ai > -1; ai--) delete[] am0m[ai];

    //for (int ai = cil->nlemt - 1; ai > -1; ai--) delete[] am0m[ai];

    delete[] vec_am0m;

    delete[] am0m;

    delete[] am0m;

  }

  }

  if (amd != NULL) {

  if (amd != NULL) {