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

Extend migration of LFFFT to C++

parent 1e80f6da

src/include/tra_subs.h

+1 −0
Original line number Diff line number Diff line
@@ -22,6 +22,7 @@ extern void sphar( 
	   double cosrth, double sinrth, double cosrph, double sinrph,

	   int ll, std::complex<double> *ylm

);

extern void thdps(int lm, double ****zpv);

//End of externally defined subroutines



/*! C++ porting of PWMALP

src/trapping/lffft.cpp

+197 −9
Original line number Diff line number Diff line
@@ -21,6 +21,7 @@ using namespace std; 
 *  \param output_path: `string` Directory to write the output files in.

 */

int main() {

  const complex<double> uim(0.0, 1.0);

  string data_file = "../../test_data/trapping/DLFFFT";

  string output_path = ".";

  double ****zpv = NULL, *ac = NULL;
@@ -32,10 +33,12 @@ int main() { 
  complex<double> **am0m = NULL;

  complex<double> **amd = NULL;

  int **indam = NULL;

  complex<double> *tmsm = NULL, *tmse = NULL, *tms = NULL;

  complex<double> *tmsm = NULL, *tmse = NULL, **tms = NULL;

  int jft, jss, jtw;

  int is, le;

  int is, le, nvam = 0;

  double vks, exris;

  cil *p_cil = new cil();

  ccr *p_ccr = new ccr();



  int num_lines = 0;

  string *file_lines = load_file(data_file, &num_lines);
@@ -57,15 +60,172 @@ int main() { 
    ttms.read(reinterpret_cast<char *>(&le), sizeof(int));

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

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

    const int nlem = le * (le + 2);

    const int nlemt = nlem + nlem;

    printf("DEBUG: is = %d\n", is);

    printf("DEBUG: le = %d\n", le);

    printf("DEBUG: vks = %lE\n", vks);

    printf("DEBUG: exris = %lE\n", exris);

    p_cil->le = le;

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

    p_cil->nlemt = p_cil->nlem + p_cil->nlem;

    if (is > 2222) { // label 120

      tms = new complex<double>*[le];

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

      // QUESTION|WARNING: original code uses LM without defining it. Where does it come from?

      int lm = le;

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

	double vreal, vimag;

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

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

	tms[i][0] = complex<double>(vreal, vimag);

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

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

	tms[i][1] = complex<double>(vreal, vimag);

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

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

	tms[i][2] = complex<double>(vreal, vimag);

      } // i loop

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

      tmsm = new complex<double>[le]();

      tmse = new complex<double>[le]();

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

	double vreal, vimag;

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

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

	tmsm[i] = complex<double>(vreal, vimag);

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

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

	tmse[i] = complex<double>(vreal, vimag);

      } // i loop

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

      am0m = new complex<double>*[p_cil->nlemt];

      for (int ai = 0; ai < p_cil->nlemt; ai++) am0m[ai] = new complex<double>[p_cil->nlemt]();

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

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

	  double vreal, vimag;

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

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

	  am0m[i][j] = complex<double>(vreal, vimag);

	} // j loop

      } // i loop

    } else if (is < 0) {

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

      amd = new complex<double>*[nvam];

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

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

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

	  double vreal, vimag;

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

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

	  amd[i][j] = complex<double>(vreal, vimag);

	} // j loop

      } // i loop

      indam = new int*[le];

      int vint;

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

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

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

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

	  indam[i][j] = vint;

	} // j loop

      } // i loop

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

      p_cil->mxmpo = vint;

      p_cil->mxim = vint * 2 - 1;

    }

    // label 150

    ttms.close();

    fstream binary_input;

    string binary_name;

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

    else binary_name = output_path + "c_TWS";

    binary_input.open(binary_name, ios::in | ios::binary);

    if (binary_input.is_open()) {

      int lmode, lm, nkv, nxv, nyv, nzv;

      double vk, exri, an, ff, tra;

      double spd, frsh, exril;

      binary_input.read(reinterpret_cast<char *>(&lmode), sizeof(int));

      binary_input.read(reinterpret_cast<char *>(&lm), sizeof(int));

      binary_input.read(reinterpret_cast<char *>(&nkv), sizeof(int));

      binary_input.read(reinterpret_cast<char *>(&nxv), sizeof(int));

      binary_input.read(reinterpret_cast<char *>(&nyv), sizeof(int));

      binary_input.read(reinterpret_cast<char *>(&nzv), sizeof(int));

      if (lm >= le) {

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

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

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

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

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

	if (vk == vks && exri == exris) {

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

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

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

	  xv = new double[nxv];

	  for (int i = 0; i < nxv; i++) binary_input.read(reinterpret_cast<char *>(&(xv[i])), sizeof(double));

	  yv = new double[nyv];

	  for (int i = 0; i < nyv; i++) binary_input.read(reinterpret_cast<char *>(&(yv[i])), sizeof(double));

	  zv = new double[nzv];

	  for (int i = 0; i < nzv; i++) binary_input.read(reinterpret_cast<char *>(&(zv[i])), sizeof(double));

	  bool goto160 = false;

	  if (jft <= 0) {

	    zpv = new double***[le];

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

	      zpv[zi] = new double**[3];

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

		zpv[zi][zj] = new double*[2];

		for (int zk = 0; zk < 2; zk++) zpv[zi][zj][zk] = new double[2]();

	      } // zj loop

	    } // zi loop

	    thdps(le, zpv);

	    double exdc = exri * exri;

	    double sqk = vk * vk * exdc;

	    p_ccr->cof = 1.0 / sqk;

	    p_ccr->cimu = p_ccr->cof / sqrt(2.0);

	    if (jss != 1) {

	      string tlfff_name = output_path + "/c_TLFFF";

	      fstream tlfff;

	      tlfff.open(tlfff_name.c_str(), ios::out | ios::binary);

	      if (tlfff.is_open()) {

		tlfff.write(reinterpret_cast<char *>(&lmode), sizeof(int));

		tlfff.write(reinterpret_cast<char *>(&le), sizeof(int));

		tlfff.write(reinterpret_cast<char *>(&nkv), sizeof(int));

		tlfff.write(reinterpret_cast<char *>(&nxv), sizeof(int));

		tlfff.write(reinterpret_cast<char *>(&nyv), sizeof(int));

		tlfff.write(reinterpret_cast<char *>(&nzv), sizeof(int));

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

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

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

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

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

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

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

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

		for (int i = 0; i < nxv; i++)

		  tlfff.write(reinterpret_cast<char *>(&(xv[i])), sizeof(double));

		for (int i = 0; i < nyv; i++)

		  tlfff.write(reinterpret_cast<char *>(&(yv[i])), sizeof(double));

		for (int i = 0; i < nzv; i++)

		  tlfff.write(reinterpret_cast<char *>(&(zv[i])), sizeof(double));

		if (jft < 0) goto160 = true;

	      } else { // Should never happen.

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

	      }

	    }

	  }

	  // label 155

	  if (!goto160) {

	    const double sq2i = 1.0 / sqrt(2.0);

	    const complex<double> sq2iti = sq2i * uim;

	    if (jss != 1) {

	      // Would open the ITT file.

	    }

	  }

	  // label 160

	}

      }

      binary_input.close();

    } else {

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

    }

  } else {

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

  }

  

  // Clean up memory

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

  if (ws != NULL) delete[] ws;
@@ -79,7 +239,35 @@ int main() { 
  if (wsl != NULL) delete[] wsl;

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

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

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

  if (tms != NULL) {

    for (int ti = le - 1; ti > -1; ti--) delete[] tms[ti];

    delete[] tms;

  }

  if (am0m != NULL) {

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

    delete[] am0m;

  }

  if (amd != NULL) {

    for (int ai = nvam - 1; ai > -1; ai--) delete[] amd[ai];

    delete[] amd;

  }

  if (indam != NULL) {

    for (int ii = le - 1; ii > -1; ii--) delete[] indam[ii];

    delete[] indam;

  }

  if (zpv != NULL) {

    for (int zi = le - 1; zi > -1; zi--) {

      for (int zj = 2; zj > -1; zj--) {

	for (int zk = 1; zk > -1; zk--) delete[] zpv[zi][zj][zk];

	delete[] zpv[zi][zj];

      } // zj loop

      delete[] zpv[zi];

    } // zi loop

    delete[] zpv;

  }

  delete p_cil;

  delete p_ccr;

  delete[] file_lines;

  printf("Done.\n");

  return 0;

}