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

Implement step-by-step conversion of FRFME

parent fb7f935e

src/libnptm/tra_subs.cpp

+14 −4
Original line number Diff line number Diff line
@@ -297,7 +297,7 @@ void pwmalp(complex<double> **w, double *up, double *un, complex<double> *ylm, i 
    int lf = l20 + 1;

    int lftl = lf * l20;

    double x = 1.0 * lftl;

    complex<double> cl = std::pow((four_pi / sqrt(x)) * uim, 1.0 * l20);

    complex<double> cl = (four_pi / sqrt(x)) * std::pow(uim, 1.0 * l20);

    int mv = l20 + lf;

    int m = -lf;

    for (int mf20 = 1; mf20 <= mv; mf20++) {
@@ -387,7 +387,7 @@ void wamff( 
  if (!(onz && onx && ony)) {

    rho = sqrt(x * x + y * y);

    cph = x / rho;

    sph = y = rho;

    sph = y / rho;

  } else {

    if (onz) { // label 10

      cph = 1.0;
@@ -429,6 +429,7 @@ void wamff( 
    }

  }

  if (lmode == 0 || sth != 0.0) { // label 25

    bool skip62 = false;

    ylm[nlmp - 1] = cc0;

    sphar(cth, sth, cph, sph, lm, ylm);

    up[0] = cph * cth;
@@ -461,16 +462,19 @@ void wamff( 
	cfam *= (sth * pmf * ftc1);

	for (int i52 = 0; i52 < nlmmt; i52++) wk[i52] = w[i52][0] * cfam;

	// returns

	skip62 = true;

      } else if (lmode == 3) { // label 53

	ftc2 = 2.0 * cthl / (cthl  + cth * rir);

	cfam *= (sth * pmf * ftc2);

	for (int i55 = 0; i55 < nlmmt; i55++) wk[i55] = w[i55][1] * cfam;

	// returns

	skip62 = true;

      }

      if (lmode == 0 || lmode == 1) { //label 48

	cf1 = cph * ftc1 * cfam;

	cf2 = -sph * ftc2 * cfam;

	// goes to 62

	skip62 = false;

      }

    } else { // label 57

      double fam = tra * std::exp(-apfa * apfa) / sqrt(cth);
@@ -479,25 +483,31 @@ void wamff( 
	double f2 = -sph * fam;

	for (int i58 = 0; i58 < nlmmt; i58++) wk[i58] = w[i58][0] * f1 + w[i58][1] * f2;

	// returns

	skip62 = true;

      } else if (lmode == 1) { // label 60

	cfam = (pmf * sth * fam) * (cph * uim * sph);

	cf1 = cph * cfam;

	cf2 = -sph * cfam;

	// follows on to 62

	skip62 = false;

      } else if (lmode == 2) { // label 65

	fam *= (pmf * sth);

	for (int i67 = 0; i67 < nlmmt; i67++) wk[i67] = w[i67][0] * fam;

	// returns

	skip62 = true;

      } else if (lmode == 3) { // label 68

	fam *= (pmf * sth);

	for (int i70 = 0; i70 < nlmmt; i70++) wk[i70] = w[i70][1] * fam;

	// returns

	skip62 = true;

      }

    }

    if (!skip62) {

      if (lmode == 0 || lmode == 1) { // label 62

	for (int i63 = 0; i63 < nlmmt; i63++) wk[i63] = w[i63][0] * cf1 + w[i63][1] * cf2;

      }

    }

  }

  // Clean up memory

  delete[] up;

  delete[] un;

src/trapping/cfrfme.cpp

+48 −37
Original line number Diff line number Diff line
@@ -49,14 +49,14 @@ void frfme(string data_file, string output_path) { 
  str_target = m.suffix().str();

  regex_search(str_target, m, re);

  int jlml = stoi(m.str());

  int lmode, lm, nks, nkv;

  double vk, exri, an, ff, tra;

  double exdc, wp, xip, xi;

  int idfc, nxi;

  double apfafa, pmf, spd, rir, ftcn, fshmx;

  double vxyzmx, delxyz, vknmx, delk, delks;

  double frsh, exril;

  int nlmmt, nrvc;

  int lmode = 0, lm = 0, nks = 0, nkv = 0;

  double vk = 0.0, exri = 0.0, an = 0.0, ff = 0.0, tra = 0.0;

  double exdc = 0.0, wp = 0.0, xip = 0.0, xi = 0.0;

  int idfc = 0, nxi = 0;

  double apfafa = 0.0, pmf = 0.0, spd = 0.0, rir = 0.0, ftcn = 0.0, fshmx = 0.0;

  double vxyzmx = 0.0, delxyz = 0.0, vknmx = 0.0, delk = 0.0, delks = 0.0;

  double frsh = 0.0, exril = 0.0;

  int nlmmt = 0, nrvc = 0;

  // Vector size variables

  int vkzm_size, wsum_size;

  // End of vector size variables
@@ -78,6 +78,7 @@ void frfme(string data_file, string output_path) { 
      tfrfme.read(reinterpret_cast<char *>(&spd),  sizeof(double));

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

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

      vkv = new double[nkv]();

      xv = new double[nxv]();

      yv = new double[nyv]();

      zv = new double[nzv]();
@@ -85,10 +86,9 @@ void frfme(string data_file, string output_path) { 
      for (int yi = 0; yi < nxv; yi++) tfrfme.read(reinterpret_cast<char *>(&(yv[yi])), sizeof(double));

      for (int zi = 0; zi < nxv; zi++) tfrfme.read(reinterpret_cast<char *>(&(zv[zi])), sizeof(double));

      fstream temptape2;

      string tempname2 = output_path + "c_TEMPTAPE2";

      string tempname2 = output_path + "/c_TEMPTAPE2";

      temptape2.open(tempname2.c_str(), ios::in | ios::binary);

      if (temptape2.is_open()) {

	//vkv = new double[nkv]();

	for (int jx = 0; jx < nkv; jx++) temptape2.read(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));

	vkzm = new double*[nkv];

	vkzm_size = nkv;
@@ -228,6 +228,11 @@ void frfme(string data_file, string output_path) { 
	nlmmt = lm * (lm + 2) * 2;

	nks = nksh * 2;

	nkv = nks + 1;

	// Array initialization

	vkv = new double[nkv]();

	vkzm = new double*[nkv];

	for (int vi = 0; vi < nkv; vi++) vkzm[vi] = new double[nkv];

	// End of array initialization

	double vkm = vk * exri;

	vknmx = vk * an;

	delk = vknmx / nksh;
@@ -239,7 +244,6 @@ void frfme(string data_file, string output_path) { 
	int nxv = nxs + 1;

	int nxshpo = nxsh + 1;

	xv = new double[nxv]();

	//xv[nxsh] = 0.0;

	for (int i24 = nxshpo; i24 <= nxs; i24++) {

	  xv[i24] = xv[i24 - 1] + delxyz;

	  xv[nxv - i24 - 1] = -xv[i24];
@@ -248,7 +252,6 @@ void frfme(string data_file, string output_path) { 
	int nyv = nys + 1;

	int nyshpo = nysh + 1;

	yv = new double[nyv]();

	//yv[nysh] = 0.0;

	for (int i25 = nyshpo; i25 <= nys; i25++) {

	  yv[i25] = yv[i25 - 1] + delxyz;

	  yv[nyv - i25 - 1] = -yv[i25];
@@ -257,18 +260,12 @@ void frfme(string data_file, string output_path) { 
	int nzv = nzs + 1;

	int nzshpo = nzsh + 1;

	zv = new double[nzv]();

	//zv[nysh] = 0.0;

	for (int i27 = nzshpo; i27 <= nzs; i27++) {

	  zv[i27] = zv[i27 - 1] + delxyz;

	  zv[nzv - i27 - 1] = -zv[i27];

	} // i27 loop

	int nrvc = nxv * nyv * nzv;

	int nkshpo = nksh + 1;

	wsum = new complex<double>*[nlmmt];

	wsum_size = nlmmt;

	for (int wsi = 0; wsi < nlmmt; wsi++) wsum[wsi] = new complex<double>[nrvc]();

	vkv = new double[nkv]();

	// vkv[nksh] = 0.0;

	for (int i28 = nkshpo; i28 <= nks; i28++) {

	  vkv[i28] = vkv[i28 - 1] + delk;

	  vkv[nkv - i28 - 1] = -vkv[i28];
@@ -319,21 +316,34 @@ void frfme(string data_file, string output_path) { 
	  temptape2.write(reinterpret_cast<char *>(&nrvc), sizeof(int));

	  temptape2.close();

	  temptape2.open("c_TEMPTAPE2", ios::in | ios::binary);

	  vkv = new double[nkv]();

	  for (int jx = 0; jx < nkv; jx++)

	    temptape2.read(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));

	  vkzm = new double*[nkv];

	  vkzm_size = nkv;

	  for (int vki = 0; vki < nkv; vki++) vkzm[vki] = new double[nkv]();

	  for (int jx = 0; jx < nkv; jx++) {

	    double value = 0.0;

	    temptape2.read(reinterpret_cast<char *>(&value), sizeof(double));

	    vkv[jx] = value;

	  }

	  for (int jy40 = 0; jy40 < nkv; jy40++) {

	    for (int jx40 = 0; jx40 < nkv; jx40++)

	      temptape2.read(reinterpret_cast<char *>(&(vkzm[jx40][jy40])), sizeof(double));

	    for (int jx40 = 0; jx40 < nkv; jx40++) {

	      double value = 0.0;

	      temptape2.read(reinterpret_cast<char *>(&value), sizeof(double));

	      vkzm[jx40][jy40] = value;

	    }

	  } // jy40 loop

	  temptape2.close();

	  wk = new complex<double>[nlmmt];

	  if (wsum != NULL) {

	    for (int wsi = wsum_size - 1; wsi > -1; wsi--) delete[] wsum[wsi];

	    delete[] wsum;

	  }

	  wsum = new complex<double>*[nlmmt];

	  for (int j80 = jlmf - 1; j80 < jlml; j80++) {

	    // w matrix

	    if (w != NULL) {

	      for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];

	      delete[] w;

	    }

	    w = new complex<double>*[nkv];

	    for (int wi = 0; wi < nkv; wi++) w[wi] = new complex<double>[nkv]();

	  for (int j80 = jlmf - 1; j80 < jlml; j80++) {

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

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

	    temptape1.open(temp_name1.c_str(), ios::in | ios::binary);

	    for (int jy50 = 0; jy50 < nkv; jy50++) {

	      for (int jx50 = 0; jx50 < nkv; jx50++) {
@@ -348,6 +358,7 @@ void frfme(string data_file, string output_path) { 
	    } // jy50 loop

	    temptape1.close();

	    int ixyz = 0;

	    wsum[j80] = new complex<double>[nrvc]();

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

	      double z = zv[iz75]  + frsh;

	      for (int iy70 = 0; iy70 < nyv; iy70++) {
@@ -364,11 +375,11 @@ void frfme(string data_file, string output_path) { 
		    double vkzl = vkzm[nkv - 1][jy60];

		    complex<double> phasl = exp(uim * (vkx * x + vky * y + vkzl * z));

		    complex<double> sumx = 0.5 * (w[0][jy60] * phasf + w[nkv - 1][jy60] * phasl);

		    for (int jx55 = 1; jx55 < nks; jx55++) {

		      vkx = vkv[jx55];

		      double vkz = vkzm[jx55][jy60];

		    for (int jx55 = 2; jx55 <= nks; jx55++) {

		      vkx = vkv[jx55 - 1];

		      double vkz = vkzm[jx55 - 1][jy60];

		      complex<double> phas = exp(uim * (vkx * x + vky * y + vkz * z));

		      sumx += (w[jx55][jy60] * phas);

		      sumx += (w[jx55 - 1][jy60] * phas);

		    } // jx55 loop

		    if (jy60 == 0 || jy60 == nkv - 1) sumx *= 0.5;

		    sumy += sumx;
@@ -439,14 +450,14 @@ void frfme(string data_file, string output_path) { 
    for (int vki = vkzm_size - 1; vki > -1; vki--) delete[] vkzm[vki];

    delete[] vkzm;

  }

  if (w != NULL) {

    for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];

    delete[] w;

  }

  if (wsum != NULL) {

    for (int wsi = wsum_size - 1; wsi > -1; wsi--) delete[] wsum[wsi];

    delete[] wsum;

  }

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

  if (w != NULL) {

    for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];

    delete[] w;

  }

  printf("Done.\n");

}