Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
ptwXY_functions.c File Reference
#include <math.h>
#include "ptwXY.h"

Go to the source code of this file.

Functions

nfu_status ptwXY_pow (statusMessageReporting *smr, ptwXYPoints *ptwXY, double v)
nfu_status ptwXY_exp (statusMessageReporting *smr, ptwXYPoints *ptwXY, double a)
ptwXYPoints * ptwXY_convolution (statusMessageReporting *smr, ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, int mode)
ptwXYPoints * ptwXY_inverse (statusMessageReporting *smr, ptwXYPoints *ptwXY)

Function Documentation

◆ ptwXY_convolution()

ptwXYPoints * ptwXY_convolution ( statusMessageReporting * smr,
ptwXYPoints * ptwXY1,
ptwXYPoints * ptwXY2,
int mode )

Definition at line 111 of file ptwXY_functions.c.

111 {
112/*
113* Currently, only supports linear-linear interpolation.
114*
115* This function calculates c(y) = integral dx f1(x) * f2(y-x)
116*
117*/
118 int64_t i1, i2, n1, n2, n;
119 ptwXYPoints *f1 = ptwXY1, *f2 = ptwXY2, *convolute;
120 double accuracy = ptwXY1->accuracy, rangeMin, rangeMax, c, y, dy;
121
122 if( ptwXY_simpleCoalescePoints( smr, ptwXY1 ) != nfu_Okay ) {
123 smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_Error, "Via (source1)." );
124 return( NULL );
125 }
126 if( ptwXY_simpleCoalescePoints( smr, ptwXY2 ) != nfu_Okay ) {
127 smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_Error, "Via (source2)." );
128 return( NULL );
129 }
130
131 if( ptwXY1->interpolation != ptwXY_interpolationLinLin ) {
132 smr_setReportError2( smr, nfu_SMR_libraryID, nfu_unsupportedInterpolation,
133 "Source1: unsupported interpolation = '%s'", ptwXY1->interpolationString );
134 return( NULL );
135 }
136 if( ptwXY2->interpolation != ptwXY_interpolationLinLin ) {
137 smr_setReportError2( smr, nfu_SMR_libraryID, nfu_unsupportedInterpolation,
138 "Source2: unsupported interpolation = '%s'", ptwXY2->interpolationString );
139 return( NULL );
140 }
141
142 n1 = f1->length;
143 n2 = f2->length;
144
145 if( ( n1 == 0 ) || ( n2 == 0 ) ) {
146 convolute = ptwXY_new( smr, ptwXY_interpolationLinLin, NULL, 1., accuracy, 0, 0, 0 );
147 if( convolute == NULL ) smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_Error, "Via." );
148 return( convolute );
149 }
150
151 if( ( n1 == 1 ) || ( n2 == 1 ) ) {
152 smr_setReportError2( smr, nfu_SMR_libraryID, nfu_tooFewPoints,
153 "Too few points: len( source1 ) = %d, len( source1 ) = %d.", (int) n1, (int) n2 );
154 return( NULL );
155 }
156
157 if( accuracy < ptwXY2->accuracy ) accuracy = ptwXY2->accuracy;
158 n = n1 * n2;
159 if( mode == 0 ) {
160 mode = 1;
161 if( n > 10000 ) mode = -1;
162 }
163 if( n > 100000 ) mode = -1;
164 if( ( convolute = ptwXY_new( smr, ptwXY_interpolationLinLin, NULL, 1., accuracy, 400, 40, 0 ) ) == NULL ) {
165 smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_Error, "Via." );
166 return( NULL );
167 }
168
169 rangeMin = f1->points[0].x + f2->points[0].x;
170 rangeMax = f1->points[n1 - 1].x + f2->points[n2 - 1].x;
171
172 if( ptwXY_setValueAtX( smr, convolute, rangeMin, 0. ) != nfu_Okay ) goto Err;
173
174 if( mode < 0 ) {
175 dy = ( rangeMax - rangeMin ) / 2000;
176 for( y = rangeMin + dy; y < rangeMax; y += dy ) {
177 if( ptwXY_convolution2( smr, f1, f2, y, rangeMin, &c ) != nfu_Okay ) goto Err;
178 if( ptwXY_setValueAtX( smr, convolute, y, c ) != nfu_Okay ) goto Err;
179 } }
180 else {
181 for( i1 = 0; i1 < n1; i1++ ) {
182 for( i2 = 0; i2 < n2; i2++ ) {
183 y = rangeMin + ( f1->points[i1].x - f1->points[0].x ) + ( f2->points[i2].x - f2->points[0].x );
184 if( y <= rangeMin ) continue;
185 if( y >= rangeMax ) continue;
186 if( ptwXY_convolution2( smr, f1, f2, y, rangeMin, &c ) != nfu_Okay ) goto Err;
187 if( ptwXY_setValueAtX( smr, convolute, y, c ) != nfu_Okay ) goto Err;
188 }
189 }
190 }
191 if( ptwXY_setValueAtX( smr, convolute, rangeMax, 0. ) != nfu_Okay ) goto Err;
192 if( ptwXY_simpleCoalescePoints( smr, convolute ) != nfu_Okay ) goto Err;
193 for( i1 = convolute->length - 1; i1 > 0; i1-- ) {
194 if( ptwXY_convolution3( smr, convolute, f1, f2, convolute->points[i1 - 1].x, convolute->points[i1 - 1].y,
195 convolute->points[i1].x, convolute->points[i1].y, rangeMin ) != nfu_Okay ) goto Err;
196 }
197
198 return( convolute );
199
200Err:
201 ptwXY_free( convolute );
202 smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_Error, "Via." );
203 return( NULL );
204}
nfu_unsupportedInterpolation
@ nfu_unsupportedInterpolation
Definition nf_utilities.h:51
nfu_Okay
@ nfu_Okay
Definition nf_utilities.h:39
nfu_tooFewPoints
@ nfu_tooFewPoints
Definition nf_utilities.h:53
nfu_Error
@ nfu_Error
Definition nf_utilities.h:40
nfu_SMR_libraryID
int nfu_SMR_libraryID
Definition nf_utilities.c:46
ptwXY_setValueAtX
nfu_status ptwXY_setValueAtX(statusMessageReporting *smr, ptwXYPoints *ptwXY, double x, double y)
Definition ptwXY_core.c:1225
ptwXY_interpolationLinLin
@ ptwXY_interpolationLinLin
Definition ptwXY.h:37
ptwXYPoints
struct ptwXYPoints_s ptwXYPoints
ptwXY_new
ptwXYPoints * ptwXY_new(statusMessageReporting *smr, ptwXY_interpolation interpolation, char const *interpolationString, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, int userFlag)
Definition ptwXY_core.c:28
ptwXY_free
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition ptwXY_core.c:782
ptwXY_simpleCoalescePoints
nfu_status ptwXY_simpleCoalescePoints(statusMessageReporting *smr, ptwXYPoints *ptwXY)
Definition ptwXY_core.c:734
smr_setReportError2
#define smr_setReportError2(smr, libraryID, code, fmt,...)
Definition statusMessageReporting.h:128
smr_setReportError2p
#define smr_setReportError2p(smr, libraryID, code, fmt)
Definition statusMessageReporting.h:129
ptwXYPoint_s::x
double x
Definition ptwXY.h:64
ptwXYPoints_s::points
ptwXYPoint * points
Definition ptwXY.h:93
ptwXYPoints_s::interpolation
ptwXY_interpolation interpolation
Definition ptwXY.h:81
ptwXYPoints_s::accuracy
double accuracy
Definition ptwXY.h:85
ptwXYPoints_s::length
int64_t length
Definition ptwXY.h:87
ptwXYPoints_s::interpolationString
char const * interpolationString
Definition ptwXY.h:82

◆ ptwXY_exp()

nfu_status ptwXY_exp ( statusMessageReporting * smr,
ptwXYPoints * ptwXY,
double a )

Definition at line 46 of file ptwXY_functions.c.

46 {
47
48 int64_t i, length;
49 nfu_status status;
50 double x1, y1, z1, x2, y2, z2;
51
52 length = ptwXY->length;
53 if( length < 1 ) return( ptwXY->status );
54
55 if( ptwXY->interpolation == ptwXY_interpolationOther ) {
56 smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_otherInterpolation, "Other interpolation not allowed." );
57 return( ptwXY->status = nfu_otherInterpolation );
58 }
59 if( ptwXY->interpolation == ptwXY_interpolationFlat ) {
60 smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_flatInterpolation, "Flat interpolation not allowed." );
61 return( ptwXY->status = nfu_flatInterpolation );
62 }
63
64 if( ( status = ptwXY_simpleCoalescePoints( smr, ptwXY ) ) != nfu_Okay ) goto Err;
65 x2 = ptwXY->points[length-1].x;
66 y2 = a * ptwXY->points[length-1].y;
67 z2 = ptwXY->points[length-1].y = exp( y2 );
68 for( i = length - 2; i >= 0; i-- ) {
69 x1 = ptwXY->points[i].x;
70 y1 = a * ptwXY->points[i].y;
71 z1 = ptwXY->points[i].y = exp( y1 );
72 if( ( status = ptwXY_exp_s( smr, ptwXY, x1, y1, z1, x2, y2, z2, 0 ) ) != nfu_Okay ) goto Err;
73 x2 = x1;
74 y2 = y1;
75 }
76 return( status );
77
78Err:
79 if( status != nfu_Okay ) smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_Error, "Via." );
80 if( ptwXY->status != nfu_Okay ) ptwXY->status = status;
81 return( status );
82}
nfu_flatInterpolation
@ nfu_flatInterpolation
Definition nf_utilities.h:59
nfu_otherInterpolation
@ nfu_otherInterpolation
Definition nf_utilities.h:58
nfu_status
enum nfu_status_e nfu_status
ptwXY_interpolationFlat
@ ptwXY_interpolationFlat
Definition ptwXY.h:38
ptwXY_interpolationOther
@ ptwXY_interpolationOther
Definition ptwXY.h:38
ptwXYPoint_s::y
double y
Definition ptwXY.h:64
ptwXYPoints_s::status
nfu_status status
Definition ptwXY.h:80

◆ ptwXY_inverse()

ptwXYPoints * ptwXY_inverse ( statusMessageReporting * smr,
ptwXYPoints * ptwXY )

Definition at line 356 of file ptwXY_functions.c.

356 {
357
358 int64_t length;
359 ptwXY_interpolation interpolation;
360 ptwXYPoints *ptwXYInverse;
361
362 length = ptwXY_length( NULL, ptwXY );
363
364 if( ptwXY->interpolation == ptwXY_interpolationFlat ) {
365 smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_flatInterpolation, "flat interpolation not allowed." );
366 return( NULL );
367 }
368
369 if( ptwXY->interpolation == ptwXY_interpolationOther ) {
370 smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_otherInterpolation, "Other interpolation not allowed." );
371 return( NULL );
372 }
373
374 if( ptwXY_simpleCoalescePoints( smr, ptwXY ) != nfu_Okay ) {
375 smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_Error, "Via." );
376 return( NULL );
377 }
378
379 switch( ptwXY->interpolation ) {
380 case ptwXY_interpolationLogLin :
381 interpolation = ptwXY_interpolationLinLog;
382 break;
383 case ptwXY_interpolationLinLog :
384 interpolation = ptwXY_interpolationLogLin;
385 break;
386 default :
387 interpolation = ptwXY->interpolation;
388 break;
389 }
390
391 if( ( ptwXYInverse = ptwXY_new( smr, interpolation, NULL, ptwXY_getBiSectionMax( ptwXY ), ptwXY_getAccuracy( ptwXY ),
392 length, 10, 0 ) ) == NULL ) {
393 smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_Error, "Via." );
394 return( NULL );
395 }
396
397 if( length == 1 ) {
398 ptwXYInverse->points[0].x = ptwXY->points[0].y;
399 ptwXYInverse->points[0].y = ptwXY->points[0].x; }
400 else if( length > 1 ) {
401 int64_t i1, start = 0, order = 1;
402
403 if( ptwXY->points[0].y > ptwXY->points[1].y ) {
404 start = length - 1;
405 order = -1;
406 }
407 ptwXYInverse->points[0].x = ptwXY->points[start].y;
408 ptwXYInverse->points[0].y = ptwXY->points[start].x;
409 for( i1 = 1, start += order; i1 < length; ++i1, start += order ) {
410 ptwXYInverse->points[i1].x = ptwXY->points[start].y;
411 ptwXYInverse->points[i1].y = ptwXY->points[start].x;
412 if( ptwXYInverse->points[i1-1].x >= ptwXYInverse->points[i1].x ) {
413 smr_setReportError2( smr, nfu_SMR_libraryID, nfu_XNotAscending,
414 "Non-ascending domain values: x[%d] = %.17e >= x[%d] = %.17e.",
415 (int) (i1-1), ptwXYInverse->points[i1-1].x, (int) i1, ptwXYInverse->points[i1].x );
416 ptwXY_free( ptwXYInverse );
417 return( NULL );
418 }
419 }
420 }
421 ptwXYInverse->length = length;
422 return( ptwXYInverse );
423}
nfu_XNotAscending
@ nfu_XNotAscending
Definition nf_utilities.h:45
ptwXY_getAccuracy
double ptwXY_getAccuracy(ptwXYPoints *ptwXY)
Definition ptwXY_core.c:566
ptwXY_interpolation
enum ptwXY_interpolation_e ptwXY_interpolation
ptwXY_interpolationLinLog
@ ptwXY_interpolationLinLog
Definition ptwXY.h:37
ptwXY_interpolationLogLin
@ ptwXY_interpolationLogLin
Definition ptwXY.h:37
ptwXY_getBiSectionMax
double ptwXY_getBiSectionMax(ptwXYPoints *ptwXY)
Definition ptwXY_core.c:582
ptwXY_length
int64_t ptwXY_length(statusMessageReporting *smr, ptwXYPoints *ptwXY)
Definition ptwXY_core.c:793

◆ ptwXY_pow()

nfu_status ptwXY_pow ( statusMessageReporting * smr,
ptwXYPoints * ptwXY,
double v )

Definition at line 24 of file ptwXY_functions.c.

24 {
25
26 nfu_status status = ptwXY_applyFunction( smr, ptwXY, ptwXY_pow_callback, (void *) &v, 0 );
27
28 if( status != nfu_Okay ) smr_setReportError2p( smr, nfu_SMR_libraryID, nfu_Error, "Via." );
29 return( status );
30}
ptwXY_applyFunction
nfu_status ptwXY_applyFunction(statusMessageReporting *smr, ptwXYPoints *ptwXY1, ptwXY_applyFunction_callback func, void *argList, int checkForRoots)
Definition ptwXY_misc.c:165