From cc97c6d667e91455b8730bb47ddb8f37e09161e4 Mon Sep 17 00:00:00 2001 From: zoeprieto Date: Mon, 7 Oct 2024 14:14:44 -0500 Subject: [PATCH 1/3] format check test added --- .github/workflows/format-check.yml | 31 ++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) create mode 100644 .github/workflows/format-check.yml diff --git a/.github/workflows/format-check.yml b/.github/workflows/format-check.yml new file mode 100644 index 0000000..159ea72 --- /dev/null +++ b/.github/workflows/format-check.yml @@ -0,0 +1,31 @@ +name: C++ Format Check + +on: + # allow workflow to be run manually + workflow_dispatch: + + pull_request: + branches: + - develop + - master + - format_C + +jobs: + cpp-linter: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: cpp-linter/cpp-linter-action@v2 + id: linter + with: + style: file + files-changed-only: true + tidy-checks: '-*' + version: '15' # clang-format version + file-annotations: true + step-summary: true + extensions: 'c,h' + + - name: Failure Check + if: steps.linter.outputs.checks-failed > 0 + run: echo "Some files failed the formatting check! See job summary and file annotations for more info" && exit 1 \ No newline at end of file From 1e0e44407171f3481bd12208841921879ac0967e Mon Sep 17 00:00:00 2001 From: zoeprieto Date: Mon, 7 Oct 2024 14:20:05 -0500 Subject: [PATCH 2/3] test format --- src/kdsource/geom.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/kdsource/geom.c b/src/kdsource/geom.c index a7b9ec3..5aa9e01 100644 --- a/src/kdsource/geom.c +++ b/src/kdsource/geom.c @@ -158,7 +158,9 @@ int Let_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char ker return 0; } -int wl_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ +int wl_perturb + + (const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ double wl = 9.045 / sqrt(part->ekin * 1e9); double wl2 = wl + bw*metric->scaling[0] * rand_type(kernel); part->ekin = 81.82e-9 / (wl2 * wl2); From 0ebfa8bff4a0ac809cd6bc0e9432d2ce19784141 Mon Sep 17 00:00:00 2001 From: zoeprieto Date: Mon, 7 Oct 2024 14:25:05 -0500 Subject: [PATCH 3/3] format fix --- .github/workflows/format-check.yml | 2 +- src/kdsource/geom.c | 1003 ++++++++++++++++------------ src/kdsource/geom.h | 129 ++-- src/kdsource/kdsource.c | 556 +++++++-------- src/kdsource/kdsource.h | 77 +-- src/kdsource/plist.c | 183 ++--- src/kdsource/plist.h | 45 +- src/kdsource/utils.c | 255 ++++--- src/kdsource/utils.h | 12 +- src/kdtool/beamtest.c | 406 +++++------ src/kdtool/resample.c | 144 ++-- 11 files changed, 1553 insertions(+), 1259 deletions(-) diff --git a/.github/workflows/format-check.yml b/.github/workflows/format-check.yml index 159ea72..a9d3d9f 100644 --- a/.github/workflows/format-check.yml +++ b/.github/workflows/format-check.yml @@ -1,4 +1,4 @@ -name: C++ Format Check +name: C Format Check on: # allow workflow to be run manually diff --git a/src/kdsource/geom.c b/src/kdsource/geom.c index 5aa9e01..723d0ad 100644 --- a/src/kdsource/geom.c +++ b/src/kdsource/geom.c @@ -1,485 +1,606 @@ -#include -#include -#include +#include +#include +#include -#include +#include #include "kdsource.h" - -void KDS_error(const char* msg); -void KDS_end(const char* msg); - -Metric* Metric_create(int dim, const double* scaling, PerturbFun perturb, int nps, const double* params){ - Metric* metric = (Metric*)malloc(sizeof(Metric)); - int i; - metric->dim = dim; - metric->scaling = (float*)malloc(dim*sizeof(float)); - if(scaling) for(i=0; iscaling[i] = (float)scaling[i]; - else for(i=0; iscaling[i] = 0; - metric->perturb = perturb; - metric->nps = nps; - metric->params = (double*)malloc(nps * sizeof(double)); - for(i=0; iparams[i] = params[i]; - return metric; +void KDS_error(const char *msg); +void KDS_end(const char *msg); + +Metric *Metric_create(int dim, const double *scaling, PerturbFun perturb, + int nps, const double *params) { + Metric *metric = (Metric *)malloc(sizeof(Metric)); + int i; + metric->dim = dim; + metric->scaling = (float *)malloc(dim * sizeof(float)); + if (scaling) + for (i = 0; i < dim; i++) + metric->scaling[i] = (float)scaling[i]; + else + for (i = 0; i < dim; i++) + metric->scaling[i] = 0; + metric->perturb = perturb; + metric->nps = nps; + metric->params = (double *)malloc(nps * sizeof(double)); + for (i = 0; i < nps; i++) + metric->params[i] = params[i]; + return metric; } -Metric* Metric_copy(const Metric* from){ - Metric* metric = (Metric*)malloc(sizeof(Metric)); - *metric = *from; - int i; - metric->scaling = (float*)malloc(metric->dim*sizeof(double)); - for(i=0; idim; i++) metric->scaling[i] = from->scaling[i]; - metric->params = (double*)malloc(metric->nps * sizeof(double)); - for(i=0; inps; i++) metric->params[i] = from->params[i]; - return metric; +Metric *Metric_copy(const Metric *from) { + Metric *metric = (Metric *)malloc(sizeof(Metric)); + *metric = *from; + int i; + metric->scaling = (float *)malloc(metric->dim * sizeof(double)); + for (i = 0; i < metric->dim; i++) + metric->scaling[i] = from->scaling[i]; + metric->params = (double *)malloc(metric->nps * sizeof(double)); + for (i = 0; i < metric->nps; i++) + metric->params[i] = from->params[i]; + return metric; } -void Metric_destroy(Metric* metric){ - free(metric->scaling); - free(metric); +void Metric_destroy(Metric *metric) { + free(metric->scaling); + free(metric); } -Geometry* Geom_create(int ord, Metric** metrics, double bw, const char* bwfilename, char kernel, - const double* trasl, const double* rot){ - Geometry* geom = (Geometry*)malloc(sizeof(Geometry)); - geom->ord = ord; - int i; - geom->ms = (Metric**)malloc(ord * sizeof(Metric*)); - for(i=0; ims[i] = metrics[i]; - geom->bw = bw; - geom->bwfilename = NULL; - geom->bwfile = NULL; - geom->kernel = kernel; - if(bwfilename) if(strlen(bwfilename)){ - FILE* bwfile; - if((bwfile=fopen(bwfilename, "rb")) == 0){ - printf("Could not open file %s\n", bwfilename); - KDS_error("Error in Geom_create"); - } - geom->bwfilename = (char*)malloc(NAME_MAX_LEN*sizeof(char)); - strcpy(geom->bwfilename, bwfilename); - geom->bwfile = bwfile; - Geom_next(geom, 1); - } - if(trasl){ - geom->trasl = (double*)malloc(3 * sizeof(double)); - for(i=0; i<3; i++) geom->trasl[i] = trasl[i]; - } - else geom->trasl = NULL; - if(rot){ - geom->rot = (double*)malloc(3 * sizeof(double)); - for(i=0; i<3; i++) geom->rot[i] = rot[i]; - } - else geom->rot = NULL; - return geom; +Geometry *Geom_create(int ord, Metric **metrics, double bw, + const char *bwfilename, char kernel, const double *trasl, + const double *rot) { + Geometry *geom = (Geometry *)malloc(sizeof(Geometry)); + geom->ord = ord; + int i; + geom->ms = (Metric **)malloc(ord * sizeof(Metric *)); + for (i = 0; i < ord; i++) + geom->ms[i] = metrics[i]; + geom->bw = bw; + geom->bwfilename = NULL; + geom->bwfile = NULL; + geom->kernel = kernel; + if (bwfilename) + if (strlen(bwfilename)) { + FILE *bwfile; + if ((bwfile = fopen(bwfilename, "rb")) == 0) { + printf("Could not open file %s\n", bwfilename); + KDS_error("Error in Geom_create"); + } + geom->bwfilename = (char *)malloc(NAME_MAX_LEN * sizeof(char)); + strcpy(geom->bwfilename, bwfilename); + geom->bwfile = bwfile; + Geom_next(geom, 1); + } + if (trasl) { + geom->trasl = (double *)malloc(3 * sizeof(double)); + for (i = 0; i < 3; i++) + geom->trasl[i] = trasl[i]; + } else + geom->trasl = NULL; + if (rot) { + geom->rot = (double *)malloc(3 * sizeof(double)); + for (i = 0; i < 3; i++) + geom->rot[i] = rot[i]; + } else + geom->rot = NULL; + return geom; } -Geometry* Geom_copy(const Geometry* from){ - Geometry* geom = (Geometry*)malloc(sizeof(Geometry)); - *geom = *from; - int i; - geom->ms = (Metric**)malloc(geom->ord * sizeof(Metric*)); - for(i=0; iord; i++) geom->ms[i] = Metric_copy(from->ms[i]); - geom->bwfilename = NULL; - geom->bwfile = NULL; - if(geom->bwfile){ - geom->bwfilename = (char*)malloc(NAME_MAX_LEN*sizeof(char)); - strcpy(geom->bwfilename, from->bwfilename); - geom->bwfile = fopen(geom->bwfilename, "rb"); - Geom_next(geom, 1); - } - if(from->trasl){ - geom->trasl = (double*)malloc(3 * sizeof(double)); - for(i=0; i<3; i++) geom->trasl[i] = from->trasl[i]; - } - else geom->trasl = NULL; - if(from->rot){ - geom->rot = (double*)malloc(3 * sizeof(double)); - for(i=0; i<3; i++) geom->rot[i] = from->rot[i]; - } - else geom->rot = NULL; - return geom; +Geometry *Geom_copy(const Geometry *from) { + Geometry *geom = (Geometry *)malloc(sizeof(Geometry)); + *geom = *from; + int i; + geom->ms = (Metric **)malloc(geom->ord * sizeof(Metric *)); + for (i = 0; i < geom->ord; i++) + geom->ms[i] = Metric_copy(from->ms[i]); + geom->bwfilename = NULL; + geom->bwfile = NULL; + if (geom->bwfile) { + geom->bwfilename = (char *)malloc(NAME_MAX_LEN * sizeof(char)); + strcpy(geom->bwfilename, from->bwfilename); + geom->bwfile = fopen(geom->bwfilename, "rb"); + Geom_next(geom, 1); + } + if (from->trasl) { + geom->trasl = (double *)malloc(3 * sizeof(double)); + for (i = 0; i < 3; i++) + geom->trasl[i] = from->trasl[i]; + } else + geom->trasl = NULL; + if (from->rot) { + geom->rot = (double *)malloc(3 * sizeof(double)); + for (i = 0; i < 3; i++) + geom->rot[i] = from->rot[i]; + } else + geom->rot = NULL; + return geom; } -int Geom_perturb(const Geometry* geom, mcpl_particle_t* part){ - int i, ret=0; - if(geom->trasl) traslv(part->position, geom->trasl, 1); - if(geom->rot){ rotv(part->position, geom->rot, 1); rotv(part->direction, geom->rot, 1); } - for(i=0; iord; i++) - ret += geom->ms[i]->perturb(geom->ms[i], part, geom->bw, geom->kernel); - if(geom->rot){ rotv(part->position, geom->rot, 0); rotv(part->direction, geom->rot, 0); } - if(geom->trasl) traslv(part->position, geom->trasl, 0); - return ret; +int Geom_perturb(const Geometry *geom, mcpl_particle_t *part) { + int i, ret = 0; + if (geom->trasl) + traslv(part->position, geom->trasl, 1); + if (geom->rot) { + rotv(part->position, geom->rot, 1); + rotv(part->direction, geom->rot, 1); + } + for (i = 0; i < geom->ord; i++) + ret += geom->ms[i]->perturb(geom->ms[i], part, geom->bw, geom->kernel); + if (geom->rot) { + rotv(part->position, geom->rot, 0); + rotv(part->direction, geom->rot, 0); + } + if (geom->trasl) + traslv(part->position, geom->trasl, 0); + return ret; } -int Geom_next(Geometry* geom, int loop){ - if(geom->bwfile){ - float temp; - if(fread(&temp, sizeof(float), 1, geom->bwfile) == 1){ - geom->bw = temp; - return 0; - } - if(!loop) KDS_end("Geom_next: End of BW file reached."); - rewind(geom->bwfile); // After 1st failed try, rewind - if(fread(&temp, sizeof(float), 1, geom->bwfile) == 1){ - geom->bw = temp; - return 1; - } - KDS_error("Error in Geom_next: Could not read BW."); - } - return 0; +int Geom_next(Geometry *geom, int loop) { + if (geom->bwfile) { + float temp; + if (fread(&temp, sizeof(float), 1, geom->bwfile) == 1) { + geom->bw = temp; + return 0; + } + if (!loop) + KDS_end("Geom_next: End of BW file reached."); + rewind(geom->bwfile); // After 1st failed try, rewind + if (fread(&temp, sizeof(float), 1, geom->bwfile) == 1) { + geom->bw = temp; + return 1; + } + KDS_error("Error in Geom_next: Could not read BW."); + } + return 0; } -void Geom_destroy(Geometry* geom){ - int i; - for(i=0; iord; i++) Metric_destroy(geom->ms[i]); - free(geom->ms); - free(geom->trasl); - free(geom->rot); - free(geom); +void Geom_destroy(Geometry *geom) { + int i; + for (i = 0; i < geom->ord; i++) + Metric_destroy(geom->ms[i]); + free(geom->ms); + free(geom->trasl); + free(geom->rot); + free(geom); } - -int E_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - part->ekin += bw*metric->scaling[0] * rand_type(kernel); - if(part->ekin < 0) part->ekin *= -1; - return 0; +int E_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + part->ekin += bw * metric->scaling[0] * rand_type(kernel); + if (part->ekin < 0) + part->ekin *= -1; + return 0; } -int Let_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - float E = part->ekin; - E *= exp(bw*metric->scaling[0] * rand_type(kernel)); - while(E > metric->params[0]){ - E = part->ekin; - E *= exp(bw*metric->scaling[0] * rand_type(kernel)); - } - part->ekin = E; - return 0; +int Let_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + float E = part->ekin; + E *= exp(bw * metric->scaling[0] * rand_type(kernel)); + while (E > metric->params[0]) { + E = part->ekin; + E *= exp(bw * metric->scaling[0] * rand_type(kernel)); + } + part->ekin = E; + return 0; } int wl_perturb - - (const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - double wl = 9.045 / sqrt(part->ekin * 1e9); - double wl2 = wl + bw*metric->scaling[0] * rand_type(kernel); - part->ekin = 81.82e-9 / (wl2 * wl2); - if(part->ekin < 0) part->ekin *= -1; - return 0; + + (const Metric *metric, mcpl_particle_t *part, double bw, char kernel) { + double wl = 9.045 / sqrt(part->ekin * 1e9); + double wl2 = wl + bw * metric->scaling[0] * rand_type(kernel); + part->ekin = 81.82e-9 / (wl2 * wl2); + if (part->ekin < 0) + part->ekin *= -1; + return 0; } -int t_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - part->time += bw*metric->scaling[0] * rand_type(kernel); - return 0; +int t_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + part->time += bw * metric->scaling[0] * rand_type(kernel); + return 0; } -int Dec_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - part->time *= pow(10, bw*metric->scaling[0] * rand_type(kernel)); - return 0; +int Dec_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + part->time *= pow(10, bw * metric->scaling[0] * rand_type(kernel)); + return 0; } -int Vol_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - double xmin=metric->params[0], xmax=metric->params[1]; - double ymin=metric->params[2], ymax=metric->params[3]; - double zmin=metric->params[4], zmax=metric->params[5]; - part->position[0] += bw*metric->scaling[0] * rand_type(kernel); - while((part->position[0] < xmin) || (part->position[0] > xmax)){ - if(part->position[0] < xmin) part->position[0] += 2 * (xmin - part->position[0]); - else part->position[0] -= 2 * (part->position[0] - xmax); - } - part->position[1] += bw*metric->scaling[1] * rand_type(kernel); - while((part->position[1] < ymin) || (part->position[1] > ymax)){ - if(part->position[1] < ymin) part->position[1] += 2 * (ymin - part->position[1]); - else part->position[1] -= 2 * (part->position[1] - ymax); - } - part->position[2] += bw*metric->scaling[2] * rand_type(kernel); - while((part->position[2] < zmin) || (part->position[2] > zmax)){ - if(part->position[2] < zmin) part->position[2] += 2 * (zmin - part->position[2]); - else part->position[2] -= 2 * (part->position[2] - zmax); - } - return 0; +int Vol_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + double xmin = metric->params[0], xmax = metric->params[1]; + double ymin = metric->params[2], ymax = metric->params[3]; + double zmin = metric->params[4], zmax = metric->params[5]; + part->position[0] += bw * metric->scaling[0] * rand_type(kernel); + while ((part->position[0] < xmin) || (part->position[0] > xmax)) { + if (part->position[0] < xmin) + part->position[0] += 2 * (xmin - part->position[0]); + else + part->position[0] -= 2 * (part->position[0] - xmax); + } + part->position[1] += bw * metric->scaling[1] * rand_type(kernel); + while ((part->position[1] < ymin) || (part->position[1] > ymax)) { + if (part->position[1] < ymin) + part->position[1] += 2 * (ymin - part->position[1]); + else + part->position[1] -= 2 * (part->position[1] - ymax); + } + part->position[2] += bw * metric->scaling[2] * rand_type(kernel); + while ((part->position[2] < zmin) || (part->position[2] > zmax)) { + if (part->position[2] < zmin) + part->position[2] += 2 * (zmin - part->position[2]); + else + part->position[2] -= 2 * (part->position[2] - zmax); + } + return 0; } -int SurfXY_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - double xmin=metric->params[0], xmax=metric->params[1]; - double ymin=metric->params[2], ymax=metric->params[3]; - part->position[0] += bw*metric->scaling[0] * rand_type(kernel); - while((part->position[0] < xmin) || (part->position[0] > xmax)){ - if(part->position[0] < xmin) part->position[0] += 2 * (xmin - part->position[0]); - else part->position[0] -= 2 * (part->position[0] - xmax); - } - part->position[1] += bw*metric->scaling[1] * rand_type(kernel); - while((part->position[1] < ymin) || (part->position[1] > ymax)){ - if(part->position[1] < ymin) part->position[1] += 2 * (ymin - part->position[1]); - else part->position[1] -= 2 * (part->position[1] - ymax); - } - return 0; +int SurfXY_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + double xmin = metric->params[0], xmax = metric->params[1]; + double ymin = metric->params[2], ymax = metric->params[3]; + part->position[0] += bw * metric->scaling[0] * rand_type(kernel); + while ((part->position[0] < xmin) || (part->position[0] > xmax)) { + if (part->position[0] < xmin) + part->position[0] += 2 * (xmin - part->position[0]); + else + part->position[0] -= 2 * (part->position[0] - xmax); + } + part->position[1] += bw * metric->scaling[1] * rand_type(kernel); + while ((part->position[1] < ymin) || (part->position[1] > ymax)) { + if (part->position[1] < ymin) + part->position[1] += 2 * (ymin - part->position[1]); + else + part->position[1] -= 2 * (part->position[1] - ymax); + } + return 0; } -int SurfR_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - double rho_min=metric->params[0], rho_max=metric->params[1]; - double psi_min=metric->params[2], psi_max=metric->params[3]; - double rho, psi, rho2, psi2; - - rho = sqrt(part->position[0]*part->position[0]+part->position[1]*part->position[1]); - psi = atan2(part->position[1], part->position[0]); - - rho2 = rho + bw*metric->scaling[0] * rand_type(kernel); - psi2 = psi + bw*metric->scaling[1]*M_PI/180 * rand_type(kernel); - - while((rho2 < rho_min) || (rho2 > rho_max)){ - if(rho2 < rho_min) rho2 += 2 * (rho_min - rho2); - else rho2 -= 2 * (rho2 - rho_max); - } - - while((psi2 < psi_min) || (psi2>psi_max)){ - if(psi2 < psi_min) psi2 += 2 * (psi_min - psi2); - else psi2 -= 2 * (psi2 - psi_max); - } - - part->position[0] = rho2*cos(psi2); - part->position[1] = rho2*sin(psi2); - return 0; +int SurfR_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + double rho_min = metric->params[0], rho_max = metric->params[1]; + double psi_min = metric->params[2], psi_max = metric->params[3]; + double rho, psi, rho2, psi2; + + rho = sqrt(part->position[0] * part->position[0] + + part->position[1] * part->position[1]); + psi = atan2(part->position[1], part->position[0]); + + rho2 = rho + bw * metric->scaling[0] * rand_type(kernel); + psi2 = psi + bw * metric->scaling[1] * M_PI / 180 * rand_type(kernel); + + while ((rho2 < rho_min) || (rho2 > rho_max)) { + if (rho2 < rho_min) + rho2 += 2 * (rho_min - rho2); + else + rho2 -= 2 * (rho2 - rho_max); + } + + while ((psi2 < psi_min) || (psi2 > psi_max)) { + if (psi2 < psi_min) + psi2 += 2 * (psi_min - psi2); + else + psi2 -= 2 * (psi2 - psi_max); + } + + part->position[0] = rho2 * cos(psi2); + part->position[1] = rho2 * sin(psi2); + return 0; } -int SurfR2_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - double rho_min=metric->params[0], rho_max=metric->params[1]; - double psi_min=metric->params[2], psi_max=metric->params[3]; - double rho, psi, rho2, psi2; - - rho_min *= rho_min; - rho_max *= rho_max; - - rho = part->position[0]*part->position[0]+part->position[1]*part->position[1]; - psi = atan2(part->position[1], part->position[0]); - - rho2 = rho + bw*metric->scaling[0] * rand_type(kernel); - psi2 = psi + bw*metric->scaling[1]*M_PI/180 * rand_type(kernel); - - while((rho2 < rho_min) || (rho2 > rho_max)){ - if(rho2 < rho_min) rho2 += 2 * (rho_min - rho2); - else rho2 -= 2 * (rho2 - rho_max); - } - - while((psi2 < psi_min) || (psi2>psi_max)){ - if(psi2 < psi_min) psi2 += 2 * (psi_min - psi2); - else psi2 -= 2 * (psi2 - psi_max); - } - - part->position[0] = sqrt(rho2)*cos(psi2); - part->position[1] = sqrt(rho2)*sin(psi2); - return 0; +int SurfR2_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + double rho_min = metric->params[0], rho_max = metric->params[1]; + double psi_min = metric->params[2], psi_max = metric->params[3]; + double rho, psi, rho2, psi2; + + rho_min *= rho_min; + rho_max *= rho_max; + + rho = part->position[0] * part->position[0] + + part->position[1] * part->position[1]; + psi = atan2(part->position[1], part->position[0]); + + rho2 = rho + bw * metric->scaling[0] * rand_type(kernel); + psi2 = psi + bw * metric->scaling[1] * M_PI / 180 * rand_type(kernel); + + while ((rho2 < rho_min) || (rho2 > rho_max)) { + if (rho2 < rho_min) + rho2 += 2 * (rho_min - rho2); + else + rho2 -= 2 * (rho2 - rho_max); + } + + while ((psi2 < psi_min) || (psi2 > psi_max)) { + if (psi2 < psi_min) + psi2 += 2 * (psi_min - psi2); + else + psi2 -= 2 * (psi2 - psi_max); + } + + part->position[0] = sqrt(rho2) * cos(psi2); + part->position[1] = sqrt(rho2) * sin(psi2); + return 0; } -int SurfCircle_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - double rho_min=metric->params[0], rho_max=metric->params[1]; - double psi_min=metric->params[2], psi_max=metric->params[3]; - double x, y, rho2=0, psi2=0; - - x = part->position[0] + bw*metric->scaling[0] * rand_type(kernel); - y = part->position[1] + bw*metric->scaling[1] * rand_type(kernel); - - rho2 = sqrt(x * x + y * y); - psi2 = atan2(y, x); - - while((rho2 < rho_min) || (rho2 > rho_max)){ - if(rho2 < rho_min) rho2 += 2 * (rho_min - rho2); - else rho2 -= 2 * (rho2 - rho_max); - } - - while((psi2 < psi_min) || (psi2>psi_max)){ - if(psi2 < psi_min) psi2 += 2 * (psi_min - psi2); - else psi2 -= 2 * (psi2 - psi_max); - } - part->position[0] = rho2*cos(psi2); - part->position[1] = rho2*sin(psi2); - return 0; +int SurfCircle_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + double rho_min = metric->params[0], rho_max = metric->params[1]; + double psi_min = metric->params[2], psi_max = metric->params[3]; + double x, y, rho2 = 0, psi2 = 0; + + x = part->position[0] + bw * metric->scaling[0] * rand_type(kernel); + y = part->position[1] + bw * metric->scaling[1] * rand_type(kernel); + + rho2 = sqrt(x * x + y * y); + psi2 = atan2(y, x); + + while ((rho2 < rho_min) || (rho2 > rho_max)) { + if (rho2 < rho_min) + rho2 += 2 * (rho_min - rho2); + else + rho2 -= 2 * (rho2 - rho_max); + } + + while ((psi2 < psi_min) || (psi2 > psi_max)) { + if (psi2 < psi_min) + psi2 += 2 * (psi_min - psi2); + else + psi2 -= 2 * (psi2 - psi_max); + } + part->position[0] = rho2 * cos(psi2); + part->position[1] = rho2 * sin(psi2); + return 0; } -int Guide_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - double x=part->position[0], y=part->position[1], z=part->position[2]; - double dx=part->direction[0], dy=part->direction[1], dz=part->direction[2], dx2, dz2; - double xwidth=metric->params[0], yheight=metric->params[1], zmax=metric->params[2], rcurv=metric->params[3]; - double t, mu, mu2, phi; - int cont=0, mirror; - if(rcurv != 0){ // Transform to curved guide variables - double r = sqrt((rcurv+x)*(rcurv+x) + z*z); - x = copysign(1, rcurv) * r - rcurv; z = fabs(rcurv) * asin(z / r); - dx2 = dx; dz2 = dz; - dx = dx2*cos(z/rcurv) + dz2*sin(z/rcurv); dz = -dx2*sin(z/rcurv) + dz2*cos(z/rcurv); - } - // Transform from (x,y,z,dx,dy,dz) to (z,t,mu,phi) - if ((y/yheight > -x/xwidth) && (y/yheight < x/xwidth)) mirror=0; // mirror x pos - else if((y/yheight > x/xwidth) && (y/yheight > -x/xwidth)) mirror=1; // mirror y pos - else if((y/yheight < -x/xwidth) && (y/yheight > x/xwidth)) mirror=2; // mirror x neg - else mirror=3; // mirror y neg - switch(mirror){ - case 0: - t = 0.5*yheight + y; - mu = dx; phi = atan2(-dy, dz); - break; - case 1: - t = 1.0*yheight + 0.5*xwidth - x; - mu = dy; phi = atan2(dx, dz); - break; - case 2: - t = 1.5*yheight + 1.0*xwidth - y; - mu = -dx; phi = atan2(dy, dz); - break; - case 3: - t = 2.0*yheight + 1.5*xwidth + x; - mu = -dy; phi = atan2(-dx, dz); - break; - } - // Perturb - z += bw*metric->scaling[0] * rand_type(kernel); - while((z < 0) || (z > zmax)){ - if(z < 0) z *= -1; - else z -= 2 * (z - zmax); - } - t += bw*metric->scaling[1] * rand_type(kernel); - switch(mirror){ // Avoid perturbation to change mirror - case 0: - while((t < 0) || (t > yheight)){ - if(t < 0) t *= -1; - else t -= 2 * (t - yheight); - } - break; - case 1: - while((t < yheight) || (t > yheight+xwidth)){ - if(t < yheight) t += 2 * (yheight - t); - else t -= 2 * (t - (yheight+xwidth)); - } - break; - case 2: - while((t < yheight+xwidth) || (t > 2*yheight+xwidth)){ - if(t < yheight+xwidth) t += 2 * (yheight+xwidth - t); - else t -= 2 * (t - (2*yheight+xwidth)); - } - break; - case 3: - while((t < 2*yheight+xwidth) || (t > 2*yheight+2*xwidth)){ - if(t < 2*yheight+xwidth) t += 2 * (2*yheight+xwidth - t); - else t -= 2 * (t - (2*yheight+2*xwidth)); - } - break; - } - if(isinf(metric->scaling[2])) mu = -1 + 2.*rand()/RAND_MAX; - else{ - mu2 = mu + bw*metric->scaling[2] * rand_type(kernel); - if(mu >= 0){ - while((mu2 < 0) || (mu2 > 1)){ - if(mu2 < 0) mu2 *= -1; - else mu2 -= 2 * (mu2 - 1); - } - } - else{ - while((mu2 < -1) || (mu2 > 0)){ - if(mu2 < -1) mu2 += 2 * (-1 - mu2); - else mu2 *= -1; - } - } - mu = mu2; - } - if(isinf(metric->scaling[3])) phi = 2.*M_PI*rand()/RAND_MAX; - else phi += bw*metric->scaling[3]*M_PI/180 * rand_type(kernel); - // Antitransform from (z,t,theta_n,theta_t) to (x,y,z,dx,dy,dz) - switch(mirror){ - case 0: - x = xwidth/2; y = t - 0.5*yheight; - dx = mu; dz = sqrt(1-mu*mu)*cos(phi); dy = -sqrt(1-mu*mu)*sin(phi); - break; - case 1: - y = yheight/2; x = -t + yheight + 0.5*xwidth; - dy = mu; dz = sqrt(1-mu*mu)*cos(phi); dx = sqrt(1-mu*mu)*sin(phi); - break; - case 2: - x = -xwidth/2; y = -t + 1.5*yheight + xwidth; - dx = -mu; dz = sqrt(1-mu*mu)*cos(phi); dy = sqrt(1-mu*mu)*sin(phi); - break; - case 3: - y = -yheight/2; x = t - 2*yheight - 1.5*xwidth; - dy = -mu; dz = sqrt(1-mu*mu)*cos(phi); dx = -sqrt(1-mu*mu)*sin(phi); - break; - } - if(rcurv != 0){ // Antitransform curved guide variables - double r = (rcurv + x) * copysign(1, rcurv), ang = z / rcurv; - x = copysign(1, rcurv) * r * cos(ang) - rcurv; z = r * sin(fabs(ang)); - dx2 = dx; dz2 = dz; dx = dx2*cos(ang) - dz2*sin(ang); dz = dx2*sin(ang) + dz2*cos(ang); - } - part->position[0] = x; part->position[1] = y; part->position[2] = z; - part->direction[0] = dx; part->direction[1] = dy; part->direction[2] = dz; - return 0; +int Guide_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + double x = part->position[0], y = part->position[1], z = part->position[2]; + double dx = part->direction[0], dy = part->direction[1], + dz = part->direction[2], dx2, dz2; + double xwidth = metric->params[0], yheight = metric->params[1], + zmax = metric->params[2], rcurv = metric->params[3]; + double t, mu, mu2, phi; + int cont = 0, mirror; + if (rcurv != 0) { // Transform to curved guide variables + double r = sqrt((rcurv + x) * (rcurv + x) + z * z); + x = copysign(1, rcurv) * r - rcurv; + z = fabs(rcurv) * asin(z / r); + dx2 = dx; + dz2 = dz; + dx = dx2 * cos(z / rcurv) + dz2 * sin(z / rcurv); + dz = -dx2 * sin(z / rcurv) + dz2 * cos(z / rcurv); + } + // Transform from (x,y,z,dx,dy,dz) to (z,t,mu,phi) + if ((y / yheight > -x / xwidth) && (y / yheight < x / xwidth)) + mirror = 0; // mirror x pos + else if ((y / yheight > x / xwidth) && (y / yheight > -x / xwidth)) + mirror = 1; // mirror y pos + else if ((y / yheight < -x / xwidth) && (y / yheight > x / xwidth)) + mirror = 2; // mirror x neg + else + mirror = 3; // mirror y neg + switch (mirror) { + case 0: + t = 0.5 * yheight + y; + mu = dx; + phi = atan2(-dy, dz); + break; + case 1: + t = 1.0 * yheight + 0.5 * xwidth - x; + mu = dy; + phi = atan2(dx, dz); + break; + case 2: + t = 1.5 * yheight + 1.0 * xwidth - y; + mu = -dx; + phi = atan2(dy, dz); + break; + case 3: + t = 2.0 * yheight + 1.5 * xwidth + x; + mu = -dy; + phi = atan2(-dx, dz); + break; + } + // Perturb + z += bw * metric->scaling[0] * rand_type(kernel); + while ((z < 0) || (z > zmax)) { + if (z < 0) + z *= -1; + else + z -= 2 * (z - zmax); + } + t += bw * metric->scaling[1] * rand_type(kernel); + switch (mirror) { // Avoid perturbation to change mirror + case 0: + while ((t < 0) || (t > yheight)) { + if (t < 0) + t *= -1; + else + t -= 2 * (t - yheight); + } + break; + case 1: + while ((t < yheight) || (t > yheight + xwidth)) { + if (t < yheight) + t += 2 * (yheight - t); + else + t -= 2 * (t - (yheight + xwidth)); + } + break; + case 2: + while ((t < yheight + xwidth) || (t > 2 * yheight + xwidth)) { + if (t < yheight + xwidth) + t += 2 * (yheight + xwidth - t); + else + t -= 2 * (t - (2 * yheight + xwidth)); + } + break; + case 3: + while ((t < 2 * yheight + xwidth) || (t > 2 * yheight + 2 * xwidth)) { + if (t < 2 * yheight + xwidth) + t += 2 * (2 * yheight + xwidth - t); + else + t -= 2 * (t - (2 * yheight + 2 * xwidth)); + } + break; + } + if (isinf(metric->scaling[2])) + mu = -1 + 2. * rand() / RAND_MAX; + else { + mu2 = mu + bw * metric->scaling[2] * rand_type(kernel); + if (mu >= 0) { + while ((mu2 < 0) || (mu2 > 1)) { + if (mu2 < 0) + mu2 *= -1; + else + mu2 -= 2 * (mu2 - 1); + } + } else { + while ((mu2 < -1) || (mu2 > 0)) { + if (mu2 < -1) + mu2 += 2 * (-1 - mu2); + else + mu2 *= -1; + } + } + mu = mu2; + } + if (isinf(metric->scaling[3])) + phi = 2. * M_PI * rand() / RAND_MAX; + else + phi += bw * metric->scaling[3] * M_PI / 180 * rand_type(kernel); + // Antitransform from (z,t,theta_n,theta_t) to (x,y,z,dx,dy,dz) + switch (mirror) { + case 0: + x = xwidth / 2; + y = t - 0.5 * yheight; + dx = mu; + dz = sqrt(1 - mu * mu) * cos(phi); + dy = -sqrt(1 - mu * mu) * sin(phi); + break; + case 1: + y = yheight / 2; + x = -t + yheight + 0.5 * xwidth; + dy = mu; + dz = sqrt(1 - mu * mu) * cos(phi); + dx = sqrt(1 - mu * mu) * sin(phi); + break; + case 2: + x = -xwidth / 2; + y = -t + 1.5 * yheight + xwidth; + dx = -mu; + dz = sqrt(1 - mu * mu) * cos(phi); + dy = sqrt(1 - mu * mu) * sin(phi); + break; + case 3: + y = -yheight / 2; + x = t - 2 * yheight - 1.5 * xwidth; + dy = -mu; + dz = sqrt(1 - mu * mu) * cos(phi); + dx = -sqrt(1 - mu * mu) * sin(phi); + break; + } + if (rcurv != 0) { // Antitransform curved guide variables + double r = (rcurv + x) * copysign(1, rcurv), ang = z / rcurv; + x = copysign(1, rcurv) * r * cos(ang) - rcurv; + z = r * sin(fabs(ang)); + dx2 = dx; + dz2 = dz; + dx = dx2 * cos(ang) - dz2 * sin(ang); + dz = dx2 * sin(ang) + dz2 * cos(ang); + } + part->position[0] = x; + part->position[1] = y; + part->position[2] = z; + part->direction[0] = dx; + part->direction[1] = dy; + part->direction[2] = dz; + return 0; } -void _vMF_perturb(double bw, double* dx, double* dy, double* dz){ - double xi = (double)rand()/RAND_MAX; - double w = 1; - w += bw*bw * log(xi+(1-xi)*exp(-2/(bw*bw))); - double phi = 2.*M_PI*rand()/RAND_MAX; - double uv = sqrt(1-w*w), u = uv*cos(phi), v = uv*sin(phi); - double dx0=*dx, dy0=*dy, dz0=*dz; - if(dz0 > 0){ - *dx = u*dz0 + w*dx0 - (v*dx0-u*dy0)*dy0/(1+dz0); - *dy = v*dz0 + w*dy0 + (v*dx0-u*dy0)*dx0/(1+dz0); - } - else{ - *dx = u*dz0 + w*dx0 + (v*dx0-u*dy0)*dy0/(1-dz0); - *dy = v*dz0 + w*dy0 - (v*dx0-u*dy0)*dx0/(1-dz0); - } - *dz = w*dz0 - u*dx0 - v*dy0; +void _vMF_perturb(double bw, double *dx, double *dy, double *dz) { + double xi = (double)rand() / RAND_MAX; + double w = 1; + w += bw * bw * log(xi + (1 - xi) * exp(-2 / (bw * bw))); + double phi = 2. * M_PI * rand() / RAND_MAX; + double uv = sqrt(1 - w * w), u = uv * cos(phi), v = uv * sin(phi); + double dx0 = *dx, dy0 = *dy, dz0 = *dz; + if (dz0 > 0) { + *dx = u * dz0 + w * dx0 - (v * dx0 - u * dy0) * dy0 / (1 + dz0); + *dy = v * dz0 + w * dy0 + (v * dx0 - u * dy0) * dx0 / (1 + dz0); + } else { + *dx = u * dz0 + w * dx0 + (v * dx0 - u * dy0) * dy0 / (1 - dz0); + *dy = v * dz0 + w * dy0 - (v * dx0 - u * dy0) * dx0 / (1 - dz0); + } + *dz = w * dz0 - u * dx0 - v * dy0; } -int Isotrop_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - if(isinf(bw*metric->scaling[0])){ // Generate isotropic direction - part->direction[2] = -1 + 2.*rand()/RAND_MAX; - double dxy = sqrt(1-part->direction[2]*part->direction[2]); - double phi = 2.*M_PI*rand()/RAND_MAX; - part->direction[0] = dxy*cos(phi); - part->direction[1] = dxy*sin(phi); - } - else if(bw*metric->scaling[0] > 0){ // Perturb following von Mises-Fischer distribution - double dx=part->direction[0], dy=part->direction[1], dz=part->direction[2]; - _vMF_perturb(bw*metric->scaling[0], &part->direction[0], &part->direction[1], &part->direction[2]); - int keepx = (int)metric->params[0], keepy = (int)metric->params[1], keepz = (int)metric->params[2]; - if(keepx && (dx*part->direction[0] < 0)) part->direction[0] *= -1; - if(keepy && (dy*part->direction[1] < 0)) part->direction[1] *= -1; - if(keepz && (dz*part->direction[2] < 0)) part->direction[2] *= -1; - } - return 0; +int Isotrop_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + if (isinf(bw * metric->scaling[0])) { // Generate isotropic direction + part->direction[2] = -1 + 2. * rand() / RAND_MAX; + double dxy = sqrt(1 - part->direction[2] * part->direction[2]); + double phi = 2. * M_PI * rand() / RAND_MAX; + part->direction[0] = dxy * cos(phi); + part->direction[1] = dxy * sin(phi); + } else if (bw * metric->scaling[0] > + 0) { // Perturb following von Mises-Fischer distribution + double dx = part->direction[0], dy = part->direction[1], + dz = part->direction[2]; + _vMF_perturb(bw * metric->scaling[0], &part->direction[0], + &part->direction[1], &part->direction[2]); + int keepx = (int)metric->params[0], keepy = (int)metric->params[1], + keepz = (int)metric->params[2]; + if (keepx && (dx * part->direction[0] < 0)) + part->direction[0] *= -1; + if (keepy && (dy * part->direction[1] < 0)) + part->direction[1] *= -1; + if (keepz && (dz * part->direction[2] < 0)) + part->direction[2] *= -1; + } + return 0; } -int Polar_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - double theta, theta2, phi; - int cont=0; - theta = acos(part->direction[2]); - phi = atan2(part->direction[1], part->direction[0]); - theta2 = theta + bw*metric->scaling[0]*M_PI/180 * rand_type(kernel); - if(cos(theta2)*cos(theta) < 0) theta += 2 * (M_PI/2 - theta); - theta = theta2; - phi += bw*metric->scaling[1]*M_PI/180 * rand_type(kernel); - part->direction[0] = sin(theta) * cos(phi); - part->direction[1] = sin(theta) * sin(phi); - part->direction[2] = cos(theta); - return 0; +int Polar_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + double theta, theta2, phi; + int cont = 0; + theta = acos(part->direction[2]); + phi = atan2(part->direction[1], part->direction[0]); + theta2 = theta + bw * metric->scaling[0] * M_PI / 180 * rand_type(kernel); + if (cos(theta2) * cos(theta) < 0) + theta += 2 * (M_PI / 2 - theta); + theta = theta2; + phi += bw * metric->scaling[1] * M_PI / 180 * rand_type(kernel); + part->direction[0] = sin(theta) * cos(phi); + part->direction[1] = sin(theta) * sin(phi); + part->direction[2] = cos(theta); + return 0; } -int PolarMu_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel){ - double mu, mu2, phi; - int cont=0; - mu = part->direction[2]; - phi = atan2(part->direction[1], part->direction[0]); - mu2 = mu + bw*metric->scaling[0] * rand_type(kernel); - if(mu >= 0){ - while((mu2 < 0) || (mu2 > 1)){ - if(mu2 < 0) mu2 *= -1; - if(mu2 > 1) mu2 -= 2 * (mu2 - 1); - } - } - else{ - while((mu2 < -1) || (mu2 > 0)){ - if(mu2 < -1) mu2 += 2 * (-1 - mu2); - if(mu2 > 0) mu2 *= -1; - } - } - mu = mu2; - phi += bw*metric->scaling[1]*M_PI/180 * rand_type(kernel); - part->direction[0] = sqrt(1-mu*mu) * cos(phi); - part->direction[1] = sqrt(1-mu*mu) * sin(phi); - part->direction[2] = mu; - return 0; +int PolarMu_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel) { + double mu, mu2, phi; + int cont = 0; + mu = part->direction[2]; + phi = atan2(part->direction[1], part->direction[0]); + mu2 = mu + bw * metric->scaling[0] * rand_type(kernel); + if (mu >= 0) { + while ((mu2 < 0) || (mu2 > 1)) { + if (mu2 < 0) + mu2 *= -1; + if (mu2 > 1) + mu2 -= 2 * (mu2 - 1); + } + } else { + while ((mu2 < -1) || (mu2 > 0)) { + if (mu2 < -1) + mu2 += 2 * (-1 - mu2); + if (mu2 > 0) + mu2 *= -1; + } + } + mu = mu2; + phi += bw * metric->scaling[1] * M_PI / 180 * rand_type(kernel); + part->direction[0] = sqrt(1 - mu * mu) * cos(phi); + part->direction[1] = sqrt(1 - mu * mu) * sin(phi); + part->direction[2] = mu; + return 0; } diff --git a/src/kdsource/geom.h b/src/kdsource/geom.h index 56564cd..7c6cdd7 100644 --- a/src/kdsource/geom.h +++ b/src/kdsource/geom.h @@ -1,77 +1,100 @@ #ifndef METRICS_H #define METRICS_H -#include -#include +#include +#include /***********************************************************************************/ /* */ -/* Utilities for geometry and metrics handling for KDSource. */ +/* Utilities for geometry and metrics handling for KDSource. */ /* */ -/* This file can be freely used as per the terms in the LICENSE file. */ +/* This file can be freely used as per the terms in the LICENSE file. */ /* */ -/* Written by Osiris Inti Abbate, 2021. */ +/* Written by Osiris Inti Abbate, 2021. */ /* */ /***********************************************************************************/ #include "mcpl.h" - typedef struct Metric Metric; -typedef int (*PerturbFun)(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); +typedef int (*PerturbFun)(const Metric *metric, mcpl_particle_t *part, + double bw, char kernel); -struct Metric{ - int dim; // Dimension - float* scaling; // Variables scaling - PerturbFun perturb; // Perturbation function - int nps; // Number of metric parameters - double* params; // Metric parameters +struct Metric { + int dim; // Dimension + float *scaling; // Variables scaling + PerturbFun perturb; // Perturbation function + int nps; // Number of metric parameters + double *params; // Metric parameters }; -Metric* Metric_create(int dim, const double* scaling, PerturbFun perturb, int nps, const double* params); -Metric* Metric_copy(const Metric* from); -void Metric_destroy(Metric* metric); - -typedef struct Geometry{ - int ord; // Number of submetrics - Metric** ms; // Submetrics - char* bwfilename; // Bandwidth file name - FILE* bwfile; // Bandwidth file - double bw; // Normalized bandwidth - char kernel; // Kernel - - double* trasl; // Geometry translation - double* rot; // Geometry rotation +Metric *Metric_create(int dim, const double *scaling, PerturbFun perturb, + int nps, const double *params); +Metric *Metric_copy(const Metric *from); +void Metric_destroy(Metric *metric); + +typedef struct Geometry { + int ord; // Number of submetrics + Metric **ms; // Submetrics + char *bwfilename; // Bandwidth file name + FILE *bwfile; // Bandwidth file + double bw; // Normalized bandwidth + char kernel; // Kernel + + double *trasl; // Geometry translation + double *rot; // Geometry rotation } Geometry; -Geometry* Geom_create(int ord, Metric** metrics, double bw, const char* bwfilename, char kernel, - const double* trasl, const double* rot); -Geometry* Geom_copy(const Geometry* from); -int Geom_perturb(const Geometry* geom, mcpl_particle_t* part); -int Geom_next(Geometry* geom, int loop); -void Geom_destroy(Geometry* geom); - -int E_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); -int Let_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); -int wl_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); - -int t_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); -int Dec_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); - -int Vol_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); -int SurfXY_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); -int SurfR_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); -int SurfR2_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); -int SurfCircle_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); -int Guide_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); - -int Isotrop_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); -int Polar_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); -int PolarMu_perturb(const Metric* metric, mcpl_particle_t* part, double bw, char kernel); +Geometry *Geom_create(int ord, Metric **metrics, double bw, + const char *bwfilename, char kernel, const double *trasl, + const double *rot); +Geometry *Geom_copy(const Geometry *from); +int Geom_perturb(const Geometry *geom, mcpl_particle_t *part); +int Geom_next(Geometry *geom, int loop); +void Geom_destroy(Geometry *geom); + +int E_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); +int Let_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); +int wl_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); + +int t_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); +int Dec_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); + +int Vol_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); +int SurfXY_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); +int SurfR_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); +int SurfR2_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); +int SurfCircle_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); +int Guide_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); + +int Isotrop_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); +int Polar_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); +int PolarMu_perturb(const Metric *metric, mcpl_particle_t *part, double bw, + char kernel); static const int _n_metrics = 14; -static const char *_metric_names[] = {"Energy", "Lethargy", "Wavelength", "Vol", "SurfXY", "SurfR", "SurfR2", "SurfCircle", "Guide", "Isotrop", "Polar", "PolarMu", "Time", "Decade"}; -static const PerturbFun _metric_perturbs[] = {E_perturb, Let_perturb, wl_perturb, Vol_perturb, SurfXY_perturb, SurfR_perturb, SurfR2_perturb, SurfCircle_perturb, Guide_perturb, Isotrop_perturb, Polar_perturb, PolarMu_perturb, t_perturb, Dec_perturb}; +static const char *_metric_names[] = { + "Energy", "Lethargy", "Wavelength", "Vol", "SurfXY", + "SurfR", "SurfR2", "SurfCircle", "Guide", "Isotrop", + "Polar", "PolarMu", "Time", "Decade"}; +static const PerturbFun _metric_perturbs[] = { + E_perturb, Let_perturb, wl_perturb, Vol_perturb, + SurfXY_perturb, SurfR_perturb, SurfR2_perturb, SurfCircle_perturb, + Guide_perturb, Isotrop_perturb, Polar_perturb, PolarMu_perturb, + t_perturb, Dec_perturb}; #endif diff --git a/src/kdsource/kdsource.c b/src/kdsource/kdsource.c index ab2517c..a0add5c 100644 --- a/src/kdsource/kdsource.c +++ b/src/kdsource/kdsource.c @@ -1,302 +1,336 @@ -#include -#include -#include +#include +#include +#include -#include -#include +#include +#include #include "kdsource.h" - -void KDS_error(const char* msg){ - printf("KDSource error: %s\n", msg); - exit(EXIT_FAILURE); +void KDS_error(const char *msg) { + printf("KDSource error: %s\n", msg); + exit(EXIT_FAILURE); } -void KDS_end(const char* msg){ - printf("KDSource terminate: %s\n", msg); - exit(EXIT_SUCCESS); +void KDS_end(const char *msg) { + printf("KDSource terminate: %s\n", msg); + exit(EXIT_SUCCESS); } -KDSource* KDS_create(double J, char kernel, PList* plist, Geometry* geom){ - KDSource* kds = (KDSource*)malloc(sizeof(KDSource)); - kds->J = J; - kds->kernel = kernel; - kds->plist = plist; - kds->geom = geom; - return kds; +KDSource *KDS_create(double J, char kernel, PList *plist, Geometry *geom) { + KDSource *kds = (KDSource *)malloc(sizeof(KDSource)); + kds->J = J; + kds->kernel = kernel; + kds->plist = plist; + kds->geom = geom; + return kds; } -KDSource* KDS_open(const char* xmlfilename){ - xmlKeepBlanksDefault(0); +KDSource *KDS_open(const char *xmlfilename) { + xmlKeepBlanksDefault(0); - int i, j, n; - double J; - char kernel='g'; - char *buf; + int i, j, n; + double J; + char kernel = 'g'; + char *buf; - char pt; - char mcplfile[NAME_MAX_LEN]; - double *trasl_plist=NULL, *rot_plist=NULL; + char pt; + char mcplfile[NAME_MAX_LEN]; + double *trasl_plist = NULL, *rot_plist = NULL; - int order; - double *trasl_geom=NULL, *rot_geom=NULL; - int switch_x2z, variable_bw; - char* bwfilename=NULL; - double bw=0; + int order; + double *trasl_geom = NULL, *rot_geom = NULL; + int switch_x2z, variable_bw; + char *bwfilename = NULL; + double bw = 0; - // Read file - printf("Reading xmlfile %s...\n", xmlfilename); - xmlDocPtr doc = xmlReadFile(xmlfilename, NULL, 0); - if(doc == NULL){ - printf("Could not open file %s\n", xmlfilename); - KDS_error("Error in KDS_open"); - } - xmlNodePtr root = xmlDocGetRootElement(doc); - if(strcmp((char*)root->name, "KDSource") != 0){ - printf("Invalid format in source XML file %s\n", xmlfilename); - KDS_error("Error in KDS_open"); - } - xmlNodePtr node = root->children; // Node: J - sscanf((char*)xmlNodeGetContent(node), "%lf", &J); // Read J + // Read file + printf("Reading xmlfile %s...\n", xmlfilename); + xmlDocPtr doc = xmlReadFile(xmlfilename, NULL, 0); + if (doc == NULL) { + printf("Could not open file %s\n", xmlfilename); + KDS_error("Error in KDS_open"); + } + xmlNodePtr root = xmlDocGetRootElement(doc); + if (strcmp((char *)root->name, "KDSource") != 0) { + printf("Invalid format in source XML file %s\n", xmlfilename); + KDS_error("Error in KDS_open"); + } + xmlNodePtr node = root->children; // Node: J + sscanf((char *)xmlNodeGetContent(node), "%lf", &J); // Read J - node = node->next; // Node: kernel - if(strcmp((char*)node->name, "kernel") == 0){ - sscanf((char*)xmlNodeGetContent(node), "%s", &kernel); // Read kernel - node = node->next; - } - else printf("No kernel specified. Using gaussian as default.\n"); + node = node->next; // Node: kernel + if (strcmp((char *)node->name, "kernel") == 0) { + sscanf((char *)xmlNodeGetContent(node), "%s", &kernel); // Read kernel + node = node->next; + } else + printf("No kernel specified. Using gaussian as default.\n"); - // PList - xmlNodePtr pltree = node; - node = pltree->children; // Node: pt - sscanf((char*)xmlNodeGetContent(node), "%c", &pt); // Read pt - node = node->next; // Node: mcplname - if(strlen((char*)xmlNodeGetContent(node)) > NAME_MAX_LEN){ - printf("mcpl file name %s exceeds NAME_MAX_LEN=%d", (char*)xmlNodeGetContent(node), NAME_MAX_LEN); - KDS_error("Error in KDS_open"); - } - strcpy(mcplfile, (char*)xmlNodeGetContent(node)); // Read mcpl file name - node = node->next; // Node: trasl - if(strlen((char*)xmlNodeGetContent(node)) > 1){ - trasl_plist = (double*)malloc(3 * sizeof(double)); - sscanf((char*)xmlNodeGetContent(node), "%lf %lf %lf", &trasl_plist[0], &trasl_plist[1], &trasl_plist[2]); - } - node = node->next; // Node: rot - if(strlen((char*)xmlNodeGetContent(node)) > 1){ - rot_plist = (double*)malloc(3 * sizeof(double)); - sscanf((char*)xmlNodeGetContent(node), "%lf %lf %lf", &rot_plist[0], &rot_plist[1], &rot_plist[2]); - } - node = node->next; // Node: x2z - sscanf((char*)xmlNodeGetContent(node), "%d", &switch_x2z); // Read switch_x2z + // PList + xmlNodePtr pltree = node; + node = pltree->children; // Node: pt + sscanf((char *)xmlNodeGetContent(node), "%c", &pt); // Read pt + node = node->next; // Node: mcplname + if (strlen((char *)xmlNodeGetContent(node)) > NAME_MAX_LEN) { + printf("mcpl file name %s exceeds NAME_MAX_LEN=%d", + (char *)xmlNodeGetContent(node), NAME_MAX_LEN); + KDS_error("Error in KDS_open"); + } + strcpy(mcplfile, (char *)xmlNodeGetContent(node)); // Read mcpl file name + node = node->next; // Node: trasl + if (strlen((char *)xmlNodeGetContent(node)) > 1) { + trasl_plist = (double *)malloc(3 * sizeof(double)); + sscanf((char *)xmlNodeGetContent(node), "%lf %lf %lf", &trasl_plist[0], + &trasl_plist[1], &trasl_plist[2]); + } + node = node->next; // Node: rot + if (strlen((char *)xmlNodeGetContent(node)) > 1) { + rot_plist = (double *)malloc(3 * sizeof(double)); + sscanf((char *)xmlNodeGetContent(node), "%lf %lf %lf", &rot_plist[0], + &rot_plist[1], &rot_plist[2]); + } + node = node->next; // Node: x2z + sscanf((char *)xmlNodeGetContent(node), "%d", &switch_x2z); // Read switch_x2z - // Geometry - xmlNodePtr gtree = pltree->next; - sscanf((char*)xmlGetProp(gtree, (const xmlChar*)"order"), "%d", &order); // Read order - int dims[order], nps[order]; - char metricnames[order][NAME_MAX_LEN]; - double *params[order]; - double *scalings[order]; - PerturbFun perturbs[order]; - xmlNodePtr mtree = gtree->children; - for(i=0; iname); // Read metricname - node = mtree->children; // Node: dim - sscanf((char*)xmlNodeGetContent(node), "%d", &dims[i]); // Read dim - scalings[i] = (double*)malloc(dims[i]*sizeof(double)); - node = node->next; // Node: params - sscanf((char*)xmlGetProp(node,(const xmlChar*)"nps"), "%d", &nps[i]); // Read ngp - params[i] = (double*)malloc(nps[i]*sizeof(double)); - buf = (char*)xmlNodeGetContent(node); - for(j=0; jnext; - } - node = mtree; // Node: trasl - if(strlen((char*)xmlNodeGetContent(node)) > 1){ - trasl_geom = (double*)malloc(3 * sizeof(double)); - sscanf((char*)xmlNodeGetContent(node), "%lf %lf %lf", &trasl_geom[0], &trasl_geom[1], &trasl_geom[2]); - } - node = node->next; // Node: rot - if(strlen((char*)xmlNodeGetContent(node)) > 1){ - rot_geom = (double*)malloc(3 * sizeof(double)); - sscanf((char*)xmlNodeGetContent(node), "%lf %lf %lf", &rot_geom[0], &rot_geom[1], &rot_geom[2]); - } - node = gtree->next; // Node: scaling - buf = (char*)xmlNodeGetContent(node); - for(i=0; inext; // Node: BW - sscanf((char*)xmlGetProp(node,(const xmlChar*)"variable"), "%d", &variable_bw); // Read variable_bw - if(variable_bw){ - bwfilename = (char*)malloc(NAME_MAX_LEN*sizeof(char)); - if(strlen((char*)xmlNodeGetContent(node)) > NAME_MAX_LEN){ - printf("BW file name %s exceeds NAME_MAX_LEN=%d", (char*)xmlNodeGetContent(node), NAME_MAX_LEN); - KDS_error("Error in KDS_open"); - } - strcpy(bwfilename, (char*)xmlNodeGetContent(node)); // Read BW file name - } - else - sscanf((char*)xmlNodeGetContent(node), "%lf", &bw); // Read BW + // Geometry + xmlNodePtr gtree = pltree->next; + sscanf((char *)xmlGetProp(gtree, (const xmlChar *)"order"), "%d", + &order); // Read order + int dims[order], nps[order]; + char metricnames[order][NAME_MAX_LEN]; + double *params[order]; + double *scalings[order]; + PerturbFun perturbs[order]; + xmlNodePtr mtree = gtree->children; + for (i = 0; i < order; i++) { + strcpy(metricnames[i], (char *)mtree->name); // Read metricname + node = mtree->children; // Node: dim + sscanf((char *)xmlNodeGetContent(node), "%d", &dims[i]); // Read dim + scalings[i] = (double *)malloc(dims[i] * sizeof(double)); + node = node->next; // Node: params + sscanf((char *)xmlGetProp(node, (const xmlChar *)"nps"), "%d", + &nps[i]); // Read ngp + params[i] = (double *)malloc(nps[i] * sizeof(double)); + buf = (char *)xmlNodeGetContent(node); + for (j = 0; j < nps[i]; j++) { // Read params + sscanf(buf, "%lf %n", ¶ms[i][j], &n); + buf += n; + } + mtree = mtree->next; + } + node = mtree; // Node: trasl + if (strlen((char *)xmlNodeGetContent(node)) > 1) { + trasl_geom = (double *)malloc(3 * sizeof(double)); + sscanf((char *)xmlNodeGetContent(node), "%lf %lf %lf", &trasl_geom[0], + &trasl_geom[1], &trasl_geom[2]); + } + node = node->next; // Node: rot + if (strlen((char *)xmlNodeGetContent(node)) > 1) { + rot_geom = (double *)malloc(3 * sizeof(double)); + sscanf((char *)xmlNodeGetContent(node), "%lf %lf %lf", &rot_geom[0], + &rot_geom[1], &rot_geom[2]); + } + node = gtree->next; // Node: scaling + buf = (char *)xmlNodeGetContent(node); + for (i = 0; i < order; i++) { // Read scalings + for (j = 0; j < dims[i]; j++) { + sscanf(buf, "%lf %n", &scalings[i][j], &n); + buf += n; + } + } + node = node->next; // Node: BW + sscanf((char *)xmlGetProp(node, (const xmlChar *)"variable"), "%d", + &variable_bw); // Read variable_bw + if (variable_bw) { + bwfilename = (char *)malloc(NAME_MAX_LEN * sizeof(char)); + if (strlen((char *)xmlNodeGetContent(node)) > NAME_MAX_LEN) { + printf("BW file name %s exceeds NAME_MAX_LEN=%d", + (char *)xmlNodeGetContent(node), NAME_MAX_LEN); + KDS_error("Error in KDS_open"); + } + strcpy(bwfilename, (char *)xmlNodeGetContent(node)); // Read BW file name + } else + sscanf((char *)xmlNodeGetContent(node), "%lf", &bw); // Read BW - // Create PList - PList* plist = PList_create(pt, mcplfile, trasl_plist, rot_plist, switch_x2z); - // Create Metric - for(i=0; iplist, part); - if(perturb) Geom_perturb(kds->geom, part); - ret += PList_next(kds->plist, loop); - ret += Geom_next(kds->geom, loop); - } - else{ // Normalize w to 1 - double bs; - int resamples = 0; - while(1){ - PList_get(kds->plist, part); - if(bias) bs = bias(part); - else bs = 1; - if(part->weight*bs > w_crit){ // If w*bs>w_crit, use w_crit/w*bs as prob of going forward in list - if(perturb) Geom_perturb(kds->geom, part); - if(rand() < w_crit/(part->weight*bs)*RAND_MAX){ - ret += PList_next(kds->plist, loop); - ret += Geom_next(kds->geom, loop); - } - break; - } - else{ // If w*bsweight*bs)/w_crit*RAND_MAX){ - take = 1; - if(perturb) Geom_perturb(kds->geom, part); - } - ret += PList_next(kds->plist, loop); - ret += Geom_next(kds->geom, loop); - if(take) break; - } - if(resamples++ > MAX_RESAMPLES) - KDS_error("Error in KDS_sample: MAX_RESAMPLES reached."); - } - part->weight = 1/bs; - } - if(ret==1 && kds->geom->bwfile) - printf("Warning: Particle list and bandwidths file have different size.\n"); - return ret; +int KDS_sample2(KDSource *kds, mcpl_particle_t *part, int perturb, + double w_crit, WeightFun bias, int loop) { + int ret = 0; + if (w_crit <= 0) { + PList_get(kds->plist, part); + if (perturb) + Geom_perturb(kds->geom, part); + ret += PList_next(kds->plist, loop); + ret += Geom_next(kds->geom, loop); + } else { // Normalize w to 1 + double bs; + int resamples = 0; + while (1) { + PList_get(kds->plist, part); + if (bias) + bs = bias(part); + else + bs = 1; + if (part->weight * bs > w_crit) { // If w*bs>w_crit, use w_crit/w*bs as + // prob of going forward in list + if (perturb) + Geom_perturb(kds->geom, part); + if (rand() < w_crit / (part->weight * bs) * RAND_MAX) { + ret += PList_next(kds->plist, loop); + ret += Geom_next(kds->geom, loop); + } + break; + } else { // If w*bsweight * bs) / w_crit * RAND_MAX) { + take = 1; + if (perturb) + Geom_perturb(kds->geom, part); + } + ret += PList_next(kds->plist, loop); + ret += Geom_next(kds->geom, loop); + if (take) + break; + } + if (resamples++ > MAX_RESAMPLES) + KDS_error("Error in KDS_sample: MAX_RESAMPLES reached."); + } + part->weight = 1 / bs; + } + if (ret == 1 && kds->geom->bwfile) + printf("Warning: Particle list and bandwidths file have different size.\n"); + return ret; } -int KDS_sample(KDSource* kds, mcpl_particle_t* part){ - return KDS_sample2(kds, part, 1, 1, NULL, 1); +int KDS_sample(KDSource *kds, mcpl_particle_t *part) { + return KDS_sample2(kds, part, 1, 1, NULL, 1); } -double KDS_w_mean(KDSource* kds, int N, WeightFun bias){ - int i; - char pt; - mcpl_particle_t part; - double w_mean=0;; - for(i=0; iplist); - Geom_destroy(kds->geom); - free(kds); +void KDS_destroy(KDSource *kds) { + PList_destroy(kds->plist); + Geom_destroy(kds->geom); + free(kds); } - -MultiSource* MS_create(int len, KDSource** s, const double* ws){ - MultiSource* ms = (MultiSource*)malloc(sizeof(MultiSource)); - ms->len = len; - ms->s = (KDSource**)malloc(len*sizeof(KDSource*)); - ms->ws = (double*)malloc(len*sizeof(double)); - ms->J = 0; - int i; - for(i=0; is[i] = s[i]; - ms->ws[i] = ws[i]; - ms->J += s[i]->J; - } - ms->cdf = (double*)malloc(ms->len*sizeof(double)); - for(i=0; ilen; i++) ms->cdf[i] = ms->ws[i]; - for(i=1; ilen; i++) ms->cdf[i] += ms->cdf[i-1]; - return ms; +MultiSource *MS_create(int len, KDSource **s, const double *ws) { + MultiSource *ms = (MultiSource *)malloc(sizeof(MultiSource)); + ms->len = len; + ms->s = (KDSource **)malloc(len * sizeof(KDSource *)); + ms->ws = (double *)malloc(len * sizeof(double)); + ms->J = 0; + int i; + for (i = 0; i < len; i++) { + ms->s[i] = s[i]; + ms->ws[i] = ws[i]; + ms->J += s[i]->J; + } + ms->cdf = (double *)malloc(ms->len * sizeof(double)); + for (i = 0; i < ms->len; i++) + ms->cdf[i] = ms->ws[i]; + for (i = 1; i < ms->len; i++) + ms->cdf[i] += ms->cdf[i - 1]; + return ms; } -MultiSource* MS_open(int len, const char** xmlfilenames, const double* ws){ - KDSource* s[len]; - int i; - for(i=0; icdf[ms->len-1] <= 0) i = (int)(y*ms->len); - else for(i=0; y*ms->cdf[ms->len-1]>ms->cdf[i]; i++); - ret = KDS_sample2(ms->s[i], part, perturb, w_crit, bias, loop); - if(ms->cdf[ms->len-1] > 0) part->weight *= (ms->s[i]->J/ms->J) / (ms->ws[i]/ms->cdf[ms->len-1]); - else part->weight *= (ms->s[i]->J/ms->J) * ms->len; - return ret; +int MS_sample2(MultiSource *ms, mcpl_particle_t *part, int perturb, + double w_crit, WeightFun bias, int loop) { + double y = rand() / ((double)RAND_MAX + 1); + int i, ret; + if (ms->cdf[ms->len - 1] <= 0) + i = (int)(y * ms->len); + else + for (i = 0; y * ms->cdf[ms->len - 1] > ms->cdf[i]; i++) + ; + ret = KDS_sample2(ms->s[i], part, perturb, w_crit, bias, loop); + if (ms->cdf[ms->len - 1] > 0) + part->weight *= (ms->s[i]->J / ms->J) / (ms->ws[i] / ms->cdf[ms->len - 1]); + else + part->weight *= (ms->s[i]->J / ms->J) * ms->len; + return ret; } -int MS_sample(MultiSource* ms, mcpl_particle_t* part){ - return MS_sample2(ms, part, 1, 1, NULL, 1); +int MS_sample(MultiSource *ms, mcpl_particle_t *part) { + return MS_sample2(ms, part, 1, 1, NULL, 1); } -double MS_w_mean(MultiSource* ms, int N, WeightFun bias){ - double w_mean=0; - int i; - for(i=0; ilen; i++) w_mean += ms->ws[i] * KDS_w_mean(ms->s[i], N, bias); - return w_mean / ms->cdf[ms->len-1]; +double MS_w_mean(MultiSource *ms, int N, WeightFun bias) { + double w_mean = 0; + int i; + for (i = 0; i < ms->len; i++) + w_mean += ms->ws[i] * KDS_w_mean(ms->s[i], N, bias); + return w_mean / ms->cdf[ms->len - 1]; } -void MS_destroy(MultiSource* ms){ - int i; - for(i=0; ilen; i++) KDS_destroy(ms->s[i]); - free(ms->s); - free(ms->ws); - free(ms->cdf); - free(ms); +void MS_destroy(MultiSource *ms) { + int i; + for (i = 0; i < ms->len; i++) + KDS_destroy(ms->s[i]); + free(ms->s); + free(ms->ws); + free(ms->cdf); + free(ms); } diff --git a/src/kdsource/kdsource.h b/src/kdsource/kdsource.h index e314fe0..0b5553f 100644 --- a/src/kdsource/kdsource.h +++ b/src/kdsource/kdsource.h @@ -1,68 +1,69 @@ #ifndef KDSOURCE_H #define KDSOURCE_H -#include -#include +#include +#include /***********************************************************************************/ /* */ -/* KDSource: KDE particle sources */ +/* KDSource: KDE particle sources */ /* */ -/* Utilities for sampling particles from a KDE source. KDSource sources use */ -/* particle lists in MCPL format, and apply on them the Kernel Density Estimation */ -/* (KDE) method. This allows coupling different Monte Carlo radiation transport */ -/* simulations, and gives variance reduction. */ +/* Utilities for sampling particles from a KDE source. KDSource sources use */ +/* particle lists in MCPL format, and apply on them the Kernel Density + * Estimation */ +/* (KDE) method. This allows coupling different Monte Carlo radiation transport + */ +/* simulations, and gives variance reduction. */ /* */ -/* Find more information and updates at https://github.com/KDSource/KDSource */ +/* Find more information and updates at https://github.com/KDSource/KDSource */ /* */ -/* This file can be freely used as per the terms in the LICENSE file. */ +/* This file can be freely used as per the terms in the LICENSE file. */ /* */ -/* Written by Osiris Inti Abbate, 2021. */ +/* Written by Osiris Inti Abbate, 2021. */ /* */ /***********************************************************************************/ #include "mcpl.h" #include "KDSourceConfig.h" -#include "utils.h" #include "geom.h" #include "plist.h" - +#include "utils.h" #define MAX_RESAMPLES 1000000 #define NAME_MAX_LEN 256 +typedef double (*WeightFun)(const mcpl_particle_t *part); -typedef double (*WeightFun)(const mcpl_particle_t* part); - -typedef struct KDSource{ - double J; // Total current [1/s] - char kernel; // Kernel - PList* plist; // Particle list - Geometry* geom; // Geometry defining variable treatment +typedef struct KDSource { + double J; // Total current [1/s] + char kernel; // Kernel + PList *plist; // Particle list + Geometry *geom; // Geometry defining variable treatment } KDSource; -KDSource* KDS_create(double J, char kernel, PList* plist, Geometry* geom); -KDSource* KDS_open(const char* xmlfilename); -int KDS_sample2(KDSource* kds, mcpl_particle_t* part, int perturb, double w_crit, WeightFun bias, int loop); -int KDS_sample(KDSource* kds, mcpl_particle_t* part); -double KDS_w_mean(KDSource* kds, int N, WeightFun bias); -void KDS_destroy(KDSource* kds); +KDSource *KDS_create(double J, char kernel, PList *plist, Geometry *geom); +KDSource *KDS_open(const char *xmlfilename); +int KDS_sample2(KDSource *kds, mcpl_particle_t *part, int perturb, + double w_crit, WeightFun bias, int loop); +int KDS_sample(KDSource *kds, mcpl_particle_t *part); +double KDS_w_mean(KDSource *kds, int N, WeightFun bias); +void KDS_destroy(KDSource *kds); -typedef struct MultiSource{ - int len; // Number of sources - KDSource** s; // Array of sources - double J; // Total current [1/s] - double* ws; // Frequency weights of sources - double* cdf; // cdf of sources weights +typedef struct MultiSource { + int len; // Number of sources + KDSource **s; // Array of sources + double J; // Total current [1/s] + double *ws; // Frequency weights of sources + double *cdf; // cdf of sources weights } MultiSource; -MultiSource* MS_create(int len, KDSource** s, const double* ws); -MultiSource* MS_open(int len, const char** xmlfilenames, const double* ws); -int MS_sample2(MultiSource* ms, mcpl_particle_t* part, int perturb, double w_crit, WeightFun bias, int loop); -int MS_sample(MultiSource* ms, mcpl_particle_t* part); -double MS_w_mean(MultiSource* ms, int N, WeightFun bias); -void MS_destroy(MultiSource* ms); - +MultiSource *MS_create(int len, KDSource **s, const double *ws); +MultiSource *MS_open(int len, const char **xmlfilenames, const double *ws); +int MS_sample2(MultiSource *ms, mcpl_particle_t *part, int perturb, + double w_crit, WeightFun bias, int loop); +int MS_sample(MultiSource *ms, mcpl_particle_t *part); +double MS_w_mean(MultiSource *ms, int N, WeightFun bias); +void MS_destroy(MultiSource *ms); #endif diff --git a/src/kdsource/plist.c b/src/kdsource/plist.c index 55154cb..cbbdffb 100644 --- a/src/kdsource/plist.c +++ b/src/kdsource/plist.c @@ -1,100 +1,111 @@ -#include -#include -#include +#include +#include +#include -#include +#include #include "kdsource.h" +void KDS_error(const char *msg); +void KDS_end(const char *msg); -void KDS_error(const char* msg); -void KDS_end(const char* msg); - -PList* PList_create(char pt, const char* filename, const double* trasl, const double* rot, int switch_x2z){ - PList* plist = (PList*)malloc(sizeof(PList)); - plist->pt = pt; - plist->pdgcode = pt2pdg(pt); - mcpl_file_t file = mcpl_open_file(filename); - plist->npts = mcpl_hdr_nparticles(file); - plist->filename = (char*)malloc(NAME_MAX_LEN*sizeof(char)); - strcpy(plist->filename, filename); - plist->file = file; - int i; - if(trasl){ - plist->trasl = (double*)malloc(3 * sizeof(double)); - for(i=0; i<3; i++) plist->trasl[i] = trasl[i]; - } - else plist->trasl = NULL; - if(rot){ - plist->rot = (double*)malloc(3 * sizeof(double)); - for(i=0; i<3; i++) plist->rot[i] = rot[i]; - } - else plist->rot = NULL; - plist->x2z = switch_x2z; - plist->part = NULL; - PList_next(plist, 1); - return plist; +PList *PList_create(char pt, const char *filename, const double *trasl, + const double *rot, int switch_x2z) { + PList *plist = (PList *)malloc(sizeof(PList)); + plist->pt = pt; + plist->pdgcode = pt2pdg(pt); + mcpl_file_t file = mcpl_open_file(filename); + plist->npts = mcpl_hdr_nparticles(file); + plist->filename = (char *)malloc(NAME_MAX_LEN * sizeof(char)); + strcpy(plist->filename, filename); + plist->file = file; + int i; + if (trasl) { + plist->trasl = (double *)malloc(3 * sizeof(double)); + for (i = 0; i < 3; i++) + plist->trasl[i] = trasl[i]; + } else + plist->trasl = NULL; + if (rot) { + plist->rot = (double *)malloc(3 * sizeof(double)); + for (i = 0; i < 3; i++) + plist->rot[i] = rot[i]; + } else + plist->rot = NULL; + plist->x2z = switch_x2z; + plist->part = NULL; + PList_next(plist, 1); + return plist; } -PList* PList_copy(const PList* from){ - PList* plist = (PList*)malloc(sizeof(PList)); - *plist = *from; - plist->filename = (char*)malloc(NAME_MAX_LEN*sizeof(char)); - strcpy(plist->filename, from->filename); - plist->file = mcpl_open_file(plist->filename); - int i; - if(from->trasl){ - plist->trasl = (double*)malloc(3 * sizeof(double)); - for(i=0; i<3; i++) plist->trasl[i] = from->trasl[i]; - } - if(from->rot){ - plist->rot = (double*)malloc(3 * sizeof(double)); - for(i=0; i<3; i++) plist->rot[i] = from->rot[i]; - } - plist->part = NULL; - PList_next(plist, 1); - return plist; +PList *PList_copy(const PList *from) { + PList *plist = (PList *)malloc(sizeof(PList)); + *plist = *from; + plist->filename = (char *)malloc(NAME_MAX_LEN * sizeof(char)); + strcpy(plist->filename, from->filename); + plist->file = mcpl_open_file(plist->filename); + int i; + if (from->trasl) { + plist->trasl = (double *)malloc(3 * sizeof(double)); + for (i = 0; i < 3; i++) + plist->trasl[i] = from->trasl[i]; + } + if (from->rot) { + plist->rot = (double *)malloc(3 * sizeof(double)); + for (i = 0; i < 3; i++) + plist->rot[i] = from->rot[i]; + } + plist->part = NULL; + PList_next(plist, 1); + return plist; } -int PList_get(const PList* plist, mcpl_particle_t* part){ - *part = *plist->part; - if(plist->trasl) - traslv(part->position, plist->trasl, 0); - if(plist->rot){ - rotv(part->position, plist->rot, 0); - rotv(part->direction, plist->rot, 0); - } - if(plist->x2z){ - double x=part->position[0], y=part->position[1], z=part->position[2]; - part->position[0] = y; part->position[1] = z; part->position[2] = x; - double dx=part->direction[0], dy=part->direction[1], dz=part->direction[2]; - part->direction[0] = dy; part->direction[1] = dz; part->direction[2] = dx; - } - return 0; +int PList_get(const PList *plist, mcpl_particle_t *part) { + *part = *plist->part; + if (plist->trasl) + traslv(part->position, plist->trasl, 0); + if (plist->rot) { + rotv(part->position, plist->rot, 0); + rotv(part->direction, plist->rot, 0); + } + if (plist->x2z) { + double x = part->position[0], y = part->position[1], z = part->position[2]; + part->position[0] = y; + part->position[1] = z; + part->position[2] = x; + double dx = part->direction[0], dy = part->direction[1], + dz = part->direction[2]; + part->direction[0] = dy; + part->direction[1] = dz; + part->direction[2] = dx; + } + return 0; } -int PList_next(PList* plist, int loop){ - int ret=0, resamples=0; - while(resamples++ < MAX_RESAMPLES){ - plist->part = mcpl_read(plist->file); - if(plist->part == NULL){ // After 1st failed try, rewind - if(!loop) KDS_end("PList_next: End of particle list reached."); - ret = 1; - mcpl_rewind(plist->file); - plist->part = mcpl_read(plist->file); - if(plist->part == NULL) - KDS_error("Error in PList_next: Could not get particle"); - } - if(plist->part->weight>0 && plist->part->pdgcode==plist->pdgcode) return ret; - } - KDS_error("Error in PList_next: MAX_RESAMPLES reached."); - return 1; +int PList_next(PList *plist, int loop) { + int ret = 0, resamples = 0; + while (resamples++ < MAX_RESAMPLES) { + plist->part = mcpl_read(plist->file); + if (plist->part == NULL) { // After 1st failed try, rewind + if (!loop) + KDS_end("PList_next: End of particle list reached."); + ret = 1; + mcpl_rewind(plist->file); + plist->part = mcpl_read(plist->file); + if (plist->part == NULL) + KDS_error("Error in PList_next: Could not get particle"); + } + if (plist->part->weight > 0 && plist->part->pdgcode == plist->pdgcode) + return ret; + } + KDS_error("Error in PList_next: MAX_RESAMPLES reached."); + return 1; } -void PList_destroy(PList* plist){ - free(plist->filename); - mcpl_close_file(plist->file); - free(plist->trasl); - free(plist->rot); - free(plist); +void PList_destroy(PList *plist) { + free(plist->filename); + mcpl_close_file(plist->file); + free(plist->trasl); + free(plist->rot); + free(plist); } diff --git a/src/kdsource/plist.h b/src/kdsource/plist.h index e392bdb..726ff53 100644 --- a/src/kdsource/plist.h +++ b/src/kdsource/plist.h @@ -1,44 +1,43 @@ #ifndef PLISTS_H #define PLISTS_H -#include -#include +#include +#include -#include +#include /***********************************************************************************/ /* */ -/* Utilities for particle lists handling for KDSource. */ +/* Utilities for particle lists handling for KDSource. */ /* */ -/* This file can be freely used as per the terms in the LICENSE file. */ +/* This file can be freely used as per the terms in the LICENSE file. */ /* */ -/* Written by Osiris Inti Abbate, 2021. */ +/* Written by Osiris Inti Abbate, 2021. */ /* */ /***********************************************************************************/ #include "mcpl.h" +typedef struct PList { + char pt; // Particle type ("n", "p", "e", ...) + int pdgcode; // PDG code for particle type -typedef struct PList{ - char pt; // Particle type ("n", "p", "e", ...) - int pdgcode; // PDG code for particle type + char *filename; // Name of MCPL file + mcpl_file_t file; // MCPL file - char* filename; // Name of MCPL file - mcpl_file_t file; // MCPL file + double *trasl; // PList translation + double *rot; // PList rotation + int x2z; // If true, apply permutation x,y,z -> y,z,x - double* trasl; // PList translation - double* rot; // PList rotation - int x2z; // If true, apply permutation x,y,z -> y,z,x - - const mcpl_particle_t* part; // Pointer to selected particle - long long npts; + const mcpl_particle_t *part; // Pointer to selected particle + long long npts; } PList; -PList* PList_create(char pt, const char* filename, const double* trasl, const double* rot, int switch_x2z); -PList* PList_copy(const PList* from); -int PList_get(const PList* plist, mcpl_particle_t* part); -int PList_next(PList* plist, int loop); -void PList_destroy(PList* plist); - +PList *PList_create(char pt, const char *filename, const double *trasl, + const double *rot, int switch_x2z); +PList *PList_copy(const PList *from); +int PList_get(const PList *plist, mcpl_particle_t *part); +int PList_next(PList *plist, int loop); +void PList_destroy(PList *plist); #endif diff --git a/src/kdsource/utils.c b/src/kdsource/utils.c index bf593db..88b3e48 100644 --- a/src/kdsource/utils.c +++ b/src/kdsource/utils.c @@ -1,131 +1,206 @@ -#include -#include -#include +#include +#include +#include #include "kdsource.h" - // Sample with normal distribution (x0=0, s=1) -double rand_norm(){ - double y1 = (double)(rand()+1) / ((double)RAND_MAX+1), y2 = rand() / (double)RAND_MAX; - return sqrt(-2 * log(y1)) * cos(2 * M_PI * y2); +double rand_norm() { + double y1 = (double)(rand() + 1) / ((double)RAND_MAX + 1), + y2 = rand() / (double)RAND_MAX; + return sqrt(-2 * log(y1)) * cos(2 * M_PI * y2); } // Sample with Epanechnikov distribution (x0=0, s=1) -double rand_epan(){ - double x1=(double)rand()/((double)RAND_MAX/2.0) - 1.0; - double x2=(double)rand()/((double)RAND_MAX/2.0) - 1.0; - double x3=(double)rand()/((double)RAND_MAX/2.0) - 1.0; - if(abs(x3)>=abs(x2) && abs(x3)>=abs(x1)) return x2; - else return x3; +double rand_epan() { + double x1 = (double)rand() / ((double)RAND_MAX / 2.0) - 1.0; + double x2 = (double)rand() / ((double)RAND_MAX / 2.0) - 1.0; + double x3 = (double)rand() / ((double)RAND_MAX / 2.0) - 1.0; + if (abs(x3) >= abs(x2) && abs(x3) >= abs(x1)) + return x2; + else + return x3; } // Sample with Tophat distribution (x0=0, s=1) -double rand_box(){ - double y1 = (double)rand()/((double)RAND_MAX/2.0) - 1.0; - return y1; +double rand_box() { + double y1 = (double)rand() / ((double)RAND_MAX / 2.0) - 1.0; + return y1; } // Sample depending on kernel selected -double rand_type(char kernel){ - if(kernel=='g') return rand_norm(); - if(kernel=='e') return rand_epan(); - if(kernel=='b') return rand_box(); - else { - printf("Cannot perturbate with current kernel. \n"); - // KDS_error("Cannot perturbate with current kernel."); - return 0; - } - return 0; +double rand_type(char kernel) { + if (kernel == 'g') + return rand_norm(); + if (kernel == 'e') + return rand_epan(); + if (kernel == 'b') + return rand_box(); + else { + printf("Cannot perturbate with current kernel. \n"); + // KDS_error("Cannot perturbate with current kernel."); + return 0; + } + return 0; } // Translation -double *traslv(double *vect, const double *trasl, int inverse){ - int i; - if(!inverse) for(i=0; i<3; i++) vect[i] += trasl[i]; - else for(i=0; i<3; i++) vect[i] -= trasl[i]; - return vect; +double *traslv(double *vect, const double *trasl, int inverse) { + int i; + if (!inverse) + for (i = 0; i < 3; i++) + vect[i] += trasl[i]; + else + for (i = 0; i < 3; i++) + vect[i] -= trasl[i]; + return vect; } // Rotation (rotvec in axis-angle format) -double *rotv(double *vect, const double *rotvec, int inverse){ - double theta, kdotv, kcrossv[3], vrot[3]; - int i; - theta = sqrt(rotvec[0]*rotvec[0]+rotvec[1]*rotvec[1]+rotvec[2]*rotvec[2]); - if(theta == 0) return vect; - if(inverse) theta = -theta; - - kdotv = rotvec[0]*vect[0]+rotvec[1]*vect[1]+rotvec[2]*vect[2]; - kcrossv[0] = rotvec[1]*vect[2]-rotvec[2]*vect[1]; - kcrossv[1] = rotvec[2]*vect[0]-rotvec[0]*vect[2]; - kcrossv[2] = rotvec[0]*vect[1]-rotvec[1]*vect[0]; - - for(i=0; i<3; i++) vrot[i] = vect[i]*cos(theta) + (kcrossv[i]/fabs(theta))*sin(theta) + rotvec[i]*kdotv/(theta*theta)*(1-cos(theta)); - for(i=0; i<3; i++) vect[i] = vrot[i]; - return vect; +double *rotv(double *vect, const double *rotvec, int inverse) { + double theta, kdotv, kcrossv[3], vrot[3]; + int i; + theta = sqrt(rotvec[0] * rotvec[0] + rotvec[1] * rotvec[1] + + rotvec[2] * rotvec[2]); + if (theta == 0) + return vect; + if (inverse) + theta = -theta; + + kdotv = rotvec[0] * vect[0] + rotvec[1] * vect[1] + rotvec[2] * vect[2]; + kcrossv[0] = rotvec[1] * vect[2] - rotvec[2] * vect[1]; + kcrossv[1] = rotvec[2] * vect[0] - rotvec[0] * vect[2]; + kcrossv[2] = rotvec[0] * vect[1] - rotvec[1] * vect[0]; + + for (i = 0; i < 3; i++) + vrot[i] = vect[i] * cos(theta) + (kcrossv[i] / fabs(theta)) * sin(theta) + + rotvec[i] * kdotv / (theta * theta) * (1 - cos(theta)); + for (i = 0; i < 3; i++) + vect[i] = vrot[i]; + return vect; } // Conversion between char pt and pdg code for particle type -int pt2pdg(char pt){ - if(pt == 'n') return 2112; - if(pt == 'p') return 22; - if(pt == 'e') return 11; - return 0; +int pt2pdg(char pt) { + if (pt == 'n') + return 2112; + if (pt == 'p') + return 22; + if (pt == 'e') + return 11; + return 0; } -char pdg2pt(int pdgcode){ - if(pdgcode == 2112) return 'n'; - if(pdgcode == 22) return 'p'; - if(pdgcode == 11) return 'e'; - return '0'; +char pdg2pt(int pdgcode) { + if (pdgcode == 2112) + return 'n'; + if (pdgcode == 22) + return 'p'; + if (pdgcode == 11) + return 'e'; + return '0'; } // Interpolation -double interp(double x, const double *xs, const double *ys, int N){ - if(xxs[N-1]){ - printf("Error in interp: x outside range.\n"); - return 0; - } - int i = 0; - while(x>xs[i+1]) i++; - return ys[i] + (ys[i+1]-ys[i])*(x-xs[i])/(xs[i+1]-xs[i]); +double interp(double x, const double *xs, const double *ys, int N) { + if (x < xs[0] || x > xs[N - 1]) { + printf("Error in interp: x outside range.\n"); + return 0; + } + int i = 0; + while (x > xs[i + 1]) + i++; + return ys[i] + (ys[i + 1] - ys[i]) * (x - xs[i]) / (xs[i + 1] - xs[i]); } // Dosimetric factors [pSv cm2] #define N_n_ARN 15 -const double log_E_n_ARN[N_n_ARN] = {-1e3,-1.75e+01,-6.21e+00,-3.69e+00,-1.94e+00,-1.39e+00,-5.71e-01,1.82e-01,9.16e-01,1.03e+00,1.16e+00,1.61e+00,2.69e+00,2.94e+00,3.91e+00}; -const double log_f_n_ARN[N_n_ARN] = {-1e3,2.36e+00,2.04e+00,2.96e+00,4.84e+00,5.31e+00,5.84e+00,6.05e+00,6.03e+00,6.02e+00,6.02e+00,6.00e+00,6.28e+00,6.37e+00,6.37e+00}; +const double log_E_n_ARN[N_n_ARN] = {-1e3, -1.75e+01, -6.21e+00, -3.69e+00, + -1.94e+00, -1.39e+00, -5.71e-01, 1.82e-01, + 9.16e-01, 1.03e+00, 1.16e+00, 1.61e+00, + 2.69e+00, 2.94e+00, 3.91e+00}; +const double log_f_n_ARN[N_n_ARN] = {-1e3, 2.36e+00, 2.04e+00, 2.96e+00, + 4.84e+00, 5.31e+00, 5.84e+00, 6.05e+00, + 6.03e+00, 6.02e+00, 6.02e+00, 6.00e+00, + 6.28e+00, 6.37e+00, 6.37e+00}; #define N_p_ARN 27 -const double log_E_p_ARN[N_p_ARN] = {-1e3,-4.61e+00,-4.20e+00,-3.91e+00,-3.51e+00,-3.22e+00,-3.00e+00,-2.81e+00,-2.53e+00,-2.30e+00,-1.90e+00,-1.61e+00,-1.20e+00,-9.16e-01,-6.93e-01,-5.11e-01,-2.23e-01,0.00e+00,4.05e-01,6.93e-01,1.10e+00,1.39e+00,1.61e+00,1.79e+00,2.08e+00,2.30e+00,3.91e+00}; -const double log_f_p_ARN[N_p_ARN] = {-1e3,-2.80e+00,-1.86e-01,4.88e-02,-2.11e-01,-4.46e-01,-5.98e-01,-6.73e-01,-6.35e-01,-4.94e-01,-1.17e-01,1.82e-01,5.88e-01,8.67e-01,1.08e+00,1.24e+00,1.48e+00,1.65e+00,1.93e+00,2.15e+00,2.41e+00,2.60e+00,2.74e+00,2.87e+00,3.07e+00,3.24e+00,3.24e+00}; +const double log_E_p_ARN[N_p_ARN] = { + -1e3, -4.61e+00, -4.20e+00, -3.91e+00, -3.51e+00, -3.22e+00, -3.00e+00, + -2.81e+00, -2.53e+00, -2.30e+00, -1.90e+00, -1.61e+00, -1.20e+00, -9.16e-01, + -6.93e-01, -5.11e-01, -2.23e-01, 0.00e+00, 4.05e-01, 6.93e-01, 1.10e+00, + 1.39e+00, 1.61e+00, 1.79e+00, 2.08e+00, 2.30e+00, 3.91e+00}; +const double log_f_p_ARN[N_p_ARN] = { + -1e3, -2.80e+00, -1.86e-01, 4.88e-02, -2.11e-01, -4.46e-01, -5.98e-01, + -6.73e-01, -6.35e-01, -4.94e-01, -1.17e-01, 1.82e-01, 5.88e-01, 8.67e-01, + 1.08e+00, 1.24e+00, 1.48e+00, 1.65e+00, 1.93e+00, 2.15e+00, 2.41e+00, + 2.60e+00, 2.74e+00, 2.87e+00, 3.07e+00, 3.24e+00, 3.24e+00}; #define N_n_ICRP 69 -const double log_E_n_ICRP[N_n_ICRP] = {-1e3,-2.07e+01,-1.84e+01,-1.75e+01,-1.61e+01,-1.54e+01,-1.45e+01,-1.38e+01,-1.31e+01,-1.22e+01,-1.15e+01,-1.08e+01,-9.90e+00,-9.21e+00,-8.52e+00,-7.60e+00,-6.91e+00,-6.21e+00,-5.30e+00,-4.61e+00,-3.91e+00,-3.51e+00,-3.00e+00,-2.66e+00,-2.30e+00,-1.90e+00,-1.61e+00,-1.20e+00,-6.93e-01,-3.57e-01,-1.05e-01,0.00e+00,1.82e-01,4.05e-01,6.93e-01,1.10e+00,1.39e+00,1.61e+00,1.79e+00,1.95e+00,2.08e+00,2.20e+00,2.30e+00,2.48e+00,2.64e+00,2.71e+00,2.77e+00,2.89e+00,3.00e+00,3.04e+00,3.40e+00,3.91e+00,4.32e+00,4.61e+00,4.87e+00,5.01e+00,5.19e+00,5.30e+00,5.70e+00,5.99e+00,6.21e+00,6.40e+00,6.55e+00,6.68e+00,6.80e+00,6.91e+00,7.60e+00,8.52e+00,9.21e+00}; -const double log_f_n_ICRP[N_n_ICRP] = {-1e3,1.13e+00,1.27e+00,1.39e+00,1.65e+00,1.77e+00,1.89e+00,1.95e+00,2.00e+00,2.04e+00,2.06e+00,2.06e+00,2.06e+00,2.05e+00,2.05e+00,2.02e+00,2.02e+00,2.03e+00,2.08e+00,2.21e+00,2.50e+00,2.75e+00,3.14e+00,3.42e+00,3.74e+00,4.10e+00,4.37e+00,4.74e+00,5.18e+00,5.45e+00,5.63e+00,5.71e+00,5.80e+00,5.90e+00,6.01e+00,6.13e+00,6.18e+00,6.20e+00,6.21e+00,6.21e+00,6.21e+00,6.21e+00,6.21e+00,6.21e+00,6.20e+00,6.20e+00,6.19e+00,6.18e+00,6.17e+00,6.16e+00,6.12e+00,6.07e+00,6.08e+00,6.10e+00,6.10e+00,6.10e+00,6.10e+00,6.10e+00,6.16e+00,6.24e+00,6.28e+00,6.34e+00,6.44e+00,6.46e+00,6.47e+00,6.50e+00,6.65e+00,6.95e+00,7.24e+00}; +const double log_E_n_ICRP[N_n_ICRP] = { + -1e3, -2.07e+01, -1.84e+01, -1.75e+01, -1.61e+01, -1.54e+01, -1.45e+01, + -1.38e+01, -1.31e+01, -1.22e+01, -1.15e+01, -1.08e+01, -9.90e+00, -9.21e+00, + -8.52e+00, -7.60e+00, -6.91e+00, -6.21e+00, -5.30e+00, -4.61e+00, -3.91e+00, + -3.51e+00, -3.00e+00, -2.66e+00, -2.30e+00, -1.90e+00, -1.61e+00, -1.20e+00, + -6.93e-01, -3.57e-01, -1.05e-01, 0.00e+00, 1.82e-01, 4.05e-01, 6.93e-01, + 1.10e+00, 1.39e+00, 1.61e+00, 1.79e+00, 1.95e+00, 2.08e+00, 2.20e+00, + 2.30e+00, 2.48e+00, 2.64e+00, 2.71e+00, 2.77e+00, 2.89e+00, 3.00e+00, + 3.04e+00, 3.40e+00, 3.91e+00, 4.32e+00, 4.61e+00, 4.87e+00, 5.01e+00, + 5.19e+00, 5.30e+00, 5.70e+00, 5.99e+00, 6.21e+00, 6.40e+00, 6.55e+00, + 6.68e+00, 6.80e+00, 6.91e+00, 7.60e+00, 8.52e+00, 9.21e+00}; +const double log_f_n_ICRP[N_n_ICRP] = { + -1e3, 1.13e+00, 1.27e+00, 1.39e+00, 1.65e+00, 1.77e+00, 1.89e+00, + 1.95e+00, 2.00e+00, 2.04e+00, 2.06e+00, 2.06e+00, 2.06e+00, 2.05e+00, + 2.05e+00, 2.02e+00, 2.02e+00, 2.03e+00, 2.08e+00, 2.21e+00, 2.50e+00, + 2.75e+00, 3.14e+00, 3.42e+00, 3.74e+00, 4.10e+00, 4.37e+00, 4.74e+00, + 5.18e+00, 5.45e+00, 5.63e+00, 5.71e+00, 5.80e+00, 5.90e+00, 6.01e+00, + 6.13e+00, 6.18e+00, 6.20e+00, 6.21e+00, 6.21e+00, 6.21e+00, 6.21e+00, + 6.21e+00, 6.21e+00, 6.20e+00, 6.20e+00, 6.19e+00, 6.18e+00, 6.17e+00, + 6.16e+00, 6.12e+00, 6.07e+00, 6.08e+00, 6.10e+00, 6.10e+00, 6.10e+00, + 6.10e+00, 6.10e+00, 6.16e+00, 6.24e+00, 6.28e+00, 6.34e+00, 6.44e+00, + 6.46e+00, 6.47e+00, 6.50e+00, 6.65e+00, 6.95e+00, 7.24e+00}; #define N_p_ICRP 65 -const double log_E_p_ICRP[N_p_ICRP] = {-1e3,-5.30e+00,-5.12e+00,-4.96e+00,-4.83e+00,-4.71e+00,-4.61e+00,-4.42e+00,-4.34e+00,-4.20e+00,-4.07e+00,-3.91e+00,-3.69e+00,-3.51e+00,-3.22e+00,-3.00e+00,-2.81e+00,-2.66e+00,-2.53e+00,-2.30e+00,-1.90e+00,-1.61e+00,-1.20e+00,-9.16e-01,-6.93e-01,-6.71e-01,-5.11e-01,-4.12e-01,-2.23e-01,0.00e+00,1.13e-01,2.85e-01,4.05e-01,6.93e-01,1.10e+00,1.39e+00,1.61e+00,1.79e+00,1.81e+00,2.08e+00,2.30e+00,2.71e+00,3.00e+00,3.40e+00,3.69e+00,3.91e+00,4.09e+00,4.38e+00,4.61e+00,5.01e+00,5.30e+00,5.70e+00,5.99e+00,6.21e+00,6.40e+00,6.68e+00,6.91e+00,7.31e+00,7.60e+00,8.01e+00,8.29e+00,8.52e+00,8.70e+00,8.99e+00,9.21e+00}; -const double log_f_p_ICRP[N_p_ICRP] = {-1e3,-4.31e+00,-4.10e+00,-3.79e+00,-3.40e+00,-3.02e+00,-2.68e+00,-2.25e+00,-2.10e+00,-1.86e+00,-1.71e+00,-1.49e+00,-1.29e+00,-1.16e+00,-1.05e+00,-9.97e-01,-9.44e-01,-8.89e-01,-8.14e-01,-6.58e-01,-2.92e-01,0.00e+00,4.12e-01,6.93e-01,9.04e-01,9.24e-01,1.07e+00,1.15e+00,1.32e+00,1.50e+00,1.59e+00,1.72e+00,1.81e+00,2.01e+00,2.28e+00,2.46e+00,2.60e+00,2.71e+00,2.72e+00,2.92e+00,3.10e+00,3.41e+00,3.64e+00,3.94e+00,4.12e+00,4.28e+00,4.41e+00,4.59e+00,4.70e+00,4.87e+00,4.97e+00,5.08e+00,5.15e+00,5.20e+00,5.23e+00,5.28e+00,5.33e+00,5.36e+00,5.46e+00,5.53e+00,5.59e+00,5.62e+00,5.65e+00,5.70e+00,5.73e+00}; - -double H10_n_ARN(double E){ - double log_E = log(E); - double log_f = interp(log_E, log_E_n_ARN, log_f_n_ARN, N_n_ARN); - return exp(log_f); +const double log_E_p_ICRP[N_p_ICRP] = { + -1e3, -5.30e+00, -5.12e+00, -4.96e+00, -4.83e+00, -4.71e+00, -4.61e+00, + -4.42e+00, -4.34e+00, -4.20e+00, -4.07e+00, -3.91e+00, -3.69e+00, -3.51e+00, + -3.22e+00, -3.00e+00, -2.81e+00, -2.66e+00, -2.53e+00, -2.30e+00, -1.90e+00, + -1.61e+00, -1.20e+00, -9.16e-01, -6.93e-01, -6.71e-01, -5.11e-01, -4.12e-01, + -2.23e-01, 0.00e+00, 1.13e-01, 2.85e-01, 4.05e-01, 6.93e-01, 1.10e+00, + 1.39e+00, 1.61e+00, 1.79e+00, 1.81e+00, 2.08e+00, 2.30e+00, 2.71e+00, + 3.00e+00, 3.40e+00, 3.69e+00, 3.91e+00, 4.09e+00, 4.38e+00, 4.61e+00, + 5.01e+00, 5.30e+00, 5.70e+00, 5.99e+00, 6.21e+00, 6.40e+00, 6.68e+00, + 6.91e+00, 7.31e+00, 7.60e+00, 8.01e+00, 8.29e+00, 8.52e+00, 8.70e+00, + 8.99e+00, 9.21e+00}; +const double log_f_p_ICRP[N_p_ICRP] = { + -1e3, -4.31e+00, -4.10e+00, -3.79e+00, -3.40e+00, -3.02e+00, -2.68e+00, + -2.25e+00, -2.10e+00, -1.86e+00, -1.71e+00, -1.49e+00, -1.29e+00, -1.16e+00, + -1.05e+00, -9.97e-01, -9.44e-01, -8.89e-01, -8.14e-01, -6.58e-01, -2.92e-01, + 0.00e+00, 4.12e-01, 6.93e-01, 9.04e-01, 9.24e-01, 1.07e+00, 1.15e+00, + 1.32e+00, 1.50e+00, 1.59e+00, 1.72e+00, 1.81e+00, 2.01e+00, 2.28e+00, + 2.46e+00, 2.60e+00, 2.71e+00, 2.72e+00, 2.92e+00, 3.10e+00, 3.41e+00, + 3.64e+00, 3.94e+00, 4.12e+00, 4.28e+00, 4.41e+00, 4.59e+00, 4.70e+00, + 4.87e+00, 4.97e+00, 5.08e+00, 5.15e+00, 5.20e+00, 5.23e+00, 5.28e+00, + 5.33e+00, 5.36e+00, 5.46e+00, 5.53e+00, 5.59e+00, 5.62e+00, 5.65e+00, + 5.70e+00, 5.73e+00}; + +double H10_n_ARN(double E) { + double log_E = log(E); + double log_f = interp(log_E, log_E_n_ARN, log_f_n_ARN, N_n_ARN); + return exp(log_f); } -double H10_p_ARN(double E){ - double log_E = log(E); - double log_f = interp(log_E, log_E_p_ARN, log_f_p_ARN, N_p_ARN); - return exp(log_f); +double H10_p_ARN(double E) { + double log_E = log(E); + double log_f = interp(log_E, log_E_p_ARN, log_f_p_ARN, N_p_ARN); + return exp(log_f); } -double H10_n_ICRP(double E){ - double log_E = log(E); - double log_f = interp(log_E, log_E_n_ICRP, log_f_n_ICRP, N_n_ICRP); - return exp(log_f); +double H10_n_ICRP(double E) { + double log_E = log(E); + double log_f = interp(log_E, log_E_n_ICRP, log_f_n_ICRP, N_n_ICRP); + return exp(log_f); } -double H10_p_ICRP(double E){ - double log_E = log(E); - double log_f = interp(log_E, log_E_p_ICRP, log_f_p_ICRP, N_p_ICRP); - return exp(log_f); +double H10_p_ICRP(double E) { + double log_E = log(E); + double log_f = interp(log_E, log_E_p_ICRP, log_f_p_ICRP, N_p_ICRP); + return exp(log_f); } - diff --git a/src/kdsource/utils.h b/src/kdsource/utils.h index 42a472f..ac89aa5 100644 --- a/src/kdsource/utils.h +++ b/src/kdsource/utils.h @@ -1,20 +1,19 @@ #ifndef AUX_H #define AUX_H -#include -#include +#include +#include /***********************************************************************************/ /* */ -/* General utilities for KDSource. */ +/* General utilities for KDSource. */ /* */ -/* This file can be freely used as per the terms in the LICENSE file. */ +/* This file can be freely used as per the terms in the LICENSE file. */ /* */ -/* Written by Osiris Inti Abbate, 2021. */ +/* Written by Osiris Inti Abbate, 2021. */ /* */ /***********************************************************************************/ - double rand_norm(); double rand_epan(); double rand_box(); @@ -33,5 +32,4 @@ double H10_p_ARN(double E); double H10_n_ICRP(double E); double H10_p_ICRP(double E); - #endif diff --git a/src/kdtool/beamtest.c b/src/kdtool/beamtest.c index b038608..1193d37 100644 --- a/src/kdtool/beamtest.c +++ b/src/kdtool/beamtest.c @@ -1,206 +1,236 @@ -#include -#include -#include +#include +#include +#include -#include +#include #include "kdsource.h" - -void display_usage(){ - printf("Usage: kdtool beamtest sourcefile [options]\n\n"); - printf("Executes a simple simulation with source defined in XML file sourcefile, in\n"); - printf("which calculates the number of particles passing thru a rectangular collimator.\n"); - printf("The simulation is repeated using the particle list directly as source, to\n"); - printf("compare the results.\n\n"); - printf("Results are computed for 4 energy groups, and stored in a results file which\n"); - printf("can be imported from a spreadsheet.\n\n"); - printf("This tool is designed to be used with flat neutron sources with particles\n"); - printf("propagating towards z direction.\n\n"); - printf("Options:\n"); - printf("\t-n N: Number of source particles (default: 1E6).\n"); - printf("\t-o results: Name of file to store results.\n"); - printf("\t-xwidth value: Width of the collimator, in cm (default: 7).\n"); - printf("\t-yheight value: Height of the collimator, in cm (default: 20).\n"); - printf("\t-z value: Position of the collimator along z axis, in cm\n"); - printf("\t (default: 500).\n"); - printf("\t-xshift: Horizontal shift of the center of the collimator,\n"); - printf("\t in cm (default: 0)\n"); - printf("\t-yshift: Vertical shift of the center of the collimator,\n"); - printf("\t in cm (default: 0)\n"); - printf("\t-h, --help: Display usage instructions.\n"); - +void display_usage() { + printf("Usage: kdtool beamtest sourcefile [options]\n\n"); + printf("Executes a simple simulation with source defined in XML file " + "sourcefile, in\n"); + printf("which calculates the number of particles passing thru a rectangular " + "collimator.\n"); + printf("The simulation is repeated using the particle list directly as " + "source, to\n"); + printf("compare the results.\n\n"); + printf("Results are computed for 4 energy groups, and stored in a results " + "file which\n"); + printf("can be imported from a spreadsheet.\n\n"); + printf("This tool is designed to be used with flat neutron sources with " + "particles\n"); + printf("propagating towards z direction.\n\n"); + printf("Options:\n"); + printf("\t-n N: Number of source particles (default: 1E6).\n"); + printf("\t-o results: Name of file to store results.\n"); + printf("\t-xwidth value: Width of the collimator, in cm (default: 7).\n"); + printf("\t-yheight value: Height of the collimator, in cm (default: 20).\n"); + printf("\t-z value: Position of the collimator along z axis, in cm\n"); + printf("\t (default: 500).\n"); + printf( + "\t-xshift: Horizontal shift of the center of the collimator,\n"); + printf("\t in cm (default: 0)\n"); + printf("\t-yshift: Vertical shift of the center of the collimator,\n"); + printf("\t in cm (default: 0)\n"); + printf("\t-h, --help: Display usage instructions.\n"); } -int beamtest_parse_args(int argc, char **argv, const char** sourcefile, - long* N, const char** results, double* xwidth, double* yheight, double* z, double* xshift, double* yshift){ - *sourcefile = 0; - *N = 1E6; - *results = "results_beamtest.txt"; - *xwidth = 7; - *yheight = 20; - *z = 500; - *xshift = 0; - *yshift = 0; - int i; - for(i=1; iposition[2] - z) * part->direction[2] > 0) - return 0; - if(part->ekin <= 0) - return 0; - double t = (z-part->position[2]) / part->direction[2]; // Time to reach collimator assuming v=1 - double x = part->position[0] + t * part->direction[0]; - double y = part->position[1] + t * part->direction[1]; - return (fabs(x-xshift)position[2] - z) * part->direction[2] > 0) + return 0; + if (part->ekin <= 0) + return 0; + double t = (z - part->position[2]) / + part->direction[2]; // Time to reach collimator assuming v=1 + double x = part->position[0] + t * part->direction[0]; + double y = part->position[1] + t * part->direction[1]; + return (fabs(x - xshift) < xwidth / 2 && fabs(y - yshift) < yheight / 2); } -int main(int argc, char *argv[]){ - // User defined options - const char *sourcefile, *results; - double xwidth, yheight, z, xshift, yshift; - long int N; - if(beamtest_parse_args(argc, argv, &sourcefile, &N, &results, &xwidth, &yheight, &z, &xshift, &yshift)) return 1; +int main(int argc, char *argv[]) { + // User defined options + const char *sourcefile, *results; + double xwidth, yheight, z, xshift, yshift; + long int N; + if (beamtest_parse_args(argc, argv, &sourcefile, &N, &results, &xwidth, + &yheight, &z, &xshift, &yshift)) + return 1; + + // KDSource object + KDSource *kds = KDS_open(sourcefile); + mcpl_particle_t part; + double w_crit = KDS_w_mean(kds, 1000, NULL); - // KDSource object - KDSource* kds = KDS_open(sourcefile); - mcpl_particle_t part; - double w_crit = KDS_w_mean(kds, 1000, NULL); + // Variables to store results + long int N_source_KDE = 0, N_source_tracks = 0; + double I_source_KDE = 0, I_source_tracks = 0, I_KDE_tot = 0, I_tracks_tot = 0, + err_KDE_tot = 0, err_tracks_tot = 0; + int ngroups = 4; + double Ecuts[] = {0, 1E-8, 3E-7, 0.5, 100}; // Energy groups cuts + double I_KDE[] = {0, 0, 0, 0}, I_tracks[] = {0, 0, 0, 0}, + err_KDE[] = {0, 0, 0, 0}, err_tracks[] = {0, 0, 0, 0}; + long int i; + int j; - // Variables to store results - long int N_source_KDE=0, N_source_tracks=0; - double I_source_KDE=0, I_source_tracks=0, I_KDE_tot=0, I_tracks_tot=0, err_KDE_tot=0, err_tracks_tot=0; - int ngroups = 4; - double Ecuts[] = {0, 1E-8, 3E-7, 0.5, 100}; // Energy groups cuts - double I_KDE[] = {0,0,0,0}, I_tracks[] = {0,0,0,0}, err_KDE[] = {0,0,0,0}, err_tracks[] = {0,0,0,0}; - long int i; - int j; + // KDE simulation + mcpl_rewind(kds->plist->file); + if (kds->geom->bwfile) + rewind(kds->geom->bwfile); + printf("Executing simulation with KDE source...\n"); + for (i = 0; i < N; i++) { + KDS_sample2(kds, &part, 1, w_crit, NULL, 1); + N_source_KDE++; + I_source_KDE += part.weight; + if (propagate_collimator(&part, xwidth, yheight, z, xshift, yshift)) { + for (j = 0; j < ngroups; j++) { + if (part.ekin > Ecuts[j] && part.ekin < Ecuts[j + 1]) { + I_KDE[j] += part.weight; + err_KDE[j] += part.weight * part.weight; + } + } + } + } + printf("Finished. Particles simulated: N_source = %ld, I_source = %lf\n\n", + N_source_KDE, I_source_KDE); - // KDE simulation - mcpl_rewind(kds->plist->file); - if(kds->geom->bwfile) rewind(kds->geom->bwfile); - printf("Executing simulation with KDE source...\n"); - for(i=0; i Ecuts[j] && part.ekin < Ecuts[j+1]){ - I_KDE[j] += part.weight; - err_KDE[j] += part.weight * part.weight; - } - } - } - } - printf("Finished. Particles simulated: N_source = %ld, I_source = %lf\n\n", N_source_KDE, I_source_KDE); + // Tracks simulation + mcpl_rewind(kds->plist->file); + if (kds->geom->bwfile) + rewind(kds->geom->bwfile); + printf("Executing simulation with tracks source...\n"); + for (i = 0; i < N; i++) { + if (KDS_sample2(kds, &part, 0, w_crit, NULL, 1)) + break; + N_source_tracks++; + I_source_tracks += part.weight; + if (propagate_collimator(&part, xwidth, yheight, z, xshift, yshift)) { + for (j = 0; j < ngroups; j++) { + if (part.ekin > Ecuts[j] && part.ekin < Ecuts[j + 1]) { + I_tracks[j] += part.weight; + err_tracks[j] += part.weight * part.weight; + } + } + } + } + printf("Finished. Particles simulated: N_source = %ld, I_source = %lf\n\n", + N_source_tracks, I_source_tracks); - // Tracks simulation - mcpl_rewind(kds->plist->file); - if(kds->geom->bwfile) rewind(kds->geom->bwfile); - printf("Executing simulation with tracks source...\n"); - for(i=0; i Ecuts[j] && part.ekin < Ecuts[j+1]){ - I_tracks[j] += part.weight; - err_tracks[j] += part.weight * part.weight; - } - } - } - } - printf("Finished. Particles simulated: N_source = %ld, I_source = %lf\n\n", N_source_tracks, I_source_tracks); + // Normalize results + for (j = 0; j < ngroups; j++) { + I_KDE[j] /= I_source_KDE; + I_tracks[j] /= I_source_tracks; + I_KDE_tot += I_KDE[j]; + I_tracks_tot += I_tracks[j]; + err_KDE_tot += err_KDE[j]; + err_tracks_tot += err_tracks[j]; + err_KDE[j] = sqrt(err_KDE[j]) / I_source_KDE; + err_tracks[j] = sqrt(err_tracks[j]) / I_source_tracks; + } + err_KDE_tot = sqrt(err_KDE_tot) / I_source_KDE; + err_tracks_tot = sqrt(err_tracks_tot) / I_source_tracks; - // Normalize results - for(j=0; j -#include -#include +#include +#include +#include -#include +#include #include "kdsource.h" - -void display_usage(){ - printf("Usage: kdtool resample sourcefile [options]\n\n"); - printf("Resample particles from source defined in XML file sourcefile, and save them in\n"); - printf("a MCPL file.\n\n"); - printf("Options:\n"); - printf("\t-o outfile: Name of MCPL file with new samples\n"); - printf("\t (default: \"resampled.mcpl\").\n"); - printf("\t-n N: Number of new samples (default: 1E5).\n"); - printf("\t-h, --help: Display usage instructions.\n"); - +void display_usage() { + printf("Usage: kdtool resample sourcefile [options]\n\n"); + printf("Resample particles from source defined in XML file sourcefile, and " + "save them in\n"); + printf("a MCPL file.\n\n"); + printf("Options:\n"); + printf("\t-o outfile: Name of MCPL file with new samples\n"); + printf("\t (default: \"resampled.mcpl\").\n"); + printf("\t-n N: Number of new samples (default: 1E5).\n"); + printf("\t-h, --help: Display usage instructions.\n"); } -int resample_parse_args(int argc, char **argv, const char** filename, const char** outfilename, long int* N){ - *filename = 0; - *outfilename = 0; - *N = 1E5; - int i; - for(i=1; i