test_htrdr_combustion_laser.c (5144B)
1 /* Copyright (C) 2018-2019, 2022-2025 Centre National de la Recherche Scientifique 2 * Copyright (C) 2020-2022 Institut Mines Télécom Albi-Carmaux 3 * Copyright (C) 2022-2025 Institut Pierre-Simon Laplace 4 * Copyright (C) 2022-2025 Institut de Physique du Globe de Paris 5 * Copyright (C) 2018-2025 |Méso|Star> (contact@meso-star.com) 6 * Copyright (C) 2022-2025 Observatoire de Paris 7 * Copyright (C) 2022-2025 Université de Reims Champagne-Ardenne 8 * Copyright (C) 2022-2025 Université de Versaille Saint-Quentin 9 * Copyright (C) 2018-2019, 2022-2025 Université Paul Sabatier 10 * 11 * This program is free software: you can redistribute it and/or modify 12 * it under the terms of the GNU General Public License as published by 13 * the Free Software Foundation, either version 3 of the License, or 14 * (at your option) any later version. 15 * 16 * This program is distributed in the hope that it will be useful, 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 * GNU General Public License for more details. 20 * 21 * You should have received a copy of the GNU General Public License 22 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 23 24 #include "combustion/htrdr_combustion_laser.h" 25 26 #include <rsys/double2.h> 27 #include <rsys/double3.h> 28 29 static void 30 dump_obj(const struct htrdr_combustion_laser_mesh* mesh, FILE* stream) 31 { 32 unsigned i; 33 ASSERT(mesh && stream); 34 35 FOR_EACH(i, 0, mesh->nvertices) { 36 fprintf(stream, "v %g %g %g\n", 37 mesh->vertices[i*3+0], 38 mesh->vertices[i*3+1], 39 mesh->vertices[i*3+2]); 40 } 41 FOR_EACH(i, 0, mesh->ntriangles) { 42 fprintf(stream, "f %u %u %u\n", 43 mesh->triangles[i*3+0]+1, 44 mesh->triangles[i*3+1]+1, 45 mesh->triangles[i*3+2]+1); 46 } 47 } 48 49 int 50 main(int argc, char** argv) 51 { 52 struct htrdr_args args = HTRDR_ARGS_DEFAULT; 53 struct htrdr_combustion_laser_mesh mesh; 54 struct htrdr_combustion_laser_create_args laser_args = 55 HTRDR_COMBUSTION_LASER_CREATE_ARGS_DEFAULT; 56 struct htrdr* htrdr = NULL; 57 struct htrdr_combustion_laser* laser = NULL; 58 double org[3]; 59 double dir[3]; 60 double range[2]; 61 double hit_range[2]; 62 double t[2]; 63 double pt[3]; 64 double x[3], y[3], z[3]; 65 double plane0[4]; 66 double plane1[4]; 67 FILE* fp = NULL; 68 69 args.verbose = 1; 70 htrdr_mpi_init(argc, argv); 71 CHK(htrdr_create(NULL, &args, &htrdr) == RES_OK); 72 73 /* Setup the laser sheet */ 74 d3(laser_args.surface.position, 0, 0, 0); 75 d3(laser_args.surface.target, 10, 10, 0); 76 d3(laser_args.surface.up, 0, 1, 0); 77 d2(laser_args.surface.size, 100, 50); 78 laser_args.wavelength = 300; 79 laser_args.flux_density = 1; 80 CHK(htrdr_combustion_laser_create(htrdr, &laser_args, &laser) == RES_OK); 81 htrdr_combustion_laser_get_mesh(laser, 100/*arbitrary extend*/, &mesh); 82 83 /* Write the laser geometry */ 84 CHK(fp = fopen("laser.obj", "w")); 85 dump_obj(&mesh, fp); 86 fclose(fp); 87 88 /* Compute the frame of the surface emission */ 89 d3_sub(z, laser_args.surface.target, laser_args.surface.position); 90 d3_cross(x, z, laser_args.surface.up); 91 d3_cross(y, x, z); 92 CHK(d3_normalize(y, y) != 0); 93 94 /* Compute the bottom plane equation of the laser sheet */ 95 pt[0] = laser_args.surface.position[0] - y[0]*laser_args.surface.size[1]*0.5; 96 pt[1] = laser_args.surface.position[1] - y[1]*laser_args.surface.size[1]*0.5; 97 pt[2] = laser_args.surface.position[2] - y[2]*laser_args.surface.size[1]*0.5; 98 plane0[0] = y[0]; 99 plane0[1] = y[1]; 100 plane0[2] = y[2]; 101 plane0[3] = -d3_dot(y, pt); 102 103 /* Compute the top plane equation of the laser sheet */ 104 pt[0] = laser_args.surface.position[0] + y[0]*laser_args.surface.size[1]*0.5; 105 pt[1] = laser_args.surface.position[1] + y[1]*laser_args.surface.size[1]*0.5; 106 pt[2] = laser_args.surface.position[2] + y[2]*laser_args.surface.size[1]*0.5; 107 plane1[0] = y[0]; 108 plane1[1] = y[1]; 109 plane1[2] = y[2]; 110 plane1[3] = -d3_dot(y, pt); 111 112 /* Trace a ray that misses the laser sheet */ 113 d3(org, 50, 0, 0); 114 d3(dir, 0, -1, 0); 115 d2(range, 0, INF); 116 htrdr_combustion_laser_trace_ray(laser, org, dir, range, hit_range); 117 CHK(hit_range[0] > DBL_MAX); 118 CHK(hit_range[1] > DBL_MAX); 119 120 /* Trace a ray that intersects both bottom and top laser planes */ 121 d3(dir, 0, 1, 0); 122 htrdr_combustion_laser_trace_ray(laser, org, dir, range, hit_range); 123 CHK(hit_range[0] < hit_range[1]); 124 125 /* Compute the intersection of the ray with the bottom/top laser planes */ 126 t[0] = (-d3_dot(plane0, org) - plane0[3]) / d3_dot(plane0, dir); 127 t[1] = (-d3_dot(plane1, org) - plane1[3]) / d3_dot(plane1, dir); 128 129 /* Check the returned distances against the computed ones */ 130 CHK(eq_eps(hit_range[0], t[0], 1.e-6*hit_range[0])); 131 CHK(eq_eps(hit_range[1], t[1], 1.e-6*hit_range[1])); 132 133 /* Trace a ray that starts into the laser sheet */ 134 range[0] = 0.5*(hit_range[0] + hit_range[1]); 135 htrdr_combustion_laser_trace_ray(laser, org, dir, range, hit_range); 136 CHK(hit_range[0] < hit_range[1]); 137 CHK(hit_range[0] == range[0]); 138 CHK(eq_eps(hit_range[1], t[1], 1.e-6*hit_range[1])); 139 140 htrdr_ref_put(htrdr); 141 htrdr_combustion_laser_ref_put(laser); 142 htrdr_mpi_finalize(); 143 144 return 0; 145 } 146