ssol_mc_receiver.c (6738B)
1 /* Copyright (C) 2018-2026 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2016, 2018 CNRS 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #include "ssol.h" 18 #include "ssol_estimator_c.h" 19 #include "ssol_object_c.h" 20 21 #include <rsys/double3.h> 22 #include <star/s3d.h> 23 24 /******************************************************************************* 25 * Exported functions 26 ******************************************************************************/ 27 res_T 28 ssol_estimator_get_mc_receiver 29 (struct ssol_estimator* estimator, 30 const struct ssol_instance* instance, 31 const enum ssol_side_flag side, 32 struct ssol_mc_receiver* rcv) 33 { 34 struct mc_receiver* mc_rcv = NULL; 35 struct mc_receiver_1side* mc_rcv1 = NULL; 36 37 if(!estimator || !instance || !rcv 38 || !(instance->receiver_mask & (int)side)) 39 return RES_BAD_ARG; 40 41 memset(rcv, 0, sizeof(rcv[0])); 42 43 mc_rcv = htable_receiver_find(&estimator->mc_receivers, &instance); 44 if(!mc_rcv) { 45 /* The receiver has no MC estimation */ 46 return RES_OK; 47 } 48 49 mc_rcv1 = side == SSOL_FRONT ? &mc_rcv->front : &mc_rcv->back; 50 #define SETUP_MC_RESULT(Name) { \ 51 const double N = (double)estimator->realisation_count; \ 52 struct mc_data* data = &mc_rcv1->Name; \ 53 double weight, sqr_weight; \ 54 mc_data_get(data, &weight, &sqr_weight); \ 55 rcv->Name.E = weight / N; \ 56 rcv->Name.V = sqr_weight/N - rcv->Name.E*rcv->Name.E; \ 57 rcv->Name.V = rcv->Name.V > 0 ? rcv->Name.V : 0; \ 58 rcv->Name.SE = sqrt(rcv->Name.V / N); \ 59 } (void)0 60 #define MC_SETUP_ALL { \ 61 SETUP_MC_RESULT(incoming_flux); \ 62 SETUP_MC_RESULT(incoming_if_no_atm_loss); \ 63 SETUP_MC_RESULT(incoming_if_no_field_loss); \ 64 SETUP_MC_RESULT(incoming_lost_in_atmosphere); \ 65 SETUP_MC_RESULT(incoming_lost_in_field); \ 66 SETUP_MC_RESULT(absorbed_flux); \ 67 SETUP_MC_RESULT(absorbed_if_no_atm_loss); \ 68 SETUP_MC_RESULT(absorbed_if_no_field_loss); \ 69 SETUP_MC_RESULT(absorbed_lost_in_atmosphere); \ 70 SETUP_MC_RESULT(absorbed_lost_in_field); \ 71 } (void)0 72 MC_SETUP_ALL; 73 #undef SETUP_MC_RESULT 74 rcv->mc__ = mc_rcv1; 75 rcv->N__ = estimator->realisation_count; 76 rcv->instance__ = instance; 77 return RES_OK; 78 } 79 80 res_T 81 ssol_mc_receiver_get_mc_shape 82 (struct ssol_mc_receiver* rcv, 83 const struct ssol_shape* shape, 84 struct ssol_mc_shape* mc) 85 { 86 struct mc_receiver_1side* mc_rcv1; 87 88 if(!rcv || !shape || !mc) return RES_BAD_ARG; 89 if(!object_has_shape(rcv->instance__->object, shape)) return RES_BAD_ARG; 90 mc_rcv1 = rcv->mc__; 91 mc->N__ = rcv->N__; 92 mc->mc__ = htable_shape2mc_find(&mc_rcv1->shape2mc, &shape); 93 mc->shape__ = shape; 94 return RES_OK; 95 } 96 97 res_T 98 ssol_mc_shape_get_mc_primitive 99 (struct ssol_mc_shape* shape, 100 const unsigned i, 101 struct ssol_mc_primitive* prim) 102 { 103 struct mc_shape_1side* mc_shape1; 104 struct mc_primitive_1side* mc_prim1; 105 unsigned ntris; 106 107 if(!shape || !prim) return RES_BAD_ARG; 108 109 SSOL(shape_get_triangles_count(shape->shape__, &ntris)); 110 if(i >= ntris) return RES_BAD_ARG; 111 112 mc_shape1 = shape->mc__; 113 if(!mc_shape1 || !(mc_prim1 = htable_prim2mc_find(&mc_shape1->prim2mc, &i))) { 114 #define SETUP_MC_RESULT(Name) { \ 115 prim->Name.E = 0; \ 116 prim->Name.V = 0; \ 117 prim->Name.SE = 0; \ 118 } (void)0 119 MC_SETUP_ALL; 120 #undef SETUP_MC_RESULT 121 } else { 122 struct s3d_attrib attr; 123 struct s3d_shape* s3d_shape; 124 double v0[3], v1[3], v2[3], E0[3], E1[3], normal[3]; 125 double area; 126 unsigned ids[3]; 127 res_T res = RES_OK; 128 129 s3d_shape = shape->shape__->shape_rt; 130 131 /* Retrieve the primitive indices */ 132 res = s3d_mesh_get_triangle_indices(s3d_shape, i, ids); 133 if(res != RES_OK) return res; 134 135 /* Fetch the primitive vertices */ 136 S3D(mesh_get_vertex_attrib(s3d_shape, ids[0], S3D_POSITION, &attr)); 137 d3_set_f3(v0, attr.value); 138 S3D(mesh_get_vertex_attrib(s3d_shape, ids[1], S3D_POSITION, &attr)); 139 d3_set_f3(v1, attr.value); 140 S3D(mesh_get_vertex_attrib(s3d_shape, ids[2], S3D_POSITION, &attr)); 141 d3_set_f3(v2, attr.value); 142 143 /* Compute the primitive area */ 144 d3_sub(E0, v1, v0); 145 d3_sub(E1, v2, v0); 146 d3_cross(normal, E0, E1); 147 area = d3_len(normal) * 0.5; 148 149 #define SETUP_MC_RESULT(Name) { \ 150 const double N = (double)shape->N__; \ 151 struct mc_data* data = &mc_prim1->Name; \ 152 double weight, sqr_weight; \ 153 mc_data_get(data, &weight, &sqr_weight); \ 154 prim->Name.E = weight / N; \ 155 prim->Name.V = sqr_weight/N - prim->Name.E*prim->Name.E; \ 156 prim->Name.V = prim->Name.V > 0 ? prim->Name.V : 0; \ 157 prim->Name.SE = sqrt(prim->Name.V / N); \ 158 prim->Name.E /= area; \ 159 prim->Name.V /= area*area; \ 160 prim->Name.SE /= area; \ 161 } (void)0 162 MC_SETUP_ALL; 163 #undef SETUP_MC_RESULT 164 #undef MC_SETUP_ALL 165 } 166 167 return RES_OK; 168 }