Geant4  10.02.p02
ptwXY_convenient.cc
Go to the documentation of this file.
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5 
6 #include <stdlib.h>
7 #include <cmath>
8 #include <float.h>
9 
10 #include "ptwXY.h"
11 
12 #if defined __cplusplus
13 #include <cmath>
14 #include "G4Exp.hh"
15 #include "G4Log.hh"
16 namespace GIDI {
17 using namespace GIDI;
18 #endif
19 
20 static nfu_status ptwXY_createGaussianCenteredSigma1_2( ptwXYPoints *ptwXY, double x1, double y1, double x2, double y2, int addX1Point );
21 /*
22 ************************************************************
23 */
24 ptwXPoints *ptwXY_getXArray( ptwXYPoints *ptwXY, nfu_status *status ) {
25 
26  int64_t i, n;
27  ptwXPoints *xArray;
28 
29  if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
30  n = ptwXY->length;
31 
32  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
33  if( ( xArray = ptwX_new( n, status ) ) == NULL ) return( NULL );
34  for( i = 0; i < n; i++ ) xArray->points[i] = ptwXY->points[i].x;
35  xArray->length = n;
36 
37  return( xArray );
38 }
39 /*
40 ************************************************************
41 */
42 nfu_status ptwXY_dullEdges( ptwXYPoints *ptwXY, double lowerEps, double upperEps, int positiveXOnly ) {
43 
44 #define minEps 5e-16
45 
46  nfu_status status;
47  double xm, xp, dx, y, x1, y1, x2, y2, sign;
48  ptwXYPoint *p;
49 
50 /* This routine can only be used for linear interpolation for the y-axes since for log interpolation, y cannot be 0.
51 This needs to be fixed and documented. */
52  if( ( status = ptwXY->status ) != nfu_Okay ) return( status );
53  if( ptwXY->interpolation == ptwXY_interpolationFlat ) return( nfu_invalidInterpolation );
54  if( ptwXY->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
55 
56  if( ptwXY->length < 2 ) return( nfu_Okay );
57 
58  if( lowerEps != 0. ) {
59  if( std::fabs( lowerEps ) < minEps ) {
60  sign = 1;
61  if( lowerEps < 0. ) sign = -1;
62  lowerEps = sign * minEps;
63  }
64 
65  p = ptwXY_getPointAtIndex_Unsafely( ptwXY, 0 );
66  x1 = p->x;
67  y1 = p->y;
68  p = ptwXY_getPointAtIndex_Unsafely( ptwXY, 1 );
69  x2 = p->x;
70  y2 = p->y;
71 
72  if( y1 != 0. ) {
73  dx = std::fabs( x1 * lowerEps );
74  if( x1 == 0 ) dx = std::fabs( lowerEps );
75  xm = x1 - dx;
76  xp = x1 + dx;
77  if( ( xp + dx ) < x2 ) {
78  if( ( status = ptwXY_getValueAtX( ptwXY, xp, &y ) ) != nfu_Okay ) return( status );
79  if( ( status = ptwXY_setValueAtX( ptwXY, xp, y ) ) != nfu_Okay ) return( status ); }
80  else {
81  xp = x2;
82  y = y2;
83  }
84  if( lowerEps > 0 ) {
85  if( ( status = ptwXY_setValueAtX( ptwXY, x1, 0. ) ) != nfu_Okay ) return( status ); }
86  else {
87  if( ( xm < 0. ) && ( x1 >= 0. ) && positiveXOnly ) {
88  if( ( status = ptwXY_setValueAtX( ptwXY, x1, 0. ) ) != nfu_Okay ) return( status ); }
89  else {
90  if( ( status = ptwXY_setValueAtX( ptwXY, xm, 0. ) ) != nfu_Okay ) return( status );
91  if( ( status = ptwXY_interpolatePoint( ptwXY->interpolation, x1, &y, xm, 0., xp, y ) ) != nfu_Okay ) return( status );
92  if( ( status = ptwXY_setValueAtX( ptwXY, x1, y ) ) != nfu_Okay ) return( status );
93  }
94  }
95  }
96  }
97 
98  if( upperEps != 0. ) {
99  if( std::fabs( upperEps ) < minEps ) {
100  sign = 1;
101  if( upperEps < 0. ) sign = -1;
102  upperEps = sign * minEps;
103  }
104 
105  p = ptwXY_getPointAtIndex_Unsafely( ptwXY, ptwXY->length - 2 );
106  x1 = p->x;
107  y1 = p->y;
108  p = ptwXY_getPointAtIndex_Unsafely( ptwXY, ptwXY->length - 1 );
109  x2 = p->x;
110  y2 = p->y;
111 
112  if( y2 != 0. ) {
113  dx = std::fabs( x2 * upperEps );
114  if( x2 == 0 ) dx = std::fabs( upperEps );
115  xm = x2 - dx;
116  xp = x2 + dx;
117  if( ( xm - dx ) > x1 ) {
118  if( ( status = ptwXY_getValueAtX( ptwXY, xm, &y ) ) != nfu_Okay ) return( status );
119  if( ( status = ptwXY_setValueAtX( ptwXY, xm, y ) ) != nfu_Okay ) return( status ); }
120  else {
121  xm = x1;
122  y = y1;
123  }
124  if( upperEps < 0 ) {
125  if( ( status = ptwXY_setValueAtX( ptwXY, x2, 0. ) ) != nfu_Okay ) return( status ); }
126  else {
127  if( ( status = ptwXY_setValueAtX( ptwXY, xp, 0. ) ) != nfu_Okay ) return( status );
128  if( ( status = ptwXY_interpolatePoint( ptwXY->interpolation, x2, &y, xm, y, xp, 0. ) ) != nfu_Okay ) return( status );
129  if( ( status = ptwXY_setValueAtX( ptwXY, x2, y ) ) != nfu_Okay ) return( status );
130  }
131  }
132  }
133 
134  return( ptwXY->status );
135 
136 #undef minEps
137 }
138 /*
139 ************************************************************
140 */
141 nfu_status ptwXY_mergeClosePoints( ptwXYPoints *ptwXY, double epsilon ) {
142 
143  int64_t i, i1, j, k, n = ptwXY->length;
144  double x, y;
145  ptwXYPoint *p1, *p2;
146 
147  if( n < 2 ) return( ptwXY->status );
148  if( epsilon < 4 * DBL_EPSILON ) epsilon = 4 * DBL_EPSILON;
149  if( ptwXY_simpleCoalescePoints( ptwXY ) != nfu_Okay ) return( ptwXY->status );
150 
151  p2 = ptwXY->points;
152  x = p2->x;
153  for( i1 = 1, p2++; i1 < ( n - 1 ); i1++, p2++ ) { /* The first point shall remain the first point and all points close to it are deleted. */
154  if( ( p2->x - x ) > 0.5 * epsilon * ( std::fabs( x ) + std::fabs( p2->x ) ) ) break;
155  }
156  if( i1 != 1 ) {
157  for( i = i1, p1 = &(ptwXY->points[1]); i < n; i++, p1++, p2++ ) *p1 = *p2;
158  n = ptwXY->length = ptwXY->length - i1 + 1;
159  }
160 
161  p1 = &(ptwXY->points[n-1]);
162  x = p1->x;
163  for( i1 = n - 2, p1--; i1 > 0; i1--, p1-- ) { /* The last point shall remain the last point and all points close to it are deleted. */
164  if( x - p1->x > 0.5 * epsilon * ( std::fabs( x ) + std::fabs( p1->x ) ) ) break;
165  }
166  if( i1 != ( n - 2 ) ) {
167  ptwXY->points[i1 + 1] = ptwXY->points[n - 1];
168  n = ptwXY->length = i1 + 2;
169  }
170 
171  for( i = 1; i < n - 1; i++ ) {
172  p1 = &(ptwXY->points[i]);
173  x = p1->x;
174  y = p1->y;
175  for( j = i + 1, p2 = &(ptwXY->points[i+1]); j < n - 1; j++, p2++ ) {
176  if( ( p2->x - p1->x ) > 0.5 * epsilon * ( std::fabs( p2->x ) + std::fabs( p1->x ) ) ) break;
177  x += p2->x;
178  y += p2->y;
179  }
180  if( ( k = ( j - i ) ) > 1 ) {
181  p1->x = x / k;
182  p1->y = y / k;
183  for( p1 = &(ptwXY->points[i+1]); j < n; j++, p1++, p2++ ) *p1 = *p2;
184  n -= ( k - 1 );
185  }
186  }
187  ptwXY->length = n;
188 
189  return( ptwXY->status );
190 }
191 /*
192 ************************************************************
193 */
194 ptwXYPoints *ptwXY_intersectionWith_ptwX( ptwXYPoints *ptwXY, ptwXPoints *ptwX, nfu_status *status ) {
195 
196  int64_t i, i1, i2, lengthX = ptwX_length( ptwX );
197  double x, y, xMin, xMax;
198  ptwXYPoints *n = NULL;
199 
200  if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
201  if( ( *status = ptwX->status ) != nfu_Okay ) return( NULL );
202  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) goto Err;
203  *status = nfu_otherInterpolation;
204  if( ptwXY->interpolation == ptwXY_interpolationOther ) return( NULL );
205  if( ( n = ptwXY_clone( ptwXY, status ) ) == NULL ) return( NULL );
206  if( ptwXY->length == 0 ) return( n );
207  xMin = ptwXY->points[0].x;
208  xMax = ptwXY->points[ptwXY->length - 1].x;
209 
210  if( ( xMin >= ptwX->points[lengthX-1] ) || ( xMax <= ptwX->points[0] ) ) { /* No overlap. */
211  n->length = 0;
212  return( n );
213  }
214 
215  for( i = 0; i < lengthX; i++ ) { /* Fill in ptwXY at x-points in ptwX. */
216  x = ptwX->points[i];
217  if( x <= xMin ) continue;
218  if( x >= xMax ) break;
219  if( ( *status = ptwXY_getValueAtX( ptwXY, x, &y ) ) != nfu_Okay ) goto Err;
220  if( ( *status = ptwXY_setValueAtX( n, x, y ) ) != nfu_Okay ) goto Err;
221  }
222  if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
223 
224  i1 = 0;
225  i2 = n->length - 1;
226  if( lengthX > 0 ) {
227  x = ptwX->points[0];
228  if( x > n->points[i1].x ) {
229  for( ; i1 < n->length; i1++ ) {
230  if( n->points[i1].x == x ) break;
231  }
232  }
233 
234  x = ptwX->points[lengthX - 1];
235  if( x < n->points[i2].x ) {
236  for( ; i2 > i1; i2-- ) {
237  if( n->points[i2].x == x ) break;
238  }
239  }
240  }
241  i2++;
242 
243  if( i1 != 0 ) {
244  for( i = i1; i < i2; i++ ) n->points[i - i1] = n->points[i];
245  }
246  n->length = i2 - i1;
247 
248  return( n );
249 
250 Err:
251  ptwXY_free( n );
252  return( NULL );
253 }
254 /*
255 ************************************************************
256 */
257 nfu_status ptwXY_areDomainsMutual( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2 ) {
258 
259  nfu_status status;
260  int64_t n1 = ptwXY1->length, n2 = ptwXY2->length;
261  ptwXYPoint *xy1, *xy2;
262 
263  if( ( status = ptwXY1->status ) != nfu_Okay ) return( status );
264  if( ( status = ptwXY2->status ) != nfu_Okay ) return( status );
265  if( n1 == 0 ) return( nfu_empty );
266  if( n2 == 0 ) return( nfu_empty );
267  if( n1 < 2 ) {
268  status = nfu_tooFewPoints; }
269  else if( n2 < 2 ) {
270  status = nfu_tooFewPoints; }
271  else {
272  xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, 0 );
273  xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, 0 );
274  if( xy1->x < xy2->x ) {
275  if( xy2->y != 0. ) status = nfu_domainsNotMutual; }
276  else if( xy1->x > xy2->x ) {
277  if( xy1->y != 0. ) status = nfu_domainsNotMutual;
278  }
279 
280  if( status == nfu_Okay ) {
281  xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, n1 - 1 );
282  xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, n2 - 1 );
283  if( xy1->x < xy2->x ) {
284  if( xy1->y != 0. ) status = nfu_domainsNotMutual; }
285  else if( xy1->x > xy2->x ) {
286  if( xy2->y != 0. ) status = nfu_domainsNotMutual;
287  }
288  }
289  }
290  return( status );
291 }
292 /*
293 ************************************************************
294 */
295 nfu_status ptwXY_tweakDomainsToMutualify( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, int epsilonFactor, double epsilon ) {
296 
297  nfu_status status;
298  int64_t n1 = ptwXY1->length, n2 = ptwXY2->length;
299  double sum, diff;
300  ptwXYPoint *xy1, *xy2;
301 
302  epsilon = std::fabs( epsilon ) + std::fabs( epsilonFactor * DBL_EPSILON );
303 
304  if( ( status = ptwXY1->status ) != nfu_Okay ) return( status );
305  if( ( status = ptwXY2->status ) != nfu_Okay ) return( status );
306  if( n1 == 0 ) return( nfu_empty );
307  if( n2 == 0 ) return( nfu_empty );
308  if( n1 < 2 ) {
309  status = nfu_tooFewPoints; }
310  else if( n2 < 2 ) {
311  status = nfu_tooFewPoints; }
312  else {
313  xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, 0 );
314  xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, 0 );
315  if( xy1->x < xy2->x ) {
316  if( xy2->y != 0. ) {
317  sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
318  diff = std::fabs( xy2->x - xy1->x );
319  if( diff > epsilon * sum ) {
320  status = nfu_domainsNotMutual; }
321  else {
322  xy1->x = xy2->x;
323  }
324  } }
325  else if( xy1->x > xy2->x ) {
326  if( xy1->y != 0. ) {
327  sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
328  diff = std::fabs( xy2->x - xy1->x );
329  if( diff > epsilon * sum ) {
330  status = nfu_domainsNotMutual; }
331  else {
332  xy2->x = xy1->x;
333  }
334  }
335  }
336 
337  if( status == nfu_Okay ) {
338  xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, n1 - 1 );
339  xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, n2 - 1 );
340  if( xy1->x < xy2->x ) {
341  if( xy1->y != 0. ) {
342  sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
343  diff = std::fabs( xy2->x - xy1->x );
344  if( diff > epsilon * sum ) {
345  status = nfu_domainsNotMutual; }
346  else {
347  xy2->x = xy1->x;
348  }
349  } }
350  else if( xy1->x > xy2->x ) {
351  if( xy2->y != 0. ) {
352  sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
353  diff = std::fabs( xy2->x - xy1->x );
354  if( diff > epsilon * sum ) {
355  status = nfu_domainsNotMutual; }
356  else {
357  xy1->x = xy2->x;
358  }
359  }
360  }
361  }
362  }
363  return( status );
364 }
365 /*
366 ************************************************************
367 */
368 nfu_status ptwXY_mutualifyDomains( ptwXYPoints *ptwXY1, double lowerEps1, double upperEps1, int positiveXOnly1,
369  ptwXYPoints *ptwXY2, double lowerEps2, double upperEps2, int positiveXOnly2 ) {
370 
371  nfu_status status;
372  int64_t n1 = ptwXY1->length, n2 = ptwXY2->length;
373  ptwXYPoint *xy1, *xy2;
374 
375  switch( status = ptwXY_areDomainsMutual( ptwXY1, ptwXY2 ) ) {
376  case nfu_Okay :
377  case nfu_empty :
378  return( nfu_Okay );
379  case nfu_domainsNotMutual :
380  break;
381  default :
382  return( status );
383  }
384 
385  if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
386  if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
387  if( ptwXY1->interpolation == ptwXY_interpolationFlat ) return( nfu_invalidInterpolation );
388  if( ptwXY2->interpolation == ptwXY_interpolationFlat ) return( nfu_invalidInterpolation );
389 
390  xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, 0 );
391  xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, 0 );
392  if( xy1->x < xy2->x ) {
393  lowerEps1 = 0.;
394  if( xy2->y == 0. ) lowerEps2 = 0.; }
395  else if( xy1->x > xy2->x ) {
396  lowerEps2 = 0.;
397  if( xy1->y == 0. ) lowerEps1 = 0.; }
398  else {
399  lowerEps1 = lowerEps2 = 0.;
400  }
401 
402  xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, n1 - 1 );
403  xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, n2 - 1 );
404  if( xy1->x < xy2->x ) {
405  upperEps2 = 0.;
406  if( xy1->y == 0. ) upperEps1 = 0.; }
407  else if( xy1->x > xy2->x ) {
408  upperEps1 = 0.;
409  if( xy2->y == 0. ) upperEps2 = 0.; }
410  else {
411  upperEps1 = upperEps2 = 0.;
412  }
413 
414  if( ( lowerEps1 != 0. ) || ( upperEps1 != 0. ) )
415  if( ( status = ptwXY_dullEdges( ptwXY1, lowerEps1, upperEps1, positiveXOnly1 ) ) != nfu_Okay ) return( status );
416  if( ( lowerEps2 != 0. ) || ( upperEps2 != 0. ) )
417  if( ( status = ptwXY_dullEdges( ptwXY2, lowerEps2, upperEps2, positiveXOnly2 ) ) != nfu_Okay ) return( status );
418 
419  return( status );
420 }
421 /*
422 ************************************************************
423 */
424 nfu_status ptwXY_copyToC_XY( ptwXYPoints *ptwXY, int64_t index1, int64_t index2, int64_t allocatedSize, int64_t *numberOfPoints, double *xys ) {
425 
426  int64_t i;
427  double *d = xys;
428  nfu_status status;
429  ptwXYPoint *pointFrom;
430 
431  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
432  if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
433 
434  if( index1 < 0 ) index1 = 0;
435  if( index2 > ptwXY->length ) index2 = ptwXY->length;
436  if( index2 < index1 ) index2 = index1;
437  *numberOfPoints = index2 - index1;
438  if( allocatedSize < ( index2 - index1 ) ) return( nfu_insufficientMemory );
439 
440  for( i = index1, pointFrom = ptwXY->points; i < index2; i++, pointFrom++ ) {
441  *(d++) = pointFrom->x;
442  *(d++) = pointFrom->y;
443  }
444 
445  return( nfu_Okay );
446 }
447 /*
448 ************************************************************
449 */
450 nfu_status ptwXY_valueTo_ptwXAndY( ptwXYPoints *ptwXY, double **xs, double **ys ) {
451 
452  int64_t i1, length = ptwXY_length( ptwXY );
453  double *xps, *yps;
454  ptwXYPoint *pointFrom;
455  nfu_status status;
456 
457  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
458  if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
459 
460  if( ( *xs = (double *) malloc( length * sizeof( double ) ) ) == NULL ) return( nfu_mallocError );
461  if( ( *ys = (double *) malloc( length * sizeof( double ) ) ) == NULL ) {
462  free( *xs );
463  *xs = NULL;
464  return( nfu_mallocError );
465  }
466 
467  for( i1 = 0, pointFrom = ptwXY->points, xps = *xs, yps = *ys; i1 < length; ++i1, ++pointFrom, ++xps, ++yps ) {
468  *xps = pointFrom->x;
469  *yps = pointFrom->y;
470  }
471 
472  return( nfu_Okay );
473 }
474 /*
475 ************************************************************
476 */
477 ptwXYPoints *ptwXY_valueTo_ptwXY( double x1, double x2, double y, nfu_status *status ) {
478 
479  ptwXYPoints *n;
480 
481  *status = nfu_XNotAscending;
482  if( x1 >= x2 ) return( NULL );
483  *status = nfu_Okay;
484  if( ( n = ptwXY_new( ptwXY_interpolationLinLin, NULL, ptwXY_maxBiSectionMax, ptwXY_minAccuracy, 2, 0, status, 0 ) ) == NULL ) return( NULL );
485  ptwXY_setValueAtX( n, x1, y );
486  ptwXY_setValueAtX( n, x2, y );
487  return( n );
488 }
489 /*
490 ************************************************************
491 */
492 ptwXYPoints *ptwXY_createGaussianCenteredSigma1( double accuracy, nfu_status *status ) {
493 
494  int64_t i, n;
495  ptwXYPoint *pm, *pp;
496  double x1, y1, x2, y2, accuracy2, yMin = 1e-10;
497  ptwXYPoints *gaussian;
498 
499  if( accuracy < 1e-5 ) accuracy = 1e-5;
500  if( accuracy > 1e-1 ) accuracy = 1e-1;
501  if( ( gaussian = ptwXY_new( ptwXY_interpolationLinLin, NULL, 1., accuracy, 200, 100, status, 0 ) ) == NULL ) return( NULL );
502  accuracy2 = accuracy = gaussian->accuracy;
503  if( accuracy2 > 5e-3 ) accuracy2 = 5e-3;
504 
505  x1 = -std::sqrt( -2. * G4Log( yMin ) );
506  y1 = yMin;
507  x2 = -5.2;
508  y2 = G4Exp( -0.5 * x2 * x2 );
509  if( ( *status = ptwXY_setValueAtX( gaussian, x1, y1 ) ) != nfu_Okay ) goto Err;
510  gaussian->accuracy = 20 * accuracy2;
511  if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
512  x1 = x2;
513  y1 = y2;
514  x2 = -4.;
515  y2 = G4Exp( -0.5 * x2 * x2 );
516  gaussian->accuracy = 5 * accuracy2;
517  if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
518  x1 = x2;
519  y1 = y2;
520  x2 = -1;
521  y2 = G4Exp( -0.5 * x2 * x2 );
522  gaussian->accuracy = accuracy;
523  if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
524  x1 = x2;
525  y1 = y2;
526  x2 = 0;
527  y2 = G4Exp( -0.5 * x2 * x2 );
528  if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
529 
530  n = gaussian->length;
531  if( ( *status = ptwXY_coalescePoints( gaussian, 2 * n + 1, NULL, 0 ) ) != nfu_Okay ) goto Err;
532  if( ( *status = ptwXY_setValueAtX( gaussian, 0., 1. ) ) != nfu_Okay ) goto Err;
533  pp = &(gaussian->points[gaussian->length]);
534  for( i = 0, pm = pp - 2; i < n; i++, pp++, pm-- ) {
535  *pp = *pm;
536  pp->x *= -1;
537  }
538  gaussian->length = 2 * n + 1;
539 
540  return( gaussian );
541 
542 Err:
543  ptwXY_free( gaussian );
544  return( NULL );
545 }
546 /*
547 ************************************************************
548 */
549 static nfu_status ptwXY_createGaussianCenteredSigma1_2( ptwXYPoints *ptwXY, double x1, double y1, double x2, double y2, int addX1Point ) {
550 
551  nfu_status status = nfu_Okay;
552  int morePoints = 0;
553  double x = 0.5 * ( x1 + x2 );
554  double y = G4Exp( -x * x / 2 ), yMin = ( y1 * ( x2 - x ) + y2 * ( x - x1 ) ) / ( x2 - x1 );
555 
556  if( std::fabs( y - yMin ) > y * ptwXY->accuracy ) morePoints = 1;
557  if( morePoints && ( status = ptwXY_createGaussianCenteredSigma1_2( ptwXY, x, y, x2, y2, 0 ) ) != nfu_Okay ) return( status );
558  if( ( status = ptwXY_setValueAtX( ptwXY, x, y ) ) != nfu_Okay ) return( status );
559  if( morePoints && ( status = ptwXY_createGaussianCenteredSigma1_2( ptwXY, x1, y1, x, y, 0 ) ) != nfu_Okay ) return( status );
560  if( addX1Point ) status = ptwXY_setValueAtX( ptwXY, x1, y1 );
561  return( status );
562 }
563 /*
564 ************************************************************
565 */
566 ptwXYPoints *ptwXY_createGaussian( double accuracy, double xCenter, double sigma, double amplitude, double xMin, double xMax,
567  double /*dullEps*/, nfu_status *status ) {
568 
569  int64_t i;
570  ptwXYPoints *gaussian, *sliced;
571  ptwXYPoint *point;
572 
573  if( ( gaussian = ptwXY_createGaussianCenteredSigma1( accuracy, status ) ) == NULL ) return( NULL );
574  for( i = 0, point = gaussian->points; i < gaussian->length; i++, point++ ) {
575  point->x = point->x * sigma + xCenter;
576  point->y *= amplitude;
577  }
578  if( ( gaussian->points[0].x < xMin ) || ( gaussian->points[gaussian->length - 1].x > xMax ) ) {
579  if( ( sliced = ptwXY_xSlice( gaussian, xMin, xMax, 10, 1, status ) ) == NULL ) goto Err;
580  ptwXY_free( gaussian );
581  gaussian = sliced;
582  }
583 
584  return( gaussian );
585 
586 Err:
587  ptwXY_free( gaussian );
588  return( NULL );
589 }
590 
591 #if defined __cplusplus
592 }
593 #endif
ptwXYPoints * ptwXY_intersectionWith_ptwX(ptwXYPoints *ptwXY, ptwXPoints *ptwX, nfu_status *status)
ptwXPoints * ptwXY_getXArray(ptwXYPoints *ptwXY, nfu_status *status)
nfu_status ptwXY_mergeClosePoints(ptwXYPoints *ptwXY, double epsilon)
int64_t ptwXY_length(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:583
nfu_status ptwXY_areDomainsMutual(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2)
ptwXYPoints * ptwXY_createGaussian(double accuracy, double xCenter, double sigma, double amplitude, double xMin, double xMax, double, nfu_status *status)
nfu_status ptwXY_valueTo_ptwXAndY(ptwXYPoints *ptwXY, double **xs, double **ys)
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
nfu_status ptwXY_mutualifyDomains(ptwXYPoints *ptwXY1, double lowerEps1, double upperEps1, int positiveXOnly1, ptwXYPoints *ptwXY2, double lowerEps2, double upperEps2, int positiveXOnly2)
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
ptwXYPoints * ptwXY_valueTo_ptwXY(double x1, double x2, double y, nfu_status *status)
ptwXPoints * ptwX_new(int64_t size, nfu_status *status)
Definition: ptwX_core.cc:24
#define minEps
ptwXYPoints * ptwXY_createGaussianCenteredSigma1(double accuracy, nfu_status *status)
#define DBL_EPSILON
Definition: templates.hh:87
const G4int n
nfu_status ptwXY_coalescePoints(ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint, int forceSmallerResize)
Definition: ptwXY_core.cc:469
ptwXYPoints * ptwXY_xSlice(ptwXYPoints *ptwXY, double xMin, double xMax, int64_t secondarySize, int fill, nfu_status *status)
Definition: ptwXY_core.cc:274
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
int64_t ptwX_length(ptwXPoints *ptwX)
Definition: ptwX_core.cc:166
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:684
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
Definition: ptwXY_core.cc:844
const G4double x[NPOINTSGL]
static nfu_status ptwXY_createGaussianCenteredSigma1_2(ptwXYPoints *ptwXY, double x1, double y1, double x2, double y2, int addX1Point)
nfu_status ptwXY_tweakDomainsToMutualify(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, int epsilonFactor, double epsilon)
nfu_status ptwXY_dullEdges(ptwXYPoints *ptwXY, double lowerEps, double upperEps, int positiveXOnly)
ptwXYPoints * ptwXY_new(ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, nfu_status *status, int userFlag)
Definition: ptwXY_core.cc:29
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
G4int sign(const T t)
A simple sign function that allows us to port fortran code to c++ more easily.
double epsilon(double density, double temperature)
nfu_status ptwXY_copyToC_XY(ptwXYPoints *ptwXY, int64_t index1, int64_t index2, int64_t allocatedSize, int64_t *numberOfPoints, double *xys)