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

Throw exception errors if configuration mismatches are detected in lffft().

parent 5a4eb213

src/trapping/clffft.cpp

+11 −135
Original line number Original line Diff line number Diff line
@@ -273,23 +273,6 @@ void lffft(string data_file, string output_path) { 
	  ccr->cof = 1.0 / sqk;

	  ccr->cof = 1.0 / sqk;

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

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

	  if (jss != 1) {

	  if (jss != 1) {

	    // string tlfff_name = output_path + "/c_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));

	    toi->set_param("lmode", 1.0 * lmode);

	    toi->set_param("lmode", 1.0 * lmode);

	    toi->set_param("le", 1.0 * le);

	    toi->set_param("le", 1.0 * le);

	    toi->set_param("nkv", 1.0 * nkv);

	    toi->set_param("nkv", 1.0 * nkv);
@@ -302,86 +285,24 @@ void lffft(string data_file, string output_path) { 
	    toi->set_param("frsh", frsh);

	    toi->set_param("frsh", frsh);

	    toi->set_param("exril", exril);

	    toi->set_param("exril", exril);

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

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

	      // double x = _xv[i];

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

	      toi->vec_x[i] = _xv[i];

	      toi->vec_x[i] = _xv[i];

	    }

	    }

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

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

	      // double y = _yv[i];

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

	      toi->vec_y[i] = _yv[i];

	      toi->vec_y[i] = _yv[i];

	    }

	    }

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

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

	      // double z = _zv[i];

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

	      toi->vec_z[i] = _zv[i];

	      toi->vec_z[i] = _zv[i];

	    }

	    }

	      if (jft < 0) goto160 = true;

	      if (jft < 0) goto160 = true;

	  }

	  }

	}

	}

	// label 155

	// label 155

	// if (!goto160) {

	//   if (jss != 1) {

	//     string tlfft_name = output_path + "/c_TLFFT";

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

	//     if (tlfft.is_open()) {

	//       string outgrid_name = output_path + "/c_grid_scale.txt";

	//       FILE *outgrid = fopen(outgrid_name.c_str(), "w");

	//       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++) {

	// 	double x = _xv[i];

	// 	if (nxv - i == i + 1) {

	// 	  double ximo = _xv[i - 1];

	// 	  double xipo = _xv[i + 1];

	// 	  if (ximo < 0.0 && xipo > 0.0) {

	// 	    double limo = log10(-ximo);

	// 	    double lipo = log10(xipo);

	// 	    double logi = (x > 0.0) ? log10(x) : log10(-x);

	// 	    if (logi < (limo + lipo) / 2.0 - 10.0) x = 0.0;

	// 	  }

	// 	}

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

	// 	fprintf(outgrid, "  %24.16lE\n", x);

	//       }

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

	// 	double y = _yv[i];

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

	//       }

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

	// 	double z = _zv[i];

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

	//       }

	//       fclose(outgrid);

	//     } else { // Should never happen.

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

	//     }

	//   }

	// }

	// label 160

	// label 160

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

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

	const int nlmmt = nlmm + nlmm;

	const int nlmmt = nlmm + nlmm;

	const int nrvc = nxv * nyv * nzv;

	const int nrvc = nxv * nyv * nzv;

	ws = new dcomplex[nlmmt]();

	ws = new dcomplex[nlmmt]();

	if (lm > le) wsl = new dcomplex[nlmmt]();

	if (lm > le) wsl = new dcomplex[nlmmt]();

	// FORTRAN writes two output formatted files without opening them

	// explicitly. It is assumed thay can be opened here.

	// string force_cs_name = output_path + "/c_force_cs.txt";

	// string torque_cs_name = output_path + "/c_torque_cs.txt";

	// FILE *output_frc = fopen(force_cs_name.c_str(), "w");

	// FILE *output_trq = fopen(torque_cs_name.c_str(), "w");

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

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

	  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++) {
@@ -444,34 +365,10 @@ void lffft(string data_file, string output_path) { 
		  toi->vec_csf1[3 * out_index + 2] = fffe[2];

		  toi->vec_csf1[3 * out_index + 2] = fffe[2];

		  toi->vec_csf2[3 * out_index + 2] = fffs[2];

		  toi->vec_csf2[3 * out_index + 2] = fffs[2];

		  toi->vec_csf3[3 * out_index + 2] = fffe[2] - fffs[2];

		  toi->vec_csf3[3 * out_index + 2] = fffe[2] - fffs[2];

		  // fprintf(

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

		  // 	  fffe[0], fffs[0], fffe[0] - fffs[0]

		  // 	  );

		  // fprintf(

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

		  // 	  fffe[1], fffs[1], fffe[1] - fffs[1]

		  // 	  );

		  // fprintf(

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

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

		  // 	  );

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

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

		  toi->vec_csf1[out_index] = fffe[0] - fffs[0];

		  toi->vec_csf1[out_index] = fffe[0] - fffs[0];

		  toi->vec_csf2[out_index] = fffe[1] - fffs[1];

		  toi->vec_csf2[out_index] = fffe[1] - fffs[1];

		  toi->vec_csf3[out_index] = fffe[2] - fffs[2];

		  toi->vec_csf3[out_index] = fffe[2] - fffs[2];

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

		  //   fprintf(

		  // 	    output_frc, " %5d%4d%4d%15.4lE%15.4lE%15.4lE\n",

		  // 	    ix475 + 1, iy475 + 1, iz475 + 1,

		  // 	    fffe[0] - fffs[0], fffe[1] - fffs[1], fffe[2] - fffs[2]

		  // 	    );

		  // }

		}

		}

		if (jft < 0) goto475 = true;

		if (jft < 0) goto475 = true;

		delete[] fffe;

		delete[] fffe;
@@ -483,6 +380,7 @@ void lffft(string data_file, string output_path) { 
		double *ffts = new double[3]();

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

		ffrt(ac, ws, ffte, ffts, cil);

		ffrt(ac, ws, ffte, ffts, cil);

		if (jss == 1) {

		if (jss == 1) {

		  // Writes to 67

		  toi->vec_cst1[3 * out_index] = ffte[0];

		  toi->vec_cst1[3 * out_index] = ffte[0];

		  toi->vec_cst2[3 * out_index] = ffts[0];

		  toi->vec_cst2[3 * out_index] = ffts[0];

		  toi->vec_cst3[3 * out_index] = ffte[0] - ffts[0];

		  toi->vec_cst3[3 * out_index] = ffte[0] - ffts[0];
@@ -492,35 +390,10 @@ void lffft(string data_file, string output_path) { 
		  toi->vec_cst1[3 * out_index + 2] = ffte[2];

		  toi->vec_cst1[3 * out_index + 2] = ffte[2];

		  toi->vec_cst2[3 * out_index + 2] = ffts[2];

		  toi->vec_cst2[3 * out_index + 2] = ffts[2];

		  toi->vec_cst3[3 * out_index + 2] = ffte[2] - ffts[2];

		  toi->vec_cst3[3 * out_index + 2] = ffte[2] - ffts[2];

		  // Writes to 67

		  // fprintf(

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

		  // 	  ffte[0], ffts[0], ffte[0] - ffts[0]

		  // 	  );

		  // fprintf(

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

		  // 	  ffte[1], ffts[1], ffte[1] - ffts[1]

		  // 	  );

		  // fprintf(

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

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

		  // 	  );

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

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

		  toi->vec_cst1[out_index] = ffte[0] - ffts[0];

		  toi->vec_cst1[out_index] = ffte[0] - ffts[0];

		  toi->vec_cst2[out_index] = ffte[1] - ffts[1];

		  toi->vec_cst2[out_index] = ffte[1] - ffts[1];

		  toi->vec_cst3[out_index] = ffte[2] - ffts[2];

		  toi->vec_cst3[out_index] = ffte[2] - ffts[2];

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

		  //   fprintf(

		  // 	    output_trq, " %5d%4d%4d%15.4lE%15.4lE%15.4lE\n",

		  // 	    ix475 + 1, iy475 + 1, iz475 + 1,

		  // 	    ffte[0] - ffts[0], ffte[1] - ffts[1], ffte[2] - ffts[2]

		  // 	    );

		  // }

		}

		}

		delete[] ffte;

		delete[] ffte;

		delete[] ffts;

		delete[] ffts;
@@ -528,13 +401,16 @@ void lffft(string data_file, string output_path) { 
	    } // ix475 loop

	    } // ix475 loop

	  } // iy475 loop

	  } // iy475 loop

	} // iz475 loop

	} // iz475 loop

	// if (jss != 1) {

      } else { // vk != vks || exri != exris

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

	message = "ERROR: T-matrix mismatch!\n";

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

	if (vk != vks) message = "ERROR: calculation wavelength not matching T-matrix wavelength!\n";

	// }

	else if (exri != exris) message = "ERROR: external medium not matching T-matrix external medium!\n";

	// fclose(output_frc);

	throw(UnrecognizedConfigurationException(message));

	// fclose(output_trq);

      }

      }

    } else { // lm < le: ERROR!

      message = "ERROR: calculation order " + to_string(lm) + " smaller than T-matrix order "

	+ to_string(le) + "!\n";

      throw(UnrecognizedConfigurationException(message));

    }

    }

    toi->write(output_path + "/c_LFFFT.hd5", "HDF5");

    toi->write(output_path + "/c_LFFFT.hd5", "HDF5");

    if (jtw == 1) toi->write(output_path + "/c_", "ASCII");

    if (jtw == 1) toi->write(output_path + "/c_", "ASCII");