Skip to content
Snippets Groups Projects
Commit 434f882b authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Implement migration of subroutine FRFMER

parent bac213b4
No related branches found
No related tags found
No related merge requests found
#include <fstream>
#include <complex>
#ifndef INCLUDE_TRA_SUBS_H_
#define INCLUDE_TRA_SUBS_H_
#endif
/*! C++ porting of FRFMER
*
* \param nkv: `int` QUESTION: definition?
* \param vkm: `double` QUESTION: definition?
* \param vkv: `double *` QUESTION: definition?
* \param vknmx: `double` QUESTION: definition?
* \param apfafa: `double` QUESTION: definition?
* \param tra: `double` QUESTION: definition?
* \param spd: `double` QUESTION: definition?
* \param rir: `double` QUESTION: definition?
* \param fctn: `double` QUESTION: definition?
* \param le: `int` QUESTION: definition?
* \param lmode: `int` QUESTION: definition?
* \param pmf: `double` QUESTION: definition?
* \param tt1: `fstream &` Handle to first temporary binary file.
* \param tt2: `fstream &` Handle to second temporary binary file.
*/
void frfmer(
int nkv, double vkm, double *vkv, double vknmx, double apfafa, double tra,
double spd, double rir, double ftcn, int le, int lmode, double pmf,
std::fstream &tt1, std::fstream &tt2
) {
const int nlemt = le * (le + 2) * 2;
const std::complex<double> cc0(0.0, 0.0);
std::complex<double> *wk = new std::complex<double>[nlemt]();
double sq = vkm * vkm;
for (int jy90 = 0; jy90 < nkv; jy90++) {
double vky = vkv[jy90];
double sqy = vky * vky;
for (int jx80 = 0; jx80 < nkv; jx80++) {
double vkx = vkv[jx80];
double sqx = vkx * vkx;
double sqn = sqx + sqy;
double vkn = sqrt(sqn);
if (vkn <= vknmx) {
double vkz = sqrt(sq - sqn);
// Would call WAMFF(WK,VKX,VKY,VKZ,LE,APFAFA,TRA,SPD,RIR,FTCN,LMODE,PMF)
for (int j = 0; j < nlemt; j++) {
double vreal = wk[j].real();
double vimag = wk[j].imag();
tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
}
tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
} else { // label 50
for (int j = 0; j < nlemt; j++) {
double vreal = 0.0;
double vimag = 0.0;
tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
}
double vkz = 0.0;
tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
}
} // jx80 loop
} // jy90 loop
delete[] wk;
}
......@@ -9,6 +9,9 @@
//#ifndef INCLUDE_SPH_SUBS_H_
//#include "../include/sph_subs.h"
//#endif
#ifndef INCLUDE_TRA_SUBS_H_
#include "../include/tra_subs.h"
#endif
using namespace std;
......@@ -293,7 +296,7 @@ int main() {
temptape2.open(temp_name2.c_str(), ios::out | ios::binary);
for (int jx = 0; jx < nkv; jx++)
temptape2.write(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));
// Would call FRFMER(NKV,VKM,VKV,VKNMX,APFAFA,TRA,SPD,RIR,FTCN,LM,LMODE,PMF,ITT1,ITT2)
frfmer(nkv, vkm, vkv, vknmx, apfafa, tra, spd, rir, ftcn, lm, lmode, pmf, temptape1, temptape2);
temptape1.write(reinterpret_cast<char *>(&apfafa), sizeof(double)); // Place-holder
temptape1.close();
temptape2.write(reinterpret_cast<char *>(&apfafa), sizeof(double));
......@@ -339,5 +342,6 @@ int main() {
for (int wsi = wsum_size - 1; wsi > -1; wsi--) delete[] wsum[wsi];
delete[] wsum;
}
printf("Done.\n");
return 0;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment