Commit 1a12652c authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Convert SAMPOA subroutine to C++

parent 6ffd3797

src/include/tra_subs.h

+39 −2
Original line number Diff line number Diff line
@@ -7,11 +7,11 @@ 
#endif



//Structures for TRAPPING (that will probably move to Commons)

struct cil {

struct CIL {

  int le, nlem, nlemt, mxmpo, mxim;

};



struct ccr {

struct CCR {

  double cof, cimu;

};

//End of TRAPPING structures
@@ -286,3 +286,40 @@ void frfmer( 
  } // jy90 loop

  delete[] wk;

}



/*! C++ porting of SAMPOA

 *

 * \param ac: Vector of complex. QUESTION: definition?

 * \param tms: Matrix of complex. QUESTION: definition?

 * \param ws: Vector of complex. QUESTION: definition?

 * \param cil: `CIL *` Pointer to a CIL structure.

 *

 * This function builds the AC vector using TMS and WS.

 */

void sampoa(

	    std::complex<double> *ac, std::complex<double> **tms, std::complex<double> *ws,

	    CIL *cil

) {

  std::complex<double> **tm = new std::complex<double>*[2];

  tm[0] = new std::complex<double>[2]();

  tm[1] = new std::complex<double>[2]();

  int i = 0;

  ac = new std::complex<double>[cil->nlemt]();

  for (int l20 = 0; l20 < cil->le; l20++) {

    tm[0][0] = tms[l20][0];

    tm[0][1] = tms[l20][1];

    tm[1][1] = tms[l20][2];

    tm[1][0] = tm[0][1];

    int l = l20 + 1;

    int ltpo = l + l + 1;

    for (int im20 = 0; im20 < ltpo; im20++) {

      int ie = i + cil->nlem;

      ac[i] = tm[0][0] * ws[i] + tm[0][1] * ws[ie];

      ac[ie] = tm[1][0] * ws[i] + tm[1][1] * ws[ie];

      i++;

    } // im20 loop

  } // l20 loop

  delete[] tm[1];

  delete[] tm[0];

  delete[] tm;

}

src/trapping/lffft.cpp

+132 −28
Original line number Diff line number Diff line
@@ -22,14 +22,15 @@ using namespace std; 
 */

int main() {

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

  const double sq2i = 1.0 / sqrt(2.0);

  const complex<double> sq2iti = sq2i * uim;



  fstream tlfff, tlfft;

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

  string output_path = ".";

  double ****zpv = NULL, *ac = NULL;

  complex<double> *ws = NULL;

  double *fffe = NULL, *fffs = NULL;

  double *ffte = NULL, *ffts = NULL;

  double ****zpv = NULL;

  double *xv = NULL, *yv = NULL, *zv = NULL;

  complex<double> *wsl = NULL;

  complex<double> *ac = NULL, *ws = NULL, *wsl = NULL;

  complex<double> **am0m = NULL;

  complex<double> **amd = NULL;

  int **indam = NULL;
@@ -37,8 +38,8 @@ int main() { 
  int jft, jss, jtw;

  int is, le, nvam = 0;

  double vks, exris;

  cil *p_cil = new cil();

  ccr *p_ccr = new ccr();

  CIL *cil = new CIL();

  CCR *ccr = new CCR();



  int num_lines = 0;

  string *file_lines = load_file(data_file, &num_lines);
@@ -60,9 +61,9 @@ 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));

    p_cil->le = le;

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

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

    cil->le = le;

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

    cil->nlemt = cil->nlem + 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]();
@@ -93,10 +94,10 @@ int main() { 
	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++) {

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

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

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

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

	  double vreal, vimag;

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

	  ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
@@ -125,8 +126,8 @@ int main() { 
	} // j loop

      } // i loop

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

      p_cil->mxmpo = vint;

      p_cil->mxim = vint * 2 - 1;

      cil->mxmpo = vint;

      cil->mxim = vint * 2 - 1;

    }

    // label 150

    ttms.close();
@@ -174,11 +175,10 @@ int main() { 
	    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);

	    ccr->cof = 1.0 / sqk;

	    ccr->cimu = 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));
@@ -209,13 +209,121 @@ int main() { 
	  }

	  // 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.

	      string tlfft_name = output_path + "/c_TLFFT";

	      tlfft.open(tlfft_name.c_str(), ios::out | ios::binary);

	      if (tlfft.is_open()) {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	      } else { // Should never happen.

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

	      }

	    }

	  }

	  // label 160

	  const int nlmm = lm * (lm + 2);

	  const int nlmmt = nlmm + nlmm;

	  ws = new complex<double>[nlmmt]();

	  if (lm > le) wsl = new complex<double>[nlmmt]();

	  for (int iz475 = 0; iz475 < nzv; iz475++) {

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

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

		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));

		  if (lm <= le) {

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

		  } else { // label 170

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

		    for (int i175 = 0; i175 < cil->nlem; i175++) {

		      int ie = i175 + cil->nlem;

		      int iel = i175 + nlmm;

		      ws[i175] = wsl[i175];

		      ws[ie] = wsl[iel];

		    } // i175 loop

		  }

		} // i loop

		// label 180

		bool goto475 = false;

		if (is != 2222) {

		  if (is != 1111) {

		    if (is > 0) { // Goes to 305

		      // Would call CAMP(AC,AM0M,WS)

		    } else if (is < 0) { // Goes to 405

		      // Would call CZAMP(AC,AMD,INDAM,WS)

		    }

		  } else {

		    // Would call SAMP(AC,TMSM,TMSE,WS)

		  }

		} else {

		  sampoa(ac, tms, ws, cil);

		}

		// label 445

		if (jft <= 0) {

		  double *fffe = new double[3]();

		  double *fffs = new double[3]();

		  // Would call FFRF(ZPV,AC,WS,FFFE,FFFS)

		  if (jss == 1) {

		    // Writes to 66

		  } else { // label 450

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

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

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

		    }

		    if (jtw == 1) {

		      // Writes to 66

		    }

		  }

		  if (jft < 0) goto475 = true;

		  delete[] fffe;

		  delete[] fffs;

		}

		// label 460

		if (!goto475) {

		  double *ffte = new double[3]();

		  double *ffts = new double[3]();

		  // Would call FFRT(AC,WS,FFTE,FFTS)

		  if (jss == 1) {

		    // Writes to 67

		  } else { // label 470

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

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

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

		    }

		    if (jtw == 1) {

		      // Writes to 67

		    }

		  }

		  delete[] ffte;

		  delete[] ffts;

		}

	      } // ix475 loop

	    } // iy475 loop

	  } // iz475 loop

	  if (jss != 1) {

	    if (jft <= 0) tlfff.close();

	    if (jft >= 0) tlfft.close();

	  }

	}

      }

      binary_input.close();
@@ -229,10 +337,6 @@ int main() { 
  // Clean up memory

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

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

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

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

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

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

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

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

  if (zv != NULL) delete[] zv;
@@ -244,7 +348,7 @@ int main() { 
    delete[] tms;

  }

  if (am0m != NULL) {

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

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

    delete[] am0m;

  }

  if (amd != NULL) {
@@ -265,8 +369,8 @@ int main() { 
    } // zi loop

    delete[] zpv;

  }

  delete p_cil;

  delete p_ccr;

  delete cil;

  delete ccr;

  delete[] file_lines;

  printf("Done.\n");

  return 0;