compute_paths.c (28970B)
1 #include "../inc/compute_paths.h" /* for compute_paths */ 2 #include "../inc/scene.h" /* for Scene, Mesh, Material */ 3 #include "../inc/materials.h" /* for g_materials */ 4 #include "../inc/common.h" /* for IN, OUT, PERROR_CLEANUP_EXIT */ 5 #include "../inc/vec3.h" /* for Vec3, vec3_sub, vec3_cross, vec3_dot, vec3_normalize */ 6 #include "../inc/ray.h" /* for Ray */ 7 8 #include <stddef.h> /* for size_t */ 9 #include <stdlib.h> /* for exit, malloc, free */ 10 #include <string.h> /* for strlen, sprintf */ 11 #include <stdio.h> /* for fopen, FILE, fclose */ 12 #include <math.h> /* for sin, cos, sqrt */ 13 #include <stdint.h> /* for int32_t, uint8_t */ 14 #include <float.h> /* for __FLT_EPSILON__ */ 15 16 /* ==== CONSTANTS ==== */ 17 18 #define PI 3.14159265358979323846f /* pi */ 19 #define SPEED_OF_LIGHT 299792458.0f /* m/s */ 20 21 /* (2w, w) binomial coefficients for w = 0, 1, ..., 19 */ 22 const float BINOMIAL_2W_W[20] = { 23 1.f, 2.f, 6.f, 20.f, 70.f, 24 252.f, 924.f, 3432.f, 12870.f, 48620.f, 25 184756.f, 705432.f, 2704156.f, 10400600.f, 40116600.f, 26 155117520.f, 601080390.f, 2333606220.f, 9075135300.f, 35345263800.f 27 }; 28 /* (aplha, k) binomial coefficients for alpha = 0, 1, ..., 19 with k = 0, 1, ..., alpha */ 29 const float* BINOMIAL_ALPHA_K[20] = { 30 (float[]){ 31 1.f 32 }, 33 (float[]){ 34 1.f, 1.f 35 }, 36 (float[]){ 37 1.f, 2.f, 1.f 38 }, 39 (float[]){ 40 1.f, 3.f, 3.f, 1.f 41 }, 42 (float[]){ 43 1.f, 4.f, 6.f, 4.f, 1.f 44 }, 45 (float[]){ 46 1.f, 5.f, 10.f, 10.f, 5.f, 47 1.f 48 }, 49 (float[]){ 50 1.f, 6.f, 15.f, 20.f, 15.f, 51 6.f, 1.f 52 }, 53 (float[]){ 54 1.f, 7.f, 21.f, 35.f, 35.f, 55 21.f, 7.f, 1.f 56 }, 57 (float[]){ 58 1.f, 8.f, 28.f, 56.f, 70.f, 59 56.f, 28.f, 8.f, 1.f 60 }, 61 (float[]){ 62 1.f, 9.f, 36.f, 84.f, 126.f, 63 126.f, 84.f, 36.f, 9.f, 1.f 64 }, 65 (float[]){ 66 1.f, 10.f, 45.f, 120.f, 210.f, 67 252.f, 210.f, 120.f, 45.f, 10.f, 68 1.f 69 }, 70 (float[]){ 71 1.f, 11.f, 55.f, 165.f, 330.f, 72 462.f, 462.f, 330.f, 165.f, 55.f, 73 11.f, 1.f 74 }, 75 (float[]){ 76 1.f, 12.f, 66.f, 220.f, 495.f, 77 792.f, 924.f, 792.f, 495.f, 220.f, 78 66.f, 12.f, 1.f 79 }, 80 (float[]){ 81 1.f, 13.f, 78.f, 286.f, 82 715.f, 1287.f, 1716.f, 1716.f, 1287.f, 83 715.f, 286.f, 78.f, 13.f, 1.f 84 }, 85 (float[]){ 86 1.f, 14.f, 91.f, 364.f, 1001.f, 87 2002.f, 3003.f, 3432.f, 3003.f, 2002.f, 88 1001.f, 364.f, 91.f, 14.f, 1.f 89 }, 90 (float[]){ 91 1.f, 15.f, 105.f, 455.f, 1365.f, 92 3003.f, 5005.f, 6435.f, 6435.f, 5005.f, 93 3003.f, 1365.f, 455.f, 105.f, 15.f, 94 1.f 95 }, 96 (float[]){ 97 1.f, 16.f, 120.f, 560.f, 1820.f, 98 4368.f, 8008.f, 11440.f, 12870.f, 11440.f, 99 8008.f, 4368.f, 1820.f, 560.f, 120.f, 100 16.f, 1.f 101 }, 102 (float[]){ 103 1.f, 17.f, 136.f, 680.f, 2380.f, 104 6188.f, 12376.f, 19448.f, 24310.f, 24310.f, 105 19448.f, 12376.f, 6188.f, 2380.f, 680.f, 106 136.f, 17.f, 1.f 107 }, 108 (float[]){ 109 1.f, 18.f, 153.f, 816.f, 3060.f, 110 8568.f, 18564.f, 31824.f, 43758.f, 48620.f, 111 43758.f, 31824.f, 18564.f, 8568.f, 3060.f, 112 816.f, 153.f, 18.f, 1.f 113 }, 114 (float[]){ 115 1.f, 19.f, 171.f, 969.f, 3876.f, 116 11628.f, 27132.f, 50388.f, 75582.f, 92378.f, 117 92378.f, 75582.f, 50388.f, 27132.f, 11628.f, 118 3876.f, 969.f, 171.f, 19.f, 1.f 119 } 120 }; 121 122 /* ==== STRUCTS ==== */ 123 124 /* Radio material eta and additional parameters */ 125 typedef struct { 126 /* eta */ 127 float eta_re, eta_sqrt_re, eta_inv_re, eta_inv_sqrt_re; 128 float eta_im, eta_sqrt_im, eta_inv_im, eta_inv_sqrt_im; 129 float eta_abs, eta_abs_pow2, eta_abs_inv_sqrt; 130 /* r = 1.0 - s */ 131 float r; 132 } MaterialPrecomputed; 133 134 /* ==== COMPLEX OPERATIONS ==== */ 135 136 void csqrtf( 137 IN float re, 138 IN float im, 139 IN float abs, 140 OUT float *sqrt_re, 141 OUT float *sqrt_im 142 ) 143 { 144 *sqrt_re = sqrtf((re + abs) / 2.f); 145 if (fabsf(im) < __FLT_EPSILON__ && re >= -__FLT_EPSILON__) 146 *sqrt_im = 0; 147 else { 148 *sqrt_im = sqrtf((abs - re) / 2.f); 149 if (im < 0.f) *sqrt_im = -*sqrt_im; 150 } 151 } 152 void cdivf( 153 IN float a_re, 154 IN float a_im, 155 IN float b_re, 156 IN float b_im, 157 OUT float *c_re, 158 OUT float *c_im 159 ) 160 { 161 float d = b_re * b_re + b_im * b_im; 162 *c_re = (a_re * b_re + a_im * b_im) / d; 163 *c_im = (a_im * b_re - a_re * b_im) / d; 164 } 165 166 /* ==== PRECOMPUTED GLOBALS ==== */ 167 168 /* Materials etas, calucalted for the given carrier_frequency */ 169 MaterialPrecomputed g_materials_precomputed[NUM_G_MATERIALS]; 170 171 void precompute_materials( 172 IN Scene* scene, 173 IN float carrier_frequency 174 ) 175 { 176 uint8_t material_eta_isdone[NUM_G_MATERIALS] = {0}; 177 for (uint32_t i = 0; i != scene->num_meshes; ++i) { 178 uint32_t mat_ind = scene->meshes[i].material_index; 179 if (material_eta_isdone[mat_ind]) 180 continue; 181 Material *mat = &g_materials[mat_ind]; 182 MaterialPrecomputed *mat_precomp = &g_materials_precomputed[mat_ind]; 183 /* calculate eta */ 184 mat_precomp->eta_re = mat->a * powf(carrier_frequency, mat->b); 185 /* eq. 12 */ 186 mat_precomp->eta_im = (mat->c * powf(carrier_frequency, mat->d)) 187 / (0.0556325027352135f * carrier_frequency); 188 mat_precomp->eta_abs_pow2 = mat_precomp->eta_re * mat_precomp->eta_re 189 + mat_precomp->eta_im * mat_precomp->eta_im; 190 mat_precomp->eta_abs = sqrtf(mat_precomp->eta_abs_pow2); 191 mat_precomp->eta_abs_inv_sqrt = 1.f / sqrtf(mat_precomp->eta_abs); 192 /* calculate sqrt(eta) */ 193 csqrtf(mat_precomp->eta_re, mat_precomp->eta_im, 194 mat_precomp->eta_abs, &mat_precomp->eta_sqrt_re, 195 &mat_precomp->eta_sqrt_im); 196 /* calculate 1 / eta */ 197 mat_precomp->eta_inv_re = mat_precomp->eta_re / mat_precomp->eta_abs_pow2; 198 mat_precomp->eta_inv_im = -mat_precomp->eta_im / mat_precomp->eta_abs_pow2; 199 /* calculate 1 / sqrt(eta) */ 200 csqrtf(mat_precomp->eta_inv_re, mat_precomp->eta_inv_im, 201 1.f / mat_precomp->eta_abs, 202 &mat_precomp->eta_inv_sqrt_re, &mat_precomp->eta_inv_sqrt_im); 203 /* calculate r */ 204 mat_precomp->r = 1.f - mat->s; 205 } 206 } 207 208 void precompute_normals(Scene* scene) { 209 /* Calculate normals */ 210 for (uint32_t i = 0; i != scene->num_meshes; ++i) { 211 Mesh *mesh = &scene->meshes[i]; 212 mesh->ns = (Vec3*)malloc(mesh->num_triangles * sizeof(Vec3)); 213 for (uint32_t j = 0; j != mesh->num_triangles; ++j) { 214 Vec3* v1 = (Vec3*)&mesh->vs[mesh->is[j * 3]]; 215 Vec3* v2 = (Vec3*)&mesh->vs[mesh->is[j * 3 + 1]]; 216 Vec3* v3 = (Vec3*)&mesh->vs[mesh->is[j * 3 + 2]]; 217 Vec3 e1 = vec3_sub(v2, v1); 218 Vec3 e2 = vec3_sub(v3, v1); 219 Vec3 n = vec3_cross(&e1, &e2); 220 n = vec3_normalize(&n); 221 memcpy(&mesh->ns[j], &n, sizeof(Vec3)); 222 } 223 } 224 } 225 226 /* ==== RT ==== */ 227 228 /** Compute Moeleer-Trumbore intersection algorithm. 229 * 230 * \param ray ray to cast to the mesh 231 * \param scene the scene to cast the ray to 232 * \param t output distance to the hit point. If no hit, t is not modified. 233 * \param mesh_ind output index of the hit mesh. If no hit, i is not modified. 234 * \param face_ind output index of the hit triangle. If no hit, i is not modified. 235 * \param theta output angle of incidence 236 */ 237 void moeller_trumbore( 238 IN Ray *ray, 239 IN Scene *scene, 240 OUT float *t, 241 OUT uint32_t *mesh_ind, 242 OUT uint32_t *face_ind, 243 OUT float *theta 244 ) 245 { 246 /* TODO BVH */ 247 /* for each triangle */ 248 Vec3 *v1, *v2, *v3, *n; 249 Vec3 e1, e2, re2_cross, s, se1_cross; 250 float d, u, v; 251 float dist = 1e9; 252 float dist_tmp; 253 for (uint32_t i = 0; i != scene->num_meshes; ++i) { 254 Mesh *mesh = &scene->meshes[i]; 255 for (uint32_t j = 0; j != mesh->num_triangles; ++j) { 256 v1 = (Vec3*)&mesh->vs[mesh->is[3 * j]]; 257 v2 = (Vec3*)&mesh->vs[mesh->is[3 * j + 1]]; 258 v3 = (Vec3*)&mesh->vs[mesh->is[3 * j + 2]]; 259 e1 = vec3_sub(v2, v1); 260 e2 = vec3_sub(v3, v1); 261 re2_cross = vec3_cross(&ray->d, &e2); 262 d = vec3_dot(&e1, &re2_cross); 263 if (d > -__FLT_EPSILON__ && d < __FLT_EPSILON__) continue; 264 s = vec3_sub(&ray->o, v1); 265 u = vec3_dot(&s, &re2_cross) / d; 266 if ((u < 0. && fabs(u) > __FLT_EPSILON__) 267 || (u > 1. && fabs(u - 1.) > __FLT_EPSILON__)) 268 continue; 269 se1_cross = vec3_cross(&s, &e1); 270 v = vec3_dot(&ray->d, &se1_cross) / d; 271 if ((v < 0. && fabs(v) > __FLT_EPSILON__) 272 || (u + v > 1. && fabs(u + v - 1.) > __FLT_EPSILON__)) 273 continue; 274 dist_tmp = vec3_dot(&e2, &se1_cross) / d; 275 if (dist_tmp > __FLT_EPSILON__ && dist_tmp < dist) { 276 dist = dist_tmp; 277 *t = dist; 278 *mesh_ind = i; 279 *face_ind = j; 280 n = (Vec3*)&mesh->ns[j]; 281 *theta = acos(vec3_dot(n, &ray->d)); 282 if (*theta > PI / 2.) 283 *theta = PI - *theta; 284 } 285 } 286 } 287 } 288 289 /** Compute reflection R_{eTE} and R_{eTM} coefficients. 290 * 291 * Implements eqs. (31a)-(31b) from ITU-R P.2040-3. 292 * 293 * \param material_index the index of the material 294 * \param theta1 the angle of incidence 295 * \param r_te_re output real part of R_{eTE} 296 * \param r_te_im output imaginary part of R_{eTE} 297 * \param r_tm_re output real part of R_{eTM} 298 * \param r_tm_im output imaginary part of R_{eTM} 299 */ 300 void refl_coefs( 301 IN uint32_t material_index, 302 IN float theta1, 303 OUT float *r_te_re, 304 OUT float *r_te_im, 305 OUT float *r_tm_re, 306 OUT float *r_tm_im 307 ) 308 { 309 MaterialPrecomputed *rm = &g_materials_precomputed[material_index]; 310 float sin_theta1 = sinf(theta1); 311 if (rm->eta_abs_inv_sqrt * sin_theta1 > 1.f - __FLT_EPSILON__) { 312 *r_te_re = *r_tm_re = 1.f; 313 *r_te_im = *r_tm_im = 0.f; 314 return; 315 } 316 317 /* eq. 33 */ 318 float sin_theta1_pow2 = sin_theta1 * sin_theta1; 319 float cos_theta2_re = sqrtf(1.f + rm->eta_inv_re / rm->eta_abs_pow2 * sin_theta1_pow2); 320 float cos_theta2_im = sqrtf(1.f - rm->eta_inv_im / rm->eta_abs_pow2 * sin_theta1_pow2); 321 322 /* R_{eTE}, eq. 31a */ 323 float sqrt_eta_cos_theta2_re = rm->eta_sqrt_re * cos_theta2_re - rm->eta_sqrt_im * cos_theta2_im; 324 float sqrt_eta_cos_theta2_im = rm->eta_sqrt_re * cos_theta2_im + rm->eta_sqrt_im * cos_theta2_re; 325 float cos_theta1 = cosf(theta1); 326 cdivf(cos_theta1 - sqrt_eta_cos_theta2_re, -sqrt_eta_cos_theta2_im, 327 cos_theta1 + sqrt_eta_cos_theta2_re, sqrt_eta_cos_theta2_im, 328 r_te_re, r_te_im); 329 330 /* R_{eTM}, eq. 31b */ 331 float sqrt_eta_cos_theta1_re = rm->eta_sqrt_re * cos_theta1; 332 float sqrt_eta_cos_theta1_im = rm->eta_sqrt_im * cos_theta1; 333 cdivf(sqrt_eta_cos_theta1_re - cos_theta2_re, 334 sqrt_eta_cos_theta1_im - cos_theta2_im, 335 sqrt_eta_cos_theta1_re + cos_theta2_re, 336 sqrt_eta_cos_theta1_im + cos_theta2_im, 337 r_tm_re, r_tm_im); 338 339 /* Apply reflection reduction factor */ 340 *r_te_re *= rm->r; 341 *r_te_im *= rm->r; 342 *r_tm_re *= rm->r; 343 *r_tm_im *= rm->r; 344 } 345 346 /** Calculate scattering coefficients. 347 * 348 * Implements a directive scattering model inspired by Blaunstein et al. (DOI 10.1109/TAP.2006.888422). 349 * Models scattering from rough surfaces (e.g., vegetation) with polarization dependence. 350 * 351 * \param theta_s scattering angle (radians, between scattered direction and surface normal) 352 * \param theta_i angle of incidence (radians, between incident ray and surface normal) 353 * \param material_index index into g_materials array 354 * \param a_te_re output real part of TE scattering coefficient 355 * \param a_te_im output imaginary part of TE scattering coefficient 356 * \param a_tm_re output real part of TM scattering coefficient 357 * \param a_tm_im output imaginary part of TM scattering coefficient 358 */ 359 void scat_coefs( 360 IN float theta_s, 361 IN float theta_i, 362 IN uint32_t material_index, 363 OUT float *a_te_re, 364 OUT float *a_te_im, 365 OUT float *a_tm_re, 366 OUT float *a_tm_im 367 ) 368 { 369 Material *mat = &g_materials[material_index]; 370 371 /* Precompute trigonometric values */ 372 float cos_theta_s = cosf(theta_s); 373 float cos_theta_i = cosf(theta_i); 374 float sin_theta_i = sinf(theta_i); 375 376 /* Scattering directivity pattern: exponential decay with angle difference */ 377 /* f = s * exp(-alpha * |theta_s - theta_i|) models directive scattering */ 378 float theta_diff = fabsf(theta_s - theta_i); 379 float f = mat->s * expf(-mat->s1_alpha * theta_diff); 380 381 /* Roughness factor: reduces specular component as alpha increases */ 382 float roughness = 1.0f / (1.0f + mat->s1_alpha); // Simplified, 0 to 1 383 float specular = roughness * cos_theta_s; // Specular-like term 384 float diffuse = (1.0f - roughness) * cos_theta_s; // Diffuse term 385 386 /* TE and TM coefficients (real parts) */ 387 /* TE: perpendicular to plane of incidence, less affected by angle */ 388 float te_real = f * (specular + diffuse); 389 /* TM: parallel to plane of incidence, stronger angular dependence */ 390 float tm_real = f * (specular * cos_theta_i + diffuse); 391 392 /* Phase shift (imaginary parts) due to surface roughness or path difference */ 393 /* Simplified: assume small phase shift proportional to roughness and angle */ 394 float phase_shift = mat->s1_alpha * sin_theta_i * 0.1f; // Arbitrary scaling 395 float te_imag = te_real * sinf(phase_shift); 396 float tm_imag = tm_real * sinf(phase_shift); 397 398 /* Normalize to ensure energy conservation (optional, adjust as needed) */ 399 float norm = sqrtf(te_real * te_real + te_imag * te_imag + 400 tm_real * tm_real + tm_imag * tm_imag); 401 if (norm > 1e-6f) { // Avoid division by zero 402 te_real /= norm; 403 te_imag /= norm; 404 tm_real /= norm; 405 tm_imag /= norm; 406 } 407 408 /* Output coefficients */ 409 *a_te_re = te_real; 410 *a_te_im = te_imag; 411 *a_tm_re = tm_real; 412 *a_tm_im = tm_imag; 413 414 /* TODO: Incorporate s2 and s3 for additional roughness or frequency effects */ 415 } 416 417 /* ==== MAIN FUNCTION ==== */ 418 419 void compute_paths( 420 IN Scene *scene, /* Pointer to a loaded scene */ 421 IN Vec3 *rx_pos, /* shape (num_rx, 3) */ 422 IN Vec3 *tx_pos, /* shape (num_tx, 3) */ 423 IN Vec3 *rx_vel, /* shape (num_rx, 3) */ 424 IN Vec3 *tx_vel, /* shape (num_tx, 3) */ 425 IN float carrier_frequency_GHz, /* > 0.0 (IN GHz!) */ 426 IN size_t num_rx, /* number of receivers */ 427 IN size_t num_tx, /* number of transmitters */ 428 IN size_t num_paths, /* number of paths */ 429 IN size_t num_bounces, /* number of bounces */ 430 OUT ChannelInfo *chanInfo_los, /* output LoS channel information */ 431 OUT RaysInfo *raysInfo_los, /* output LoS rays information */ 432 OUT ChannelInfo *chanInfo_scat, /* output scatter channel information */ 433 OUT RaysInfo *raysInfo_scat /* output scatter rays information */ 434 ) 435 { 436 /* Precompute globals */ 437 precompute_materials(scene, carrier_frequency_GHz); 438 precompute_normals(scene); 439 440 /* Calculate a fibonacci sphere to init rays that leave the transmitters */ 441 /* rays shape = (num_tx, num_paths) */ 442 Ray *rays = (Ray*)malloc(num_tx * num_paths * sizeof(Ray)); 443 for (size_t path = 0; path < num_paths; ++path) { 444 float k = (float)path + .5f; 445 float phi = acos(1.f - 2.f * k / num_paths); 446 float theta = PI * (1.f + sqrtf(5.f)) * k; 447 Vec3 d = (Vec3){ 448 cos(theta) * sin(phi), 449 sin(theta) * sin(phi), 450 cos(phi) 451 }; 452 for (size_t tx = 0; tx < num_tx; ++tx) { 453 rays[tx * num_paths + path].o = tx_pos[tx]; 454 rays[tx * num_paths + path].d = d; 455 } 456 } 457 458 /* Init bouncing rays gains and delays */ 459 /* shape (num_tx, num_paths) */ 460 float *a_te_re = (float*)malloc(num_tx * num_paths * sizeof(float)); 461 float *a_te_im = (float*)calloc(num_tx * num_paths, sizeof(float)); 462 float *a_tm_re = (float*)malloc(num_tx * num_paths * sizeof(float)); 463 float *a_tm_im = (float*)calloc(num_tx * num_paths, sizeof(float)); 464 for (size_t i = 0; i < num_tx * num_paths; ++i) 465 a_te_re[i] = a_tm_re[i] = 1.f; 466 float *tau = (float*)calloc(num_tx * num_paths, sizeof(float)); 467 468 /* Active rays bitmask. 1 if the ray is still traced, 0 if it has left the scene */ 469 uint8_t *active = (uint8_t*)malloc((num_tx * num_paths / 8 + 1) * sizeof(uint8_t)); 470 for (size_t i = 0; i < num_tx * num_paths / 8 + 1; ++i) 471 raysInfo_scat->rays_active[i] = active[i] = 0xff; 472 size_t num_active = num_tx * num_paths; 473 474 /* Temporary variables */ 475 float t, r_te_re, r_te_im, r_tm_re, r_tm_im; 476 float theta; /* incidence angle */ 477 float theta_s; /* scattering angle */ 478 uint32_t mesh_ind, face_ind, material_ind; 479 Ray* r; 480 481 /* Friis free space loss multiplier. 482 Multiply this by a distance and take the second power to get a free space loss. */ 483 float carrier_frequency_hz = carrier_frequency_GHz * 1e9; 484 float free_space_loss_multiplier = 4.f * PI * carrier_frequency_hz / SPEED_OF_LIGHT; 485 float free_space_loss; 486 487 /* Doppler frequency shift carrier frequency multiplier */ 488 float doppler_shift_multiplier = carrier_frequency_hz / SPEED_OF_LIGHT; 489 490 /* Initialize the doppler frequency shift using the tx velocities */ 491 /* The scatter->freq_shift has the shape of (num_rx, num_tx, num_bounces*num_paths) */ 492 /* The first num_tx*num_paths elements are the same for each rx and bounce */ 493 /* so we can initialize them only once and memcpy them to the rest */ 494 for (size_t tx = 0; tx != num_tx; ++tx) 495 for (size_t path = 0; path != num_paths; ++path) { 496 size_t off_rays = tx * num_paths + path; 497 size_t off_scat = tx * num_paths * num_bounces + path; 498 chanInfo_scat->freq_shift[off_scat] = vec3_dot(&tx_vel[tx], &rays[off_rays].d); 499 chanInfo_scat->freq_shift[off_scat] *= doppler_shift_multiplier; 500 } 501 for (size_t bounce = 1; bounce < num_bounces; ++bounce) 502 if (!memcpy(chanInfo_scat->freq_shift + num_tx * num_paths * bounce, 503 chanInfo_scat->freq_shift, num_tx * num_paths * sizeof(float))) 504 exit(70); 505 for (size_t rx = 1; rx < num_rx; ++rx) 506 if (!memcpy(chanInfo_scat->freq_shift + num_tx * num_paths * num_bounces * rx, 507 chanInfo_scat->freq_shift, num_tx * num_paths * num_bounces * sizeof(float))) 508 exit(70); 509 510 /***************************************************************************/ 511 /* LoS paths */ 512 /***************************************************************************/ 513 514 /* Consider imaginary parts of the LoS gains to be 0 in any way */ 515 for (size_t off = 0; off != num_rx * num_tx; ++off) 516 chanInfo_los->a_te_im[off] = chanInfo_los->a_tm_im[off] = 0.f; 517 518 /* Calculate LoS paths directly from tx to rx */ 519 /* off = rx * num_tx + tx */ 520 for (size_t rx = 0, off = 0; rx != num_rx; ++rx) 521 for (size_t tx = 0; tx != num_tx; ++tx, ++off) { 522 t = -1.f; 523 mesh_ind = face_ind = -1; 524 525 /* Create a ray from the tx to the rx */ 526 r = &raysInfo_los->rays[off]; 527 r->o = tx_pos[tx]; 528 r->d = vec3_sub(&rx_pos[rx], &r->o); 529 530 /* In case the tx and the rx are at the same position */ 531 if (vec3_dot(&r->d, &r->d) < __FLT_EPSILON__) { 532 /* Assume the direction to be (1, 0, 0) */ 533 chanInfo_los->directions_rx[off] = (Vec3){1.f, 0.f, 0.f}; 534 chanInfo_los->directions_tx[off] = (Vec3){-1.f, 0.f, 0.f}; 535 /* Set gains to 1+0j */ 536 chanInfo_los->a_te_re[off] = chanInfo_los->a_tm_re[off] = 1.f; 537 /* Set delay to 0 */ 538 chanInfo_los->tau[off] = 0.f; 539 /* Set doppler shift to 0. */ 540 chanInfo_los->freq_shift[off] = 0.f; 541 /* Set the active bit to 1 */ 542 raysInfo_los->rays_active[off / 8] |= 1 << (off % 8); 543 continue; 544 } 545 546 /* Is there any abstacle between the tx and the rx? */ 547 moeller_trumbore(r, scene, &t, &mesh_ind, &face_ind, &theta); 548 if (mesh_ind != (uint32_t)-1 && t <= 1.f) { 549 /* An obstacle between the tx and the rx has been hit */ 550 chanInfo_los->a_te_re[off] = chanInfo_los->a_tm_re[off] = chanInfo_los->tau[off] = 0.f; 551 /* Set the active bit to 0 */ 552 raysInfo_los->rays_active[off / 8] &= ~(1 << (off % 8)); 553 continue; 554 } 555 556 /* No obstacle has been hit, the path is LoS */ 557 /* Calculate the distance */ 558 t = sqrtf(vec3_dot(&r->d, &r->d)); 559 /* Calculate the direction */ 560 Vec3 d = {r->d.x / t, r->d.y / t, r->d.z / t}; 561 chanInfo_los->directions_tx[off] = d; 562 chanInfo_los->directions_rx[off] = (Vec3){-d.x, -d.y, -d.z}; 563 /* Calculate the free space loss */ 564 free_space_loss = free_space_loss_multiplier * t; 565 if (free_space_loss > 1.f) { 566 chanInfo_los->a_te_re[off] = 1.f / free_space_loss; 567 chanInfo_los->a_tm_re[off] = 1.f / free_space_loss; 568 } else 569 chanInfo_los->a_te_re[off] =chanInfo_los-> a_tm_re[off] = 1.f; 570 /* Calculate the delay */ 571 chanInfo_los->tau[off] = t / SPEED_OF_LIGHT; 572 /* Calculate the doppler shift */ 573 chanInfo_los->freq_shift[off] = vec3_dot(tx_vel, &d) - vec3_dot(rx_vel, &d); 574 chanInfo_los->freq_shift[off] *= carrier_frequency_hz / SPEED_OF_LIGHT; 575 /* Set the active bit to 1 */ 576 raysInfo_los->rays_active[off / 8] |= 1 << (off % 8); 577 } 578 579 /***************************************************************************/ 580 /* Scatter paths */ 581 /***************************************************************************/ 582 583 /* Bounce the rays. On each bounce: 584 * - Add path length / SPEED_OF_LIGHT to tau 585 * - Update gains using eqs. (31a)-(31b) from ITU-R P.2040-3 586 * - Spawn scattering rays towards the rx 587 * - TODO Spawn refraction rays with gains as per eqs. (31c)-(31d) from ITU-R P.2040-3 588 */ 589 raysInfo_scat->rays = (Ray*)memcpy(raysInfo_scat->rays, rays, num_tx * num_paths * sizeof(Ray)); 590 591 for (size_t bounce = 0; bounce < num_bounces; ++bounce) { 592 /* active[active_byte_index] & (1 << active_bit_pos) */ 593 size_t active_byte_index = 0; 594 uint8_t active_bit = 1; 595 size_t off_tx_path = 0; /* tx * num_paths + path */ 596 for (size_t tx = 0; tx < num_tx; ++tx) { 597 for (size_t path = 0; path < num_paths; ++path, ++off_tx_path, active_bit <<= 1) { 598 599 /* Check if the ray is active */ 600 if (!active_bit) { 601 active_bit = 1; 602 ++active_byte_index; 603 } 604 if (!(active[active_byte_index] & active_bit)) 605 continue; 606 607 /***********************************************************************/ 608 /* Reflect the rays */ 609 /***********************************************************************/ 610 611 /* Init */ 612 t = -1.f; 613 mesh_ind = face_ind = -1; 614 /* Find the hit point and trinagle and the angle of incidence */ 615 moeller_trumbore(&rays[off_tx_path], scene, &t, &mesh_ind, &face_ind, &theta); 616 if (mesh_ind == (uint32_t)-1) { /* Ray hit nothing */ 617 active[active_byte_index] &= ~active_bit; 618 --num_active; 619 continue; 620 } 621 /* Calculate the reflection coefficients R_{eTE} and R_{eTM} */ 622 material_ind = scene->meshes[mesh_ind].material_index; 623 refl_coefs(material_ind, theta, 624 &r_te_re, &r_te_im, 625 &r_tm_re, &r_tm_im); 626 /* Calculate the free space loss */ 627 free_space_loss = free_space_loss_multiplier * t; 628 free_space_loss *= free_space_loss; 629 if (free_space_loss > 1.f) { 630 r_te_re /= free_space_loss; 631 r_te_im /= free_space_loss; 632 r_tm_re /= free_space_loss; 633 r_tm_im /= free_space_loss; 634 } 635 /* Update the gains as a' = a * refl_coefs */ 636 float a_te_re_new = a_te_re[off_tx_path] * r_te_re - a_te_im[off_tx_path] * r_te_im; 637 float a_te_im_new = a_te_re[off_tx_path] * r_te_im + a_te_im[off_tx_path] * r_te_re; 638 float a_tm_re_new = a_tm_re[off_tx_path] * r_tm_re - a_tm_im[off_tx_path] * r_tm_im; 639 float a_tm_im_new = a_tm_re[off_tx_path] * r_tm_im + a_tm_im[off_tx_path] * r_tm_re; 640 a_te_re[off_tx_path] = a_te_re_new; 641 a_te_im[off_tx_path] = a_te_im_new; 642 a_tm_re[off_tx_path] = a_tm_re_new; 643 a_tm_im[off_tx_path] = a_tm_im_new; 644 /* Update the delay */ 645 tau[off_tx_path] += t / SPEED_OF_LIGHT; 646 647 /* Move and reflect the ray */ 648 r = &rays[off_tx_path]; 649 /* Advance the ray to the hit point */ 650 Vec3 h = vec3_scale(&r->d, t); 651 r->o = vec3_add(&h, &r->o); 652 /* Reflect the ray's direction as d' = d - 2*(d \cdot n)*n */ 653 Vec3 n = scene->meshes[mesh_ind].ns[face_ind]; 654 t = vec3_dot(&r->d, &n); 655 h = vec3_scale(&n, 2.f * t); 656 r->d = vec3_sub(&r->d, &h); 657 /* Advance the ray origin a bit to avoid intersection with the same triangle */ 658 h = vec3_scale(&r->d, 1e-4f); 659 r->o = vec3_add(&r->o, &h); 660 661 /* Calculate the doppler shift caused by the reflection */ 662 Vec3* mesh_vel = &scene->meshes[mesh_ind].velocity; 663 Vec3 d_rd = vec3_sub(&r->d, &rays[off_tx_path].d); 664 chanInfo_scat->freq_shift[off_tx_path] += vec3_dot(&d_rd, mesh_vel) * doppler_shift_multiplier; 665 666 /***********************************************************************/ 667 /* Scatter the rays */ 668 /***********************************************************************/ 669 670 /* Scatter a path from the hit point towards the rx */ 671 Ray r_scat = { .o = r->o }; 672 for (size_t rx = 0; rx < num_rx; ++rx) { 673 /* [rx, tx, bounce, path] */ 674 size_t off_scat = ((rx * num_tx + tx) * num_bounces + bounce) * num_paths + path; 675 /* Create a ray from the hit point to the rx */ 676 r_scat.d = vec3_sub(&rx_pos[rx], &r_scat.o); 677 float d2rx = sqrtf(vec3_dot(&r_scat.d, &r_scat.d)); 678 r_scat.d = vec3_normalize(&r_scat.d); 679 /* Check if the ray has a LoS of the rx */ 680 t = -1.f; 681 mesh_ind = face_ind = -1; 682 moeller_trumbore(&r_scat, scene, &t, &mesh_ind, &face_ind, &theta); 683 if (mesh_ind != (uint32_t)-1 && t <= 1.f) { 684 /* An obstacle between the hit point and the rx has been hit */ 685 chanInfo_scat->a_te_re[off_scat] = chanInfo_scat->a_te_im[off_scat] 686 = chanInfo_scat->a_tm_re[off_scat] 687 = chanInfo_scat->a_tm_im[off_scat] 688 = chanInfo_scat->tau[off_scat] 689 = 0.f; 690 continue; 691 } 692 /* No obstacle has been hit */ 693 /* Calculate the scattering angle */ 694 theta_s = acosf(vec3_dot(&r_scat.d, &n)); 695 /* Calculate the scattering coefficients */ 696 scat_coefs(theta_s, theta, material_ind, 697 &a_te_re_new, &a_te_im_new, &a_tm_re_new, &a_tm_im_new); 698 chanInfo_scat->a_te_re[off_scat] = a_te_re[off_tx_path] * a_te_re_new 699 - a_te_im[off_tx_path] * a_te_im_new; 700 chanInfo_scat->a_te_im[off_scat] = a_te_re[off_tx_path] * a_te_im_new 701 + a_te_im[off_tx_path] * a_te_re_new; 702 chanInfo_scat->a_tm_re[off_scat] = a_tm_re[off_tx_path] * a_tm_re_new 703 - a_tm_im[off_tx_path] * a_tm_im_new; 704 chanInfo_scat->a_tm_im[off_scat] = a_tm_re[off_tx_path] * a_tm_im_new 705 + a_tm_im[off_tx_path] * a_tm_re_new; 706 /* Calculate the direction */ 707 chanInfo_scat->directions_rx[off_scat] = (Vec3){-r_scat.d.x, -r_scat.d.y, -r_scat.d.z}; 708 /* Calculate the delay */ 709 chanInfo_scat->tau[off_scat] = tau[off_tx_path] + d2rx / SPEED_OF_LIGHT; 710 /* Calculate the free space loss */ 711 free_space_loss = free_space_loss_multiplier * d2rx; 712 free_space_loss *= free_space_loss; 713 if (free_space_loss > 1.f) { 714 chanInfo_scat->a_te_re[off_scat] /= free_space_loss; 715 chanInfo_scat->a_te_im[off_scat] /= free_space_loss; 716 chanInfo_scat->a_tm_re[off_scat] /= free_space_loss; 717 chanInfo_scat->a_tm_im[off_scat] /= free_space_loss; 718 } 719 /* Calculate the doppler shift */ 720 Vec3 d_rd = vec3_sub(&r_scat.d, &rays[off_tx_path].d); 721 float freq_shift = vec3_dot(&d_rd, mesh_vel) * doppler_shift_multiplier; 722 chanInfo_scat->freq_shift[off_scat] -= freq_shift; 723 } 724 725 /***********************************************************************/ 726 /* Refract the rays */ 727 /***********************************************************************/ 728 /* TODO */ 729 } 730 731 /* Copy the rays to the raysInfo_scat */ 732 size_t off_scatter_rays = (tx * num_bounces + (bounce + 1)) * num_paths; 733 size_t off_active = (tx * num_bounces + (bounce + 1)) * (num_paths / 8 + 1); 734 memcpy( 735 raysInfo_scat->rays + off_scatter_rays, 736 rays + tx * num_paths, 737 num_paths * sizeof(Ray) 738 ); 739 memcpy( 740 raysInfo_scat->rays_active + off_active, 741 active, 742 (num_paths / 8 + 1) * sizeof(uint8_t) 743 ); 744 } 745 } 746 747 /* Free */ 748 free(a_te_re); 749 free(a_te_im); 750 free(a_tm_re); 751 free(a_tm_im); 752 free(tau); 753 free(active); 754 free(rays); 755 /* Free the scene */ 756 /* TODO */ 757 }