Skip to content

Commit 4d10d1c

Browse files
marinerhemantmarinerhemant
committed
Propagate higher-order distortion model (p7-p10) across FF_HEDM suite
This commit fully propagates the dipole (`p7`, `p8`) and trefoil (`p9`, `p10`) radial distortion parameters across the remaining elements of the MIDAS FF_HEDM optimization suite and forward simulation modules. The goal is to seamlessly introduce higher-order spatial harmonic corrections to the detector transform mathematically without breaking layer compatibility with older runs. It resolves the full 17-parameter geometric distortion state natively. Key structural changes include: - `PeaksFittingOMPZarrRefactor.c`: Added Zarr parsing for p7-p10 keys and extended the internal `DistortFunc` to inject the dipole and trefoil vectors. - `FitMultipleGrains.c`: Expanded global variable `OptP` arrays from 13 to 17 constraints, adjusted global index offsets, and safely updated `CorrectTiltSpatialDistortion` signatures. - `FitWedgeParallel.c` & `FitSetupParamsAllZarr.c`: Plumbed parameter block reads from `ZarrReader`, passing `p7`–`p10` correctly through `YsZsCalc` so that theoretical models conform to the non-linear detector field. - `ForwardSimulationCompressed.c`: Overhauled the spatial projection signatures inside `CorrectTiltSpatialDistortion` so forward-simulated sinograms map to the higher-harmonics accurately. - Deprecated `CalcPeakProfile.c`: This file was found to be entirely obsolete to the active pipeline (only tied to legacy `.c` codes) and has been moved to `src/archive/` to keep `FF_HEDM/src` clean. References in `CMakeLists.txt` have been safely relinked to the archive.
1 parent 19540ef commit 4d10d1c

8 files changed

Lines changed: 114 additions & 31 deletions

‎FF_HEDM/CMakeLists.txt‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -77,9 +77,9 @@ endfunction()
7777
# Most executables link against COMMON_LINK_LIBRARIES implicitly via add_ff_hedm_executable.
7878
# Any libraries listed in LINK_LIBRARIES are additional, target-specific dependencies.
7979

80-
# add_ff_hedm_executable(Calibrant SOURCES src/Calibrant.c src/CalcPeakProfile.c)
81-
# add_ff_hedm_executable(CalibrantOMP SOURCES src/CalibrantOMP.c src/CalcPeakProfile.c src/MIDAS_Math.c OMP) # Archived: CalibrantPanelShiftsOMP is the active superset
82-
# add_ff_hedm_executable(CalibrantPanelShiftsOMP SOURCES src/CalibrantPanelShiftsOMP.c src/CalcPeakProfile.c src/DetectorGeometry.c src/Panel.c src/FileReader.c src/MIDAS_Math.c src/MIDAS_ParamParser.c OMP) # Archived: CalibrantIntegratorOMP is the active replacement
80+
# add_ff_hedm_executable(Calibrant SOURCES src/archive/Calibrant.c src/archive/CalcPeakProfile.c)
81+
# add_ff_hedm_executable(CalibrantOMP SOURCES src/archive/CalibrantOMP.c src/archive/CalcPeakProfile.c src/MIDAS_Math.c OMP) # Archived: CalibrantPanelShiftsOMP is the active superset
82+
# add_ff_hedm_executable(CalibrantPanelShiftsOMP SOURCES src/archive/CalibrantPanelShiftsOMP.c src/archive/CalcPeakProfile.c src/DetectorGeometry.c src/Panel.c src/FileReader.c src/MIDAS_Math.c src/MIDAS_ParamParser.c OMP) # Archived: CalibrantIntegratorOMP is the active replacement
8383
add_ff_hedm_executable(CalibrantIntegratorOMP SOURCES src/CalibrantIntegratorOMP.c src/DetectorGeometry.c src/Panel.c src/FileReader.c src/MIDAS_Math.c src/MIDAS_ParamParser.c OMP)
8484
target_link_libraries(CalibrantIntegratorOMP PRIVATE midas_integration)
8585
# add_ff_hedm_executable(FitSetup SOURCES src/FitSetupParamsAll.c)

‎FF_HEDM/src/CalibrationCore.c‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,8 @@
99
//
1010

1111
#include "CalibrationCore.h"
12+
#include "MIDAS_Math.h"
13+
#include <nlopt.h>
1214
#include <math.h>
1315
#include <stdio.h>
1416
#include <stdlib.h>

‎FF_HEDM/src/FitMultipleGrains.c‎

Lines changed: 34 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -527,6 +527,7 @@ static inline void CorrectTiltSpatialDistortion(
527527
int nIndices, double MaxRad, double **SpotInfoAll, double px, double Lsd,
528528
double ybc, double zbc, double tx, double ty, double tz, double p0,
529529
double p1, double p2, double p3, double p4, double p5, double p6,
530+
double p7, double p8, double p9, double p10,
530531
double **SpotInfoCorr, int nPanels, Panel *panels) {
531532
double txr, tyr, tzr;
532533
txr = deg2rad * tx;
@@ -567,10 +568,13 @@ static inline void CorrectTiltSpatialDistortion(
567568
Eta = CalcEtaAngleLocal(XYZ[1], XYZ[2]);
568569
RNorm = Rad / MaxRad;
569570
EtaT = 90 - Eta;
571+
double RNorm3 = RNorm * RNorm * RNorm;
572+
double dipole = p7 * RNorm * cos(deg2rad * (EtaT + p8));
573+
double trefoil = p9 * RNorm3 * cos(deg2rad * (3.0 * EtaT + p10));
570574
DistortFunc = (p0 * (pow(RNorm, n0)) * (cos(deg2rad * (2 * EtaT + p6)))) +
571575
(p1 * (pow(RNorm, n1)) * (cos(deg2rad * (4 * EtaT + p3)))) +
572576
(panelP2 * (pow(RNorm, n2))) + p4 * pow(RNorm, 6.0) +
573-
p5 * pow(RNorm, 4.0) + 1;
577+
p5 * pow(RNorm, 4.0) + dipole + trefoil + 1;
574578
Rcorr = Rad * DistortFunc;
575579
Rcorr = Rcorr * (Lsd / panelLsd); // re-project to global Lsd plane
576580
SpotInfoCorr[i][0] = SpotInfoAll[i][0];
@@ -631,6 +635,10 @@ static void calc_grain_errors(int nGrains, int nPanels, Panel *panels,
631635
double p4 = x[global_offset + 10];
632636
double p5 = x[global_offset + 11];
633637
double p6 = x[global_offset + 12];
638+
double p7 = x[global_offset + 13];
639+
double p8 = x[global_offset + 14];
640+
double p9 = x[global_offset + 15];
641+
double p10 = x[global_offset + 16];
634642

635643
for (int g = 0; g < nGrains; g++) {
636644
struct GrainData *g_data = &grains[g];
@@ -665,7 +673,7 @@ static void calc_grain_errors(int nGrains, int nPanels, Panel *panels,
665673

666674
CorrectTiltSpatialDistortion(g_data->nSpots, g_data->RhoD, SpotInfoShifted,
667675
g_data->px, g_data->Lsd, ybc, zbc, tx, ty, tz,
668-
p0, p1, p2, p3, p4, p5, p6, SpotInfoCorr,
676+
p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, SpotInfoCorr,
669677
nPanels, panels);
670678
FreeMemMatrix(SpotInfoShifted, g_data->nSpots);
671679

@@ -708,10 +716,14 @@ static double problem_function(unsigned n, const double *x, double *grad,
708716
double p4 = x[global_offset + 10];
709717
double p5 = x[global_offset + 11];
710718
double p6 = x[global_offset + 12];
719+
double p7 = x[global_offset + 13];
720+
double p8 = x[global_offset + 14];
721+
double p9 = x[global_offset + 15];
722+
double p10 = x[global_offset + 16];
711723

712724
// Panels are at the end of x
713725
if (nPanels > 1) {
714-
int p_idx = global_offset + 13;
726+
int p_idx = global_offset + 17;
715727
for (int i = 0; i < nPanels; i++) {
716728
if (i == f_data->fixPanel) {
717729
panels[i].dY = 0;
@@ -765,7 +777,7 @@ static double problem_function(unsigned n, const double *x, double *grad,
765777

766778
CorrectTiltSpatialDistortion(g_data->nSpots, g_data->RhoD, SpotInfoShifted,
767779
g_data->px, g_data->Lsd, ybc, zbc, tx, ty, tz,
768-
p0, p1, p2, p3, p4, p5, p6, SpotInfoCorr,
780+
p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, SpotInfoCorr,
769781
f_data->nPanels, f_data->panels);
770782
FreeMemMatrix(SpotInfoShifted, g_data->nSpots);
771783

@@ -794,11 +806,11 @@ static double problem_function(unsigned n, const double *x, double *grad,
794806
return totalError;
795807
}
796808

797-
void FitMultipleGrains(struct GrainData *grains, int nGrains, double OptP[13],
798-
double NonOptP[5], double tols[26], double *Out,
809+
void FitMultipleGrains(struct GrainData *grains, int nGrains, double OptP[17],
810+
double NonOptP[5], double tols[30], double *Out,
799811
Panel *panels, int nPanels, double tolShifts,
800812
int *spotsPerPanel, int fixPanel) {
801-
unsigned n = nGrains * 12 + 13; // Grains(12 each) + Global(13)
813+
unsigned n = nGrains * 12 + 17; // Grains(12 each) + Global(17)
802814
if (nPanels > 1) {
803815
n += (nPanels - 1) * 2;
804816
}
@@ -845,16 +857,16 @@ void FitMultipleGrains(struct GrainData *grains, int nGrains, double OptP[13],
845857

846858
// 2. Set Global Parameters
847859
int global_offset = nGrains * 12;
848-
for (int i = 0; i < 13; i++) {
860+
for (int i = 0; i < 17; i++) {
849861
x[global_offset + i] = OptP[i];
850-
// Tols for global params are at indices 12-24 in 'tols' array
862+
// Tols for global params are at indices 12-28 in 'tols' array
851863
xl[global_offset + i] = x[global_offset + i] - tols[12 + i];
852864
xu[global_offset + i] = x[global_offset + i] + tols[12 + i];
853865
}
854866

855867
// 3. Set Panel Parameters
856868
if (nPanels > 1) {
857-
int p_idx = global_offset + 13; // Start of panel parameters
869+
int p_idx = global_offset + 17; // Start of panel parameters
858870
for (int i = 0; i < nPanels; i++) {
859871
if (i == fixPanel) { // Skip the fixed panel
860872
continue;
@@ -924,6 +936,10 @@ void FitMultipleGrains(struct GrainData *grains, int nGrains, double OptP[13],
924936
printf("p4 %.12f\n", gOpt[10]);
925937
printf("p5 %.12f\n", gOpt[11]);
926938
printf("p6 %.12f\n", gOpt[12]);
939+
printf("p7 %.12f\n", gOpt[13]);
940+
printf("p8 %.12f\n", gOpt[14]);
941+
printf("p9 %.12f\n", gOpt[15]);
942+
printf("p10 %.12f\n", gOpt[16]);
927943

928944
printf("\n-------------------------------------------------------------------"
929945
"----------------------------------------------\n");
@@ -951,7 +967,7 @@ void FitMultipleGrains(struct GrainData *grains, int nGrains, double OptP[13],
951967
Out[i] = x[i];
952968

953969
if (nPanels > 1) {
954-
int p_idx = global_offset + 13;
970+
int p_idx = global_offset + 17;
955971
for (int i = 0; i < nPanels; i++) {
956972
if (i == fixPanel) {
957973
panels[i].dY = 0;
@@ -1016,6 +1032,7 @@ int main(int argc, char *argv[]) {
10161032
double Lsd = cfg.Lsd, px = cfg.px, Wavelength = cfg.Wavelength;
10171033
double p0 = cfg.p0, p1 = cfg.p1, p2 = cfg.p2, p3 = cfg.p3;
10181034
double p4 = cfg.p4, p5 = cfg.p5, p6 = cfg.p6;
1035+
double p7 = cfg.p7, p8 = cfg.p8, p9 = cfg.p9, p10 = cfg.p10;
10191036
double RhoD = cfg.RhoD, yBC = cfg.ybc, zBC = cfg.zbc, wedge = cfg.Wedge;
10201037
double MinEta = cfg.MinEta;
10211038
int NrPixels = cfg.NrPixels;
@@ -1310,8 +1327,8 @@ int main(int argc, char *argv[]) {
13101327
double NonOptP[5] = {RhoD, Lsd, px, Wavelength, MinEta};
13111328
// int NonOptPInt[5] = {NrPixels, nOmeRanges, nRings, nSpots, nhkls}; //
13121329
// Removed, unused and caused errors
1313-
double OptP[13] = {tx, ty, tz, yBC, zBC, wedge, p0, p1, p2, p3, p4, p5, p6};
1314-
double tols[25] = {
1330+
double OptP[17] = {tx, ty, tz, yBC, zBC, wedge, p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10};
1331+
double tols[29] = {
13151332
250,
13161333
250,
13171334
250,
@@ -1337,6 +1354,10 @@ int main(int argc, char *argv[]) {
13371354
tolP4,
13381355
tolP4, // tolP5 (reuse tolP4)
13391356
0, // tolP6
1357+
0, // tolP7
1358+
0, // tolP8
1359+
0, // tolP9
1360+
0, // tolP10
13401361
}; // 250 microns for position, 0.0005 degrees for orient, 1 % for
13411362
// latticeParameter, 1 degree for tilts, 1 pixel for yBC,
13421363
// 1 pixel for zBC, 0.00001 degree for wedge, 1E-3 for p0,p1,p2, 45 for

‎FF_HEDM/src/FitSetupParamsAllZarr.c‎

Lines changed: 33 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -152,6 +152,10 @@ struct my_func_data {
152152
double p4;
153153
double p5;
154154
double p6;
155+
double p7;
156+
double p8;
157+
double p9;
158+
double p10;
155159
};
156160

157161
static double problem_function(unsigned n, const double *x, double *grad,
@@ -178,6 +182,10 @@ static double problem_function(unsigned n, const double *x, double *grad,
178182
p4 = f_data->p4;
179183
double p5 = f_data->p5;
180184
double p6 = f_data->p6;
185+
double p7 = f_data->p7;
186+
double p8 = f_data->p8;
187+
double p9 = f_data->p9;
188+
double p10 = f_data->p10;
181189
txr = deg2rad * tx;
182190
tyr = deg2rad * ty;
183191
tzr = deg2rad * tz;
@@ -212,10 +220,13 @@ static double problem_function(unsigned n, const double *x, double *grad,
212220
Eta = CalcEtaAngleLocal(XYZ[1], XYZ[2]);
213221
RNorm = Rad / MaxRad;
214222
EtaT = 90 - Eta;
223+
double RNorm3 = RNorm * RNorm * RNorm;
224+
double dipole = p7 * RNorm * cos(deg2rad * (EtaT + p8));
225+
double trefoil = p9 * RNorm3 * cos(deg2rad * (3.0 * EtaT + p10));
215226
DistortFunc = (p0 * (pow(RNorm, n0)) * (cos(deg2rad * (2 * EtaT + p6)))) +
216227
(p1 * (pow(RNorm, n1)) * (cos(deg2rad * (4 * EtaT + p3)))) +
217228
(p2 * (pow(RNorm, n2))) + p4 * pow(RNorm, 6.0) +
218-
p5 * pow(RNorm, 4.0) + 1;
229+
p5 * pow(RNorm, 4.0) + dipole + trefoil + 1;
219230
Rcorr = Rad * DistortFunc;
220231
RIdeal = Lsd * tan(deg2rad * IdealTtheta[i]);
221232
Diff = fabs(1 - (Rcorr / RIdeal));
@@ -248,7 +259,7 @@ void FitTiltBCLsd(int nIndices, double *YMean, double *ZMean,
248259
double zbc, double tx, double tyIn, double tzIn, double *ty,
249260
double *tz, double *LsdFit, double *ybcFit, double *zbcFit,
250261
double p0, double p1, double p2, double p3, double p4, double p5,
251-
double p6, double *MeanDiff, double tolTilts, double tolLsd,
262+
double p6, double p7, double p8, double p9, double p10, double *MeanDiff, double tolTilts, double tolLsd,
252263
double tolBC, double px) {
253264
unsigned n = 5;
254265
struct my_func_data f_data;
@@ -266,6 +277,10 @@ void FitTiltBCLsd(int nIndices, double *YMean, double *ZMean,
266277
f_data.p4 = p4;
267278
f_data.p5 = p5;
268279
f_data.p6 = p6;
280+
f_data.p7 = p7;
281+
f_data.p8 = p8;
282+
f_data.p9 = p9;
283+
f_data.p10 = p10;
269284
double x[n], xl[n], xu[n];
270285
x[0] = Lsd;
271286
xl[0] = Lsd - tolLsd;
@@ -311,7 +326,7 @@ void FitTiltBCLsd(int nIndices, double *YMean, double *ZMean,
311326
static inline void CorrectTiltSpatialDistortion(
312327
int nIndices, double MaxRad, double *YMean, double *ZMean, double px,
313328
double Lsd, double ybc, double zbc, double tx, double ty, double tz,
314-
double p0, double p1, double p2, double p3, double p4, double p5, double p6,
329+
double p0, double p1, double p2, double p3, double p4, double p5, double p6, double p7, double p8, double p9, double p10,
315330
double *YCorrected, double *ZCorrected, int nPanels, Panel *panels) {
316331
double txr, tyr, tzr;
317332
txr = deg2rad * tx;
@@ -597,7 +612,7 @@ int main(int argc, char *argv[]) {
597612
double LatticeConstant[6], Wavelength, MaxRingRad, Lsd, MaxTtheta, TthetaTol,
598613
ybc, zbc, px, tyIn, tzIn, BeamSize = 0;
599614
double tx, tolTilts = 1, tolLsd = 5000, tolBC = 1, p0, p1, p2, p3, p4 = 0, p5 = 0,
600-
p6 = 0, RhoD, wedge, MinEta;
615+
p6 = 0, p7 = 0, p8 = 0, p9 = 0, p10 = 0, RhoD, wedge, MinEta;
601616
// Panel parameters
602617
int NPanelsY = 0, NPanelsZ = 0, PanelSizeY = 0, PanelSizeZ = 0;
603618
int locPanelGapsY = -1, locPanelGapsZ = -1, nGapsY = 0, nGapsZ = 0;
@@ -944,6 +959,18 @@ int main(int argc, char *argv[]) {
944959
if (strstr(finfo->name, "analysis/process/analysis_parameters/p6/0") !=
945960
NULL)
946961
ReadZarrChunk(arch, count, &p6, sizeof(double));
962+
if (strstr(finfo->name, "analysis/process/analysis_parameters/p7/0") !=
963+
NULL)
964+
ReadZarrChunk(arch, count, &p7, sizeof(double));
965+
if (strstr(finfo->name, "analysis/process/analysis_parameters/p8/0") !=
966+
NULL)
967+
ReadZarrChunk(arch, count, &p8, sizeof(double));
968+
if (strstr(finfo->name, "analysis/process/analysis_parameters/p9/0") !=
969+
NULL)
970+
ReadZarrChunk(arch, count, &p9, sizeof(double));
971+
if (strstr(finfo->name, "analysis/process/analysis_parameters/p10/0") !=
972+
NULL)
973+
ReadZarrChunk(arch, count, &p10, sizeof(double));
947974
if (strstr(finfo->name,
948975
"analysis/process/analysis_parameters/NPanelsY/0") != NULL) {
949976
ReadZarrChunk(arch, count, &NPanelsY, sizeof(int));
@@ -1294,7 +1321,7 @@ int main(int argc, char *argv[]) {
12941321
printf("Fitting parameters.\n");
12951322
FitTiltBCLsd(nIndices, Ys, Zs, IdealTtheta, Lsd, RhoD, ybc, zbc, tx, tyIn,
12961323
tzIn, &ty, &tz, &LsdFit, &ybcFit, &zbcFit, p0, p1, p2, p3, p4, p5,
1297-
p6, &MeanDiff, tolTilts, tolLsd, tolBC, px);
1324+
p6, p7, p8, p9, p10, &MeanDiff, tolTilts, tolLsd, tolBC, px);
12981325
printf("Number of function calls: %d\n", NrCalls);
12991326
printf("LsdFit:\t\t%0.12f\nYBCFit:\t\t%0.12f\nZBCFit:\t\t%0.12f\ntyFit:"
13001327
"\t\t%0.12f\ntzFit:\t\t%0.12f\nMeanStrain:\t%0.12lf\n",
@@ -1314,7 +1341,7 @@ int main(int argc, char *argv[]) {
13141341
YCorrected = malloc(nIndices * sizeof(*YCorrected));
13151342
ZCorrected = malloc(nIndices * sizeof(*ZCorrected));
13161343
CorrectTiltSpatialDistortion(nIndices, RhoD, Ys, Zs, px, LsdFit, ybcFit,
1317-
zbcFit, tx, ty, tz, p0, p1, p2, p3, p4, p5, p6,
1344+
zbcFit, tx, ty, tz, p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10,
13181345
YCorrected, ZCorrected, nPanels, panels);
13191346
double *YCorrWedge, *ZCorrWedge, *OmegaCorrWedge, *EtaCorrWedge,
13201347
*TthetaCorrWedge, YCorrWedgeT, ZCorrWedgeT, OmegaCorrWedgeT,

‎FF_HEDM/src/FitWedgeParallel.c‎

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -90,7 +90,8 @@ struct my_func_data {
9090
};
9191

9292
void YsZsCalc(double Lsd, double Ycen, double Zcen, double p0, double p1,
93-
double p2, double p3, double p4, double p5, double p6, double MaxRad, double tx,
93+
double p2, double p3, double p4, double p5, double p6,
94+
double p7, double p8, double p9, double p10, double MaxRad, double tx,
9495
double ty, double tz, double px, double Ys, double Zs,
9596
double *YOut, double *ZOut, int nPanels, Panel *panels) {
9697
double txr = deg2rad * tx;
@@ -130,10 +131,13 @@ void YsZsCalc(double Lsd, double Ycen, double Zcen, double p0, double p1,
130131
Eta = CalcEtaAngleLocal(XYZ[1], XYZ[2]);
131132
RNorm = Rad / MaxRad;
132133
EtaT = 90 - Eta;
134+
double RNorm3 = RNorm * RNorm * RNorm;
135+
double dipole = p7 * RNorm * cos(deg2rad * (EtaT + p8));
136+
double trefoil = p9 * RNorm3 * cos(deg2rad * (3.0 * EtaT + p10));
133137
DistortFunc = (p0 * (pow(RNorm, n0)) * (cos(deg2rad * (2 * EtaT + p6)))) +
134138
(p1 * (pow(RNorm, n1)) * (cos(deg2rad * (4 * EtaT + p3)))) +
135139
(panelP2 * (pow(RNorm, n2))) + p4 * pow(RNorm, 6.0) +
136-
p5 * pow(RNorm, 4.0) + 1;
140+
p5 * pow(RNorm, 4.0) + dipole + trefoil + 1;
137141
Rcorr = Rad * DistortFunc;
138142
Rcorr = Rcorr * (Lsd / panelLsd); // re-project to global Lsd plane
139143
*YOut = -Rcorr * sin(deg2rad * Eta);
@@ -268,7 +272,7 @@ static double problem_function(unsigned n, const double *x, double *grad,
268272
// Returns 0 if success, 1 if hit bounds
269273
int FitWedgeThreadSafe(double Lsd, double Ycen, double Zcen, double p0,
270274
double p1, double p2, double p3, double p4, double p5,
271-
double p6, double MaxRad, double tx, double ty, double tz,
275+
double p6, double p7, double p8, double p9, double p10, double MaxRad, double tx, double ty, double tz,
272276
double px, double Ys, double Zs, double MinOme,
273277
double MaxOme, double WedgeIn, double *WedgeFit,
274278
double *MinFOut, double *NumUncertainty,
@@ -279,7 +283,7 @@ int FitWedgeThreadSafe(double Lsd, double Ycen, double Zcen, double p0,
279283
f_data.Ome2 = MaxOme;
280284
f_data.Wavelength = Wavelength;
281285
double Y, Z;
282-
YsZsCalc(Lsd, Ycen, Zcen, p0, p1, p2, p3, p4, p5, p6, MaxRad, tx, ty, tz, px, Ys, Zs,
286+
YsZsCalc(Lsd, Ycen, Zcen, p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, MaxRad, tx, ty, tz, px, Ys, Zs,
283287
&Y, &Z, nPanels, panels);
284288
f_data.Ys = Y;
285289
f_data.Zs = Z;
@@ -382,6 +386,7 @@ int main(int argc, char *argv[]) {
382386
double Wedge = cfg.Wedge;
383387
double p0 = cfg.p0, p1 = cfg.p1, p2 = cfg.p2, p3 = cfg.p3;
384388
double p4 = cfg.p4, p5 = cfg.p5, p6 = cfg.p6;
389+
double p7 = cfg.p7, p8 = cfg.p8, p9 = cfg.p9, p10 = cfg.p10;
385390
double MaxRad = cfg.RhoD;
386391
double Wavelength = cfg.Wavelength;
387392
double OmegaStart = cfg.OmegaStart;
@@ -578,7 +583,7 @@ int main(int argc, char *argv[]) {
578583
double resMinF = 0;
579584
double resNumUncertainty = 0;
580585
int status = FitWedgeThreadSafe(
581-
Lsd, Ycen, Zcen, p0, p1, p2, p3, p4, p5, p6, MaxRad, tx, ty, tz, px, YsUse,
586+
Lsd, Ycen, Zcen, p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, MaxRad, tx, ty, tz, px, YsUse,
582587
ZsUse, MinOme, MaxOme, Wedge, &resWedge, &resMinF, &resNumUncertainty,
583588
Wavelength, nPanels, panels);
584589

‎FF_HEDM/src/ForwardSimulationCompressed.c‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -600,7 +600,7 @@ static inline double CorrectWedge(double eta, double theta, double wl,
600600
static inline void CorrectTiltSpatialDistortion(
601601
double px, double Lsd, double ybc, double zbc, double tx, double ty,
602602
double tz, double RhoD, double p0, double p1, double p2, double p3,
603-
double p4, double p5, double p6, int NrPixels, double *yDispl, double *zDispl,
603+
double p4, double p5, double p6, double p7, double p8, double p9, double p10, int NrPixels, double *yDispl, double *zDispl,
604604
int nPanels, Panel *panels, int nCPUs) {
605605
double txr, tyr, tzr;
606606
txr = deg2rad * tx;
@@ -650,11 +650,14 @@ static inline void CorrectTiltSpatialDistortion(
650650
CalcEtaAngle(XYZ[1], XYZ[2], &Eta);
651651
double RNorm = Rad / RhoD;
652652
double EtaT = 90 - Eta;
653+
double RNorm3 = RNorm * RNorm * RNorm;
654+
double dipole = p7 * RNorm * cos(deg2rad * (EtaT + p8));
655+
double trefoil = p9 * RNorm3 * cos(deg2rad * (3.0 * EtaT + p10));
653656
double DistortFunc =
654657
(p0 * (pow(RNorm, n0)) * (cos(deg2rad * (2 * EtaT + p6)))) +
655658
(p1 * (pow(RNorm, n1)) * (cos(deg2rad * (4 * EtaT + p3)))) +
656659
(panelP2 * (pow(RNorm, n2))) + p4 * pow(RNorm, 6.0) +
657-
p5 * pow(RNorm, 4.0) + 1;
660+
p5 * pow(RNorm, 4.0) + dipole + trefoil + 1;
658661
double Rcorr = Rad * DistortFunc;
659662
double YCorr = -Rcorr * sin(Eta * deg2rad);
660663
double ZCorr = Rcorr * cos(Eta * deg2rad);
@@ -1094,6 +1097,7 @@ int main(int argc, char *argv[]) {
10941097
memcpy(LatC, cfg.LatticeConstant, sizeof(LatC));
10951098
double p0 = cfg.p0, p1 = cfg.p1, p2 = cfg.p2, p3 = cfg.p3;
10961099
double p4 = cfg.p4, p5 = cfg.p5, p6 = cfg.p6, RhoD = cfg.RhoD;
1100+
double p7 = cfg.p7, p8 = cfg.p8, p9 = cfg.p9, p10 = cfg.p10;
10971101
double GaussWidth = cfg.GaussWidth > 0 ? cfg.GaussWidth : 1;
10981102
double PeakIntensity = cfg.PeakIntensity;
10991103
int NPanelsY = cfg.NPanelsY, NPanelsZ = cfg.NPanelsZ;
@@ -1443,7 +1447,7 @@ int main(int argc, char *argv[]) {
14431447
-32100.0; // Set to a special number to check if it was set or not.
14441448
}
14451449
CorrectTiltSpatialDistortion(px, Lsd, yBC, zBC, tx, ty, tz, RhoD, p0, p1, p2,
1446-
p3, p4, p5, p6, NrPixels, yDispl, zDispl, nPanels,
1450+
p3, p4, p5, p6, p7, p8, p9, p10, NrPixels, yDispl, zDispl, nPanels,
14471451
panels, nCPUs);
14481452
end_time = omp_get_wtime();
14491453
diftotal = end_time - start0;

0 commit comments

Comments
 (0)