commit 71bb30e011e42ec9a643987d363bfb7b6eb51e2c
parent 80f965bc2c1db09fb2ab0bbf3ff6ec257db011f7
Author: Egor Achkasov <eaachkasov@edu.hse.ru>
Date: Tue, 24 Dec 2024 19:16:54 +0100
Impl reflection gains
Diffstat:
6 files changed, 301 insertions(+), 97 deletions(-)
diff --git a/Makefile b/Makefile
@@ -1,2 +1,2 @@
all:
- gcc -g -O0 compute_paths.c test.c -lm -o test
+ gcc -g -O0 -fno-builtin compute_paths.c test.c -lm -o test
diff --git a/compute_paths.c b/compute_paths.c
@@ -7,9 +7,10 @@
#include <math.h> /* for sin, cos, sqrt */
#include <stdint.h>
-#define PI 3.14159265358979323846
-#define SPEED_OF_LIGHT 299792458.0
-#define EPS 1e-6
+#define PI 3.14159265358979323846f /* pi */
+#define SPEED_OF_LIGHT 299792458.0f /* m/s */
+#define EPS 1e-6 /* small number for float comparison */
+#define EPS0 8.854187817e-12 /* vacuum permitivity */
typedef struct {
float x, y, z;
@@ -20,7 +21,9 @@ typedef struct {
} Ray;
typedef struct {
- float a, b, c, d;
+ float eta_re, eta_sqrt_re, eta_inv_re, eta_inv_sqrt_re;
+ float eta_im, eta_sqrt_im, eta_inv_im, eta_inv_sqrt_im;
+ float eta_abs, eta_abs_pow2, eta_abs_inv_sqrt;
} RadioMaterial;
typedef struct {
@@ -33,6 +36,46 @@ typedef struct {
size_t num_indices;
Vec3 *normals;
} Mesh;
+void free_mesh(Mesh *mesh)
+{
+ free(mesh->vertices);
+ free(mesh->rms);
+ free(mesh->indices);
+ free(mesh->rm_indices);
+ free(mesh->normals);
+}
+
+/* ==== COMPLEX OPERATIONS ==== */
+
+void csqrtf(
+ IN float re,
+ IN float im,
+ IN float abs,
+ OUT float *sqrt_re,
+ OUT float *sqrt_im
+)
+{
+ *sqrt_re = sqrtf((re + abs) / 2.f);
+ if (re > -EPS && im < EPS)
+ *sqrt_im = 0;
+ else {
+ *sqrt_im = sqrtf((abs - re) / 2.f);
+ if (im < 0.f) *sqrt_im = -*sqrt_im;
+ }
+}
+void cdivf(
+ IN float a_re,
+ IN float a_im,
+ IN float b_re,
+ IN float b_im,
+ OUT float *c_re,
+ OUT float *c_im
+)
+{
+ float d = b_re * b_re + b_im * b_im;
+ *c_re = (a_re * b_re + a_im * b_im) / d;
+ *c_im = (a_im * b_re - a_re * b_im) / d;
+}
/* ==== VECTOR OPERATIONS ==== */
@@ -85,9 +128,10 @@ Vec3 vec3_scale(const Vec3 *a, float s)
* All the faces must be triangles.
*
* \param mesh_filepath path to the PLY file
+ * \param carrier_frequency the carrier frequency in GHz
* \return the loaded mesh
*/
-Mesh load_mesh_ply(const char *mesh_filepath)
+Mesh load_mesh_ply(const char *mesh_filepath, float carrier_frequency)
{
FILE *f = fopen(mesh_filepath, "rb");
if (!f) {
@@ -167,12 +211,33 @@ Mesh load_mesh_ply(const char *mesh_filepath)
exit(8);
/* RADIO_MATERIALS */
+ /* read a, b, c and d
+ and calculate relative permitivity */
unsigned int name_size;
+ float abcd[4];
+ RadioMaterial *rm;
for (size_t i = 0; i < num_materials; ++i) {
/* skip name */
if (fread(&name_size, sizeof(unsigned int), 1, f) != 1) exit(8);
fseek(f, name_size, SEEK_CUR);
- if (fread(&mesh.rms[i], sizeof(RadioMaterial), 1, f) != 1) exit(8);
+ /* read a, b, c and d */
+ if (fread(abcd, sizeof(float)*4, 1, f) != 1) exit(8);
+ rm = &mesh.rms[i];
+ /* calculate eta */
+ rm->eta_re = abcd[0] * powf(carrier_frequency, abcd[1]);
+ rm->eta_im = (abcd[2] * powf(carrier_frequency, abcd[3]))
+ / (EPS0 * -2.f * PI * carrier_frequency);
+ rm->eta_abs_pow2 = rm->eta_re * rm->eta_re + rm->eta_im * rm->eta_im;
+ rm->eta_abs = sqrtf(rm->eta_abs_pow2);
+ rm->eta_abs_inv_sqrt = 1.f / sqrtf(rm->eta_abs);
+ /* calculate sqrt(eta) */
+ csqrtf(rm->eta_re, rm->eta_im, rm->eta_abs, &rm->eta_sqrt_re, &rm->eta_sqrt_im);
+ /* calculate 1 / eta */
+ rm->eta_inv_re = rm->eta_re / rm->eta_abs_pow2;
+ rm->eta_inv_im = -rm->eta_im / rm->eta_abs_pow2;
+ /* calculate 1 / sqrt(eta) */
+ csqrtf(rm->eta_inv_re, rm->eta_inv_im, 1.f / rm->eta_abs,
+ &rm->eta_inv_sqrt_re, &rm->eta_inv_sqrt_im);
}
/* INDICES */
@@ -208,8 +273,15 @@ Mesh load_mesh_ply(const char *mesh_filepath)
* \param mesh the mesh
* \param t output distance to the hit point. If no hit, t is not modified.
* \param ind output index of the hit triangle. If no hit, i is not modified.
+ * \param theta output angle of incidence
*/
-void moeller_trumbore(Ray *ray, Mesh *mesh, float *t, int32_t *ind)
+void moeller_trumbore(
+ IN Ray *ray,
+ IN Mesh *mesh,
+ OUT float *t,
+ OUT int32_t *ind,
+ OUT float *theta
+)
{
/* TODO BVH */
/* for each triangle */
@@ -240,11 +312,64 @@ void moeller_trumbore(Ray *ray, Mesh *mesh, float *t, int32_t *ind)
if (dist_tmp > EPS && dist_tmp < dist) {
dist = dist_tmp;
*t = dist;
- *ind = i;
+ *ind = i / 3;
+ *theta = acos(vec3_dot(&mesh->normals[i / 3], &ray->d));
+ if (*theta > PI / 2.)
+ *theta = PI - *theta;
}
}
}
+/** Compute reflection R_{eTE} and R_{eTM} coefficients.
+ *
+ * Implements eqs. (31a)-(31b) from ITU-R P.2040-3.
+ *
+ * \param rm the radio material of the hit triangle
+ * \param theta1 the angle of incidence
+ * \param r_te_re output real part of R_{eTE}
+ * \param r_te_im output imaginary part of R_{eTE}
+ * \param r_tm_re output real part of R_{eTM}
+ * \param r_tm_im output imaginary part of R_{eTM}
+ */
+void refl_coefs(
+ IN RadioMaterial *rm,
+ IN float theta1,
+ OUT float *r_te_re,
+ OUT float *r_te_im,
+ OUT float *r_tm_re,
+ OUT float *r_tm_im
+)
+{
+ float sin_theta1 = sinf(theta1);
+ if (rm->eta_abs_inv_sqrt * sin_theta1 > 1.f - EPS) {
+ *r_te_re = *r_tm_re = 1.f;
+ *r_te_im = *r_tm_im = 0.f;
+ return;
+ }
+
+ float sin_theta1_pow2 = sin_theta1 * sin_theta1;
+ float cos_theta2_re = sqrtf(1.f - rm->eta_inv_re / rm->eta_abs_pow2 * sin_theta1_pow2);
+ float cos_theta2_im = sqrtf(1.f + rm->eta_inv_im / rm->eta_abs_pow2 * sin_theta1_pow2);
+
+ /* R_{eTE} */
+ float sqrt_eta_cos_theta2_re = rm->eta_sqrt_re * cos_theta2_re - rm->eta_sqrt_im * cos_theta2_im;
+ float sqrt_eta_cos_theta2_im = rm->eta_sqrt_re * cos_theta2_im + rm->eta_sqrt_im * cos_theta2_re;
+ float cos_theta1 = cosf(theta1);
+ cdivf(cos_theta1 - sqrt_eta_cos_theta2_re, -sqrt_eta_cos_theta2_im,
+ cos_theta1 + sqrt_eta_cos_theta2_re, sqrt_eta_cos_theta2_im,
+ r_te_re, r_te_im);
+
+ /* R_{eTM} */
+ float sqrt_eta_cos_theta1_re = rm->eta_sqrt_re * cos_theta1;
+ float sqrt_eta_cos_theta1_im = rm->eta_sqrt_im * cos_theta1;
+ cdivf(sqrt_eta_cos_theta1_re - cos_theta2_re,
+ sqrt_eta_cos_theta1_im - cos_theta2_im,
+ sqrt_eta_cos_theta1_re + cos_theta2_re,
+ sqrt_eta_cos_theta1_im + cos_theta2_im,
+ r_tm_re, r_tm_im);
+}
+
+
/* ==== MAIN FUNCTION ==== */
void compute_paths(
@@ -258,21 +383,23 @@ void compute_paths(
IN size_t num_tx, /* number of transmitters */
IN size_t num_paths, /* number of paths */
IN size_t num_bounces, /* number of bounces */
- OUT float *a_re, /* output array real parts of gains (num_rx, num_tx, num_paths) */
- OUT float *a_im, /* output array imaginary parts of gains (num_rx, num_tx, num_paths) */
+ OUT float *a_te_re, /* output array real parts of TE gains (num_rx, num_tx, num_paths) */
+ OUT float *a_te_im, /* output array imaginary parts of TE gains (num_rx, num_tx, num_paths) */
+ OUT float *a_tm_re, /* output array real parts of TM gains (num_rx, num_tx, num_paths) */
+ OUT float *a_tm_im, /* output array imaginary parts of TM gains (num_rx, num_tx, num_paths) */
OUT float *tau /* output array of delays (num_rx, num_tx, num_paths) */
)
{
/* Load the scene */
- Mesh mesh = load_mesh_ply(mesh_filepath);
+ Mesh mesh = load_mesh_ply(mesh_filepath, carrier_frequency);
/* Calculate fibonacci sphere */
Vec3 *ray_directions = (Vec3*)malloc(num_paths * sizeof(Vec3));
float k, phi, theta;
for (size_t i = 0; i < num_paths; ++i) {
- k = (float)i + .5;
- phi = acos(1. - 2. * k / num_paths);
- theta = PI * (1. + sqrt(5.)) * k;
+ k = (float)i + .5f;
+ phi = acos(1.f - 2.f * k / num_paths);
+ theta = PI * (1.f + sqrtf(5.f)) * k;
ray_directions[i] = (Vec3){
cos(theta) * sin(phi),
sin(theta) * sin(phi),
@@ -293,50 +420,91 @@ void compute_paths(
}
free(ray_directions);
- /* Calculate the paths:
- - sum distances for each path into tau
- - a???
+ /* Bounce the rays. On each bounce:
+ * - Add path length / SPEED_OF_LIGHT to tau
+ * - Update gains using eqs. (31a)-(31b) from ITU-R P.2040-3
*/
- /* TODO calculate a and tau in moeller_trumbore */
+ /* TODO add active mask */
/* shape (num_tx, num_paths) */
- Vec3 *hits = (Vec3*)malloc(num_tx * num_paths * sizeof(Vec3));
- int32_t *hit_indices = (int32_t*)malloc(num_tx * num_paths * sizeof(int32_t));
- tau = (float*)memset(tau, 0, num_tx * num_paths * sizeof(float));
- float t;
+ float *tau_t = (float*)calloc(num_tx * num_paths, sizeof(float));
+ float *a_te_re_t = (float*)malloc(num_tx * num_paths * sizeof(float));
+ float *a_te_im_t = (float*)calloc(num_tx * num_paths, sizeof(float));
+ float *a_tm_re_t = (float*)malloc(num_tx * num_paths * sizeof(float));
+ float *a_tm_im_t = (float*)calloc(num_tx * num_paths, sizeof(float));
+ for (size_t i = 0; i < num_tx * num_paths; ++i)
+ a_te_re_t[i] = a_tm_re_t[i] = 1.f;
+ float t, r_te_re, r_te_im, r_tm_re, r_tm_im;
+ float a_te_re_new, a_te_im_new, a_tm_re_new, a_tm_im_new;
int32_t ind;
Ray *r;
- Vec3 *h;
- size_t off;
+ Vec3 *h = (Vec3*)malloc(sizeof(Vec3)); /* temp vector */
for (size_t i = 0; i < num_bounces; ++i) {
- off = 0;
- for (size_t j = 0; j < num_tx; ++j)
- for (size_t k = 0; k < num_paths; ++k) {
- t = -1.;
- ind = -1;
- r = &rays[off];
- moeller_trumbore(r, &mesh, &t, &ind);
- tau[off] += t;
- hit_indices[off] = ind;
- h = &hits[off];
- *h = vec3_scale(&r->d, t);
- *h = vec3_add(h, &r->o);
- r->o = *h;
- off += 1;
- }
+ for (size_t off = 0; off != num_tx * num_paths; ++off) {
+ /* Init */
+ t = -1.f;
+ ind = -1;
+ r = &rays[off];
+ /* Find the hit point and trinagle.
+ Calculate an angle of incidence */
+ moeller_trumbore(r, &mesh, &t, &ind, &theta);
+ /* Calcilate the reflection coefficients
+ R_{eTE} and R_{eTM} */
+ refl_coefs(&mesh.rms[mesh.rm_indices[ind]],
+ theta,
+ &r_te_re, &r_te_im,
+ &r_tm_re, &r_tm_im);
+ /* Update the gains */
+ a_te_re_new = a_te_re_t[off] * r_te_re - a_te_im_t[off] * r_te_im;
+ a_te_im_new = a_te_re_t[off] * r_te_im + a_te_im_t[off] * r_te_re;
+ a_tm_re_new = a_tm_re_t[off] * r_tm_re - a_tm_im_t[off] * r_tm_im;
+ a_tm_im_new = a_tm_re_t[off] * r_tm_im + a_tm_im_t[off] * r_tm_re;
+ a_te_re_t[off] = a_te_re_new;
+ a_te_im_t[off] = a_te_im_new;
+ a_tm_re_t[off] = a_tm_re_new;
+ a_tm_im_t[off] = a_tm_im_new;
+ /* Update the delay */
+ tau[off] += t / SPEED_OF_LIGHT;
+ /* Advance the ray to the hit point */
+ *h = vec3_scale(&r->d, t);
+ *h = vec3_add(h, &r->o);
+ r->o = *h;
+ }
}
+ free(h);
- /* Calculate tau */
- for (size_t i = 0; i < num_paths; ++i)
- tau[i] /= SPEED_OF_LIGHT;
+ /* Do the final propagations from the final hit points
+ * directly to the receivers.
+ * Also broadcast the gains and delays to rx
+ * by adding num_rx dimension to the local tau_t and a_t.
+ */
+ size_t off_a;
+ for (size_t i = 0; i < num_rx; ++i)
+ for (size_t off = 0; off != num_tx * num_paths; ++off) {
+ /* Init */
+ t = -1.f;
+ ind = -1;
+ r = &rays[off];
+ off_a = i * num_tx * num_paths + off;
+ /* See if the return ray hits any abstacle */
+ moeller_trumbore(r, &mesh, &t, &ind, &theta);
+ if (ind != -1) { /* Return ray is blocked */
+ tau[off_a] = -1.;
+ a_te_im[off_a] = a_te_re[off_a] = 0.f;
+ a_tm_im[off_a] = a_tm_re[off_a] = 0.f;
+ }
+ /* Set the gains */
+ a_te_re[off_a] = a_te_re_t[off];
+ a_te_im[off_a] = a_te_im_t[off];
+ a_tm_re[off_a] = a_tm_re_t[off];
+ a_tm_im[off_a] = a_tm_im_t[off];
+ /* Set the delays */
+ r->d = (Vec3){rx_positions[i * 3], rx_positions[i * 3 + 1], rx_positions[i * 3 + 2]};
+ r->d = vec3_sub(&r->o, &r->d);
+ t = sqrtf(vec3_dot(&r->d, &r->d));
+ tau[off_a] = tau_t[off] + t / SPEED_OF_LIGHT;
+ }
/* Free */
- free(hits);
- free(hit_indices);
free(rays);
- /* Free the mesh */
- free(mesh.vertices);
- free(mesh.rms);
- free(mesh.indices);
- free(mesh.rm_indices);
- free(mesh.normals);
+ free_mesh(&mesh);
}
diff --git a/compute_paths.h b/compute_paths.h
@@ -18,12 +18,14 @@
* \param tx_positions transmitter positions, shape (num_tx, 3)
* \param rx_velocities receiver velocities, shape (num_rx, 3)
* \param tx_velocities transmitter velocities, shape (num_tx, 3)
- * \param carrier_frequency carrier frequency in Hz. Must be > 0.0
+ * \param carrier_frequency carrier frequency in GHz. Must be > 0.0
* \param num_rx number of receivers. Must be > 0
* \param num_tx number of transmitters. Must be > 0
* \param num_paths number of paths to compute. Must be > 0
- * \param a_re output array of real parts of gains, shape (num_rx, num_tx, num_paths)
- * \param a_im output array of imaginary parts of gains, shape (num_rx, num_tx, num_paths)
+ * \param a_te_re output array of real parts of transverse electric gains, shape (num_rx, num_tx, num_paths)
+ * \param a_te_im output array of imaginary parts of transverse electric gains, shape (num_rx, num_tx, num_paths)
+ * \param a_tm_re output array of real parts of transverse magnetic gains, shape (num_rx, num_tx, num_paths)
+ * \param a_tm_im output array of imaginary parts of transverse magnetic gains, shape (num_rx, num_tx, num_paths)
* \param tau output array of delays in seconds, shape (num_rx, num_tx, num_paths)
*/
void compute_paths(
@@ -37,8 +39,10 @@ void compute_paths(
IN size_t num_tx, /* number of transmitters */
IN size_t num_paths, /* number of paths */
IN size_t num_bounces, /* number of bounces */
- OUT float *a_re, /* output array real parts of gains (num_rx, num_tx, num_paths) */
- OUT float *a_im, /* output array imaginary parts of gains (num_rx, num_tx, num_paths) */
+ OUT float *a_te_re, /* output array real parts of TE gains (num_rx, num_tx, num_paths) */
+ OUT float *a_te_im, /* output array imaginary parts of TE gains (num_rx, num_tx, num_paths) */
+ OUT float *a_tm_re, /* output array real parts of TM gains (num_rx, num_tx, num_paths) */
+ OUT float *a_tm_im, /* output array imaginary parts of TM gains (num_rx, num_tx, num_paths) */
OUT float *tau /* output array of delays (num_rx, num_tx, num_paths) */
);
diff --git a/compute_paths_pybind11.cpp b/compute_paths_pybind11.cpp
@@ -11,7 +11,10 @@ extern "C" {
namespace py = pybind11;
-std::tuple< py::array_t<float>, py::array_t<float>, py::array_t<float> >
+std::tuple<
+ py::array_t<float>, py::array_t<float>,
+ py::array_t<float>, py::array_t<float>,
+ py::array_t<float> >
compute_paths_wrapper(
const std::string &mesh_filepath,
py::array_t<float> rx_positions,
@@ -31,8 +34,10 @@ compute_paths_wrapper(
py::buffer_info tx_vel_info = tx_velocities.request();
// Output arrays
- float *a_im = new float[num_tx * num_paths];
- float *a_re = new float[num_tx * num_paths];
+ float *a_te_im = new float[num_tx * num_paths];
+ float *a_te_re = new float[num_tx * num_paths];
+ float *a_tm_im = new float[num_tx * num_paths];
+ float *a_tm_re = new float[num_tx * num_paths];
float *tau = new float[num_tx * num_paths];
// Call the C function
@@ -42,30 +47,36 @@ compute_paths_wrapper(
(const float*)tx_pos_info.ptr, // Tx positions
(const float*)rx_vel_info.ptr, // Rx velocities
(const float*)tx_vel_info.ptr, // Tx velocities
- carrier_frequency,
+ carrier_frequency / 1e9, // Carrier frequency in GHz
(size_t)num_rx,
(size_t)num_tx,
(size_t)num_paths,
(size_t)num_bounces,
- a_re,
- a_im,
- tau
+ a_te_re, a_te_im, a_tm_re, a_tm_im, // Gains
+ tau // Delays
);
// Convert output arrays into numpy arrays for easy use in Python
// TODO prevent copying data
// TODO construct np.complex64 array for a
- py::array_t<float> a_im_array = py::array_t<float>({num_tx, num_paths}, a_im);
- py::array_t<float> a_re_array = py::array_t<float>({num_tx, num_paths}, a_re);
- py::array_t<float> tau_array = py::array_t<float>({num_tx, num_paths}, tau);
+ py::array_t<float> a_te_im_array = py::array_t<float>({num_rx, num_tx, num_paths}, a_te_im);
+ py::array_t<float> a_te_re_array = py::array_t<float>({num_rx, num_tx, num_paths}, a_te_re);
+ py::array_t<float> a_tm_im_array = py::array_t<float>({num_rx, num_tx, num_paths}, a_tm_im);
+ py::array_t<float> a_tm_re_array = py::array_t<float>({num_rx, num_tx, num_paths}, a_tm_re);
+ py::array_t<float> tau_array = py::array_t<float>({num_rx, num_tx, num_paths}, tau);
// Deallocate arrays
- delete[] a_im;
- delete[] a_re;
+ delete[] a_te_im;
+ delete[] a_te_re;
+ delete[] a_tm_im;
+ delete[] a_tm_re;
delete[] tau;
// Return the results as a tuple (gains, delays)
- return std::make_tuple(a_re_array, a_im_array, tau_array);
+ return std::make_tuple(
+ a_te_re_array, a_te_im_array,
+ a_tm_re_array, a_tm_im_array,
+ tau_array);
}
PYBIND11_MODULE(rt, m) {
diff --git a/test.c b/test.c
@@ -5,30 +5,50 @@
int main(int argc, char **argv)
{
- if (argc < 2) {
- fprintf(stderr, "Usage: %s <path_to_mesh.ply>\n", argv[0]);
- return 1;
- }
+ if (argc < 2) {
+ fprintf(stderr, "Usage: %s <path_to_mesh.ply>\n", argv[0]);
+ return 1;
+ }
- float rx_positions[3] = {0.0, 0.0, 2.5};
- float tx_positions[3] = {0.0, 0.0, 2.5};
- float rx_velocities[3] = {0.0, 0.0, 0.0};
- float tx_velocities[3] = {0.0, 0.0, 0.0};
- float carrier_frequency = 2.4e9; /* 2.4 GHz */
- size_t num_rx = 1;
- size_t num_tx = 1;
- size_t num_paths = 10000;
- size_t num_bounces = 3;
- float *a_re = (float*)malloc(num_tx * num_paths * sizeof(float));
- float *a_im = (float*)malloc(num_tx * num_paths * sizeof(float));
- float *tau = (float*)malloc(num_tx * num_paths * sizeof(float));
- compute_paths(
- argv[1],
- rx_positions,
- tx_positions,
- rx_velocities,
- tx_velocities,
- carrier_frequency,
- num_rx, num_tx, num_paths, num_bounces,
- a_re, a_im, tau);
+ float rx_positions[3] = {0.0, 0.0, 2.5};
+ float tx_positions[3] = {0.0, 0.0, 2.5};
+ float rx_velocities[3] = {0.0, 0.0, 0.0};
+ float tx_velocities[3] = {0.0, 0.0, 0.0};
+ float carrier_frequency = 2.4e9; /* 2.4 GHz */
+ size_t num_rx = 1;
+ size_t num_tx = 1;
+ size_t num_paths = 10000;
+ size_t num_bounces = 3;
+ float *a_te_re = (float*)malloc(num_rx * num_tx * num_paths * sizeof(float));
+ float *a_te_im = (float*)malloc(num_rx * num_tx * num_paths * sizeof(float));
+ float *a_tm_re = (float*)malloc(num_rx * num_tx * num_paths * sizeof(float));
+ float *a_tm_im = (float*)malloc(num_rx * num_tx * num_paths * sizeof(float));
+ float *tau = (float*)malloc(num_rx * num_tx * num_paths * sizeof(float));
+ compute_paths(
+ argv[1],
+ rx_positions,
+ tx_positions,
+ rx_velocities,
+ tx_velocities,
+ carrier_frequency,
+ num_rx, num_tx, num_paths, num_bounces,
+ a_te_re, a_te_im, a_tm_re, a_tm_im,
+ tau);
+
+ /* Save the results */
+ FILE *f = fopen("a_te_re.bin", "wb");
+ fwrite(a_te_re, sizeof(float), num_rx * num_tx * num_paths, f);
+ fclose(f);
+ f = fopen("a_te_im.bin", "wb");
+ fwrite(a_te_im, sizeof(float), num_rx * num_tx * num_paths, f);
+ fclose(f);
+ f = fopen("a_tm_re.bin", "wb");
+ fwrite(a_tm_re, sizeof(float), num_rx * num_tx * num_paths, f);
+ fclose(f);
+ f = fopen("a_tm_im.bin", "wb");
+ fwrite(a_tm_im, sizeof(float), num_rx * num_tx * num_paths, f);
+ fclose(f);
+ f = fopen("tau.bin", "wb");
+ fwrite(tau, sizeof(float), num_rx * num_tx * num_paths, f);
+ fclose(f);
}
diff --git a/test.py b/test.py
@@ -14,7 +14,7 @@ num_paths = 10000
num_bounces = 3
# Call compute_paths
-a_im, a_re, tau = compute_paths(
+a_te_re, a_te_im, a_tm_re, a_tm_im, tau = compute_paths(
mesh_filepath,
rx_positions,
tx_positions,
@@ -27,5 +27,6 @@ a_im, a_re, tau = compute_paths(
num_bounces
)
-print("Gains:", a_re + 1.j*a_im)
-print("Delays:", tau)
+print(f"TE Gains: {a_te_re + 1.j*a_te_im}")
+print(f"TM Gains: {a_tm_re + 1.j*a_tm_im}")
+print(f"Delays: {tau}")
