| Geant4 Cross Reference |
1 /* 1 /*
2 # <<BEGIN-copyright>> 2 # <<BEGIN-copyright>>
3 # <<END-copyright>> 3 # <<END-copyright>>
4 */ 4 */
5 5
6 #include <cmath> 6 #include <cmath>
7 #include <float.h> 7 #include <float.h>
8 8
9 #include "ptwXY.h" 9 #include "ptwXY.h"
10 10
11 #if defined __cplusplus 11 #if defined __cplusplus
12 namespace GIDI { 12 namespace GIDI {
13 using namespace GIDI; 13 using namespace GIDI;
14 #endif 14 #endif
15 15
16 static double ptwXY_mod2( double v, double m, 16 static double ptwXY_mod2( double v, double m, int pythonMod );
17 static nfu_status ptwXY_mul2_s_ptwXY( ptwXYPoi 17 static nfu_status ptwXY_mul2_s_ptwXY( ptwXYPoints *n, ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double x1, double y1, double x2, double y2, int level );
18 static nfu_status ptwXY_div_s_ptwXY( ptwXYPoin 18 static nfu_status ptwXY_div_s_ptwXY( ptwXYPoints *n, ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double x1, double y1, double x2, double y2,
19 int level, int isNAN1, int isNAN2 ); 19 int level, int isNAN1, int isNAN2 );
20 static ptwXYPoints *ptwXY_div_ptwXY_forFlats( 20 static ptwXYPoints *ptwXY_div_ptwXY_forFlats( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int safeDivide );
21 static nfu_status ptwXY_getValueAtX_ignore_XOu 21 static nfu_status ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXYPoints *ptwXY1, double x, double *y );
22 /* 22 /*
23 ********************************************** 23 ************************************************************
24 */ 24 */
25 nfu_status ptwXY_slopeOffset( ptwXYPoints *ptw 25 nfu_status ptwXY_slopeOffset( ptwXYPoints *ptwXY, double slope, double offset ) {
26 26
27 int64_t i, nonOverflowLength = ptwXY_getNo 27 int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
28 ptwXYPoint *p; 28 ptwXYPoint *p;
29 ptwXYOverflowPoint *o, *overflowHeader = & 29 ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
30 30
31 if( ptwXY->status != nfu_Okay ) return( pt 31 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
32 32
33 for( i = 0, p = ptwXY->points; i < nonOver 33 for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = slope * p->y + offset;
34 for( o = overflowHeader->next; o != overfl 34 for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = slope * o->point.y + offset;
35 return( ptwXY->status ); 35 return( ptwXY->status );
36 } 36 }
37 /* 37 /*
38 ********************************************** 38 ************************************************************
39 */ 39 */
40 nfu_status ptwXY_add_double( ptwXYPoints *ptwX 40 nfu_status ptwXY_add_double( ptwXYPoints *ptwXY, double value ) { return( ptwXY_slopeOffset( ptwXY, 1., value ) ); }
41 nfu_status ptwXY_sub_doubleFrom( ptwXYPoints * 41 nfu_status ptwXY_sub_doubleFrom( ptwXYPoints *ptwXY, double value ) { return( ptwXY_slopeOffset( ptwXY, 1., -value ) ); }
42 nfu_status ptwXY_sub_fromDouble( ptwXYPoints * 42 nfu_status ptwXY_sub_fromDouble( ptwXYPoints *ptwXY, double value ) { return( ptwXY_slopeOffset( ptwXY, -1., value ) ); }
43 nfu_status ptwXY_mul_double( ptwXYPoints *ptwX 43 nfu_status ptwXY_mul_double( ptwXYPoints *ptwXY, double value ) { return( ptwXY_slopeOffset( ptwXY, value, 0. ) ); }
44 nfu_status ptwXY_div_doubleFrom( ptwXYPoints * 44 nfu_status ptwXY_div_doubleFrom( ptwXYPoints *ptwXY, double value ) {
45 45
46 if( value == 0. ) { 46 if( value == 0. ) {
47 ptwXY->status = nfu_divByZero; } 47 ptwXY->status = nfu_divByZero; }
48 else { 48 else {
49 ptwXY_slopeOffset( ptwXY, 1. / value, 49 ptwXY_slopeOffset( ptwXY, 1. / value, 0. );
50 } 50 }
51 return( ptwXY->status ); 51 return( ptwXY->status );
52 } 52 }
53 nfu_status ptwXY_div_fromDouble( ptwXYPoints * 53 nfu_status ptwXY_div_fromDouble( ptwXYPoints *ptwXY, double value ) {
54 /* 54 /*
55 * This does not do any infilling and it shou 55 * This does not do any infilling and it should?????????
56 */ 56 */
57 57
58 int64_t i, nonOverflowLength = ptwXY_getNo 58 int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
59 ptwXYPoint *p; 59 ptwXYPoint *p;
60 ptwXYOverflowPoint *o, *overflowHeader = & 60 ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
61 61
62 if( ptwXY->status != nfu_Okay ) return( pt 62 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
63 if( ptwXY->interpolation == ptwXY_interpol 63 if( ptwXY->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
64 64
65 for( i = 0, p = ptwXY->points; i < nonOver 65 for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) if( p->y == 0. ) ptwXY->status = nfu_divByZero;
66 for( o = overflowHeader->next; o != overfl 66 for( o = overflowHeader->next; o != overflowHeader; o = o->next ) if( o->point.y == 0. ) ptwXY->status = nfu_divByZero;
67 if( ptwXY->status != nfu_divByZero ) { 67 if( ptwXY->status != nfu_divByZero ) {
68 for( i = 0, p = ptwXY->points; i < non 68 for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = value / p->y;
69 for( o = overflowHeader->next; o != ov 69 for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = value / o->point.y;
70 } 70 }
71 return( ptwXY->status ); 71 return( ptwXY->status );
72 } 72 }
73 /* 73 /*
74 ********************************************** 74 ************************************************************
75 */ 75 */
76 nfu_status ptwXY_mod( ptwXYPoints *ptwXY, doub 76 nfu_status ptwXY_mod( ptwXYPoints *ptwXY, double m, int pythonMod ) {
77 77
78 int64_t i, nonOverflowLength = ptwXY_getNo 78 int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
79 ptwXYPoint *p; 79 ptwXYPoint *p;
80 ptwXYOverflowPoint *o, *overflowHeader = & 80 ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
81 81
82 if( ptwXY->status != nfu_Okay ) return( pt 82 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
83 if( m == 0 ) return( ptwXY->status = nfu_d 83 if( m == 0 ) return( ptwXY->status = nfu_divByZero );
84 84
85 for( i = 0, p = ptwXY->points; i < nonOver 85 for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = ptwXY_mod2( p->y, m, pythonMod );
86 for( o = overflowHeader->next; o != overfl 86 for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = ptwXY_mod2( o->point.y, m, pythonMod );
87 return( ptwXY->status ); 87 return( ptwXY->status );
88 } 88 }
89 /* 89 /*
90 ********************************************** 90 ************************************************************
91 */ 91 */
92 static double ptwXY_mod2( double v, double m, 92 static double ptwXY_mod2( double v, double m, int pythonMod ) {
93 93
94 double r = std::fmod( std::fabs( v ), std: 94 double r = std::fmod( std::fabs( v ), std::fabs( m ) );
95 95
96 if( pythonMod ) { 96 if( pythonMod ) {
97 if( ( v * m ) < 0. ) r = std::fabs( m 97 if( ( v * m ) < 0. ) r = std::fabs( m ) - std::fabs( r );
98 if( m < 0. ) r *= -1.; } 98 if( m < 0. ) r *= -1.; }
99 else { 99 else {
100 if( v < 0. ) r *= -1.; 100 if( v < 0. ) r *= -1.;
101 } 101 }
102 102
103 return( r ); 103 return( r );
104 } 104 }
105 /* 105 /*
106 ********************************************** 106 ************************************************************
107 */ 107 */
108 ptwXYPoints *ptwXY_binary_ptwXY( ptwXYPoints * 108 ptwXYPoints *ptwXY_binary_ptwXY( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double v1, double v2, double v1v2, nfu_status *status ) {
109 109
110 int64_t i; 110 int64_t i;
111 int unionOptions = ptwXY_union_fill | ptwX 111 int unionOptions = ptwXY_union_fill | ptwXY_union_mergeClosePoints;
112 double y; 112 double y;
113 ptwXYPoints *n; 113 ptwXYPoints *n;
114 ptwXYPoint *p; 114 ptwXYPoint *p;
115 115
116 *status = nfu_otherInterpolation; 116 *status = nfu_otherInterpolation;
117 if( ptwXY1->interpolation == ptwXY_interpo 117 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
118 if( ptwXY2->interpolation == ptwXY_interpo 118 if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
119 if( ( *status = ptwXY_areDomainsMutual( pt 119 if( ( *status = ptwXY_areDomainsMutual( ptwXY1, ptwXY2 ) ) != nfu_Okay ) return( NULL );
120 if( ( ptwXY1->interpolation == ptwXY_inter 120 if( ( ptwXY1->interpolation == ptwXY_interpolationFlat ) || ( ptwXY2->interpolation == ptwXY_interpolationFlat ) ) {
121 *status = nfu_invalidInterpolation; 121 *status = nfu_invalidInterpolation;
122 if( ( ptwXY1->interpolation != ptwXY2- 122 if( ( ptwXY1->interpolation != ptwXY2->interpolation ) ) return( NULL );
123 } 123 }
124 if( ( n = ptwXY_union( ptwXY1, ptwXY2, sta 124 if( ( n = ptwXY_union( ptwXY1, ptwXY2, status, unionOptions ) ) != NULL ) {
125 for( i = 0, p = n->points; i < n->leng 125 for( i = 0, p = n->points; i < n->length; i++, p++ ) {
126 if( ( *status = ptwXY_getValueAtX_ 126 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, p->x, &y ) ) != nfu_Okay ) goto Err;
127 p->y = v1 * p->y + v2 * y + v1v2 * 127 p->y = v1 * p->y + v2 * y + v1v2 * y * p->y;
128 } 128 }
129 } 129 }
130 return( n ); 130 return( n );
131 Err: 131 Err:
132 if( n ) ptwXY_free( n ); 132 if( n ) ptwXY_free( n );
133 return( NULL ); 133 return( NULL );
134 } 134 }
135 /* 135 /*
136 ********************************************** 136 ************************************************************
137 */ 137 */
138 ptwXYPoints *ptwXY_add_ptwXY( ptwXYPoints *ptw 138 ptwXYPoints *ptwXY_add_ptwXY( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status ) {
139 139
140 ptwXYPoints *sum; 140 ptwXYPoints *sum;
141 141
142 if( ptwXY1->length == 0 ) { 142 if( ptwXY1->length == 0 ) {
143 sum = ptwXY_clone( ptwXY2, status ); } 143 sum = ptwXY_clone( ptwXY2, status ); }
144 else if( ptwXY2->length == 0 ) { 144 else if( ptwXY2->length == 0 ) {
145 sum = ptwXY_clone( ptwXY1, status ); } 145 sum = ptwXY_clone( ptwXY1, status ); }
146 else { 146 else {
147 sum = ptwXY_binary_ptwXY( ptwXY1, ptwX 147 sum = ptwXY_binary_ptwXY( ptwXY1, ptwXY2, 1., 1., 0., status );
148 } 148 }
149 return( sum ); 149 return( sum );
150 } 150 }
151 /* 151 /*
152 ********************************************** 152 ************************************************************
153 */ 153 */
154 ptwXYPoints *ptwXY_sub_ptwXY( ptwXYPoints *ptw 154 ptwXYPoints *ptwXY_sub_ptwXY( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status ) {
155 155
156 ptwXYPoints *diff; 156 ptwXYPoints *diff;
157 157
158 if( ptwXY1->length == 0 ) { 158 if( ptwXY1->length == 0 ) {
159 diff = ptwXY_clone( ptwXY2, status ); 159 diff = ptwXY_clone( ptwXY2, status );
160 if( ( *status = ptwXY_neg( diff ) ) != 160 if( ( *status = ptwXY_neg( diff ) ) != nfu_Okay ) diff = ptwXY_free( diff ); }
161 else if( ptwXY2->length == 0 ) { 161 else if( ptwXY2->length == 0 ) {
162 diff = ptwXY_clone( ptwXY1, status ); 162 diff = ptwXY_clone( ptwXY1, status ); }
163 else { 163 else {
164 diff = ptwXY_binary_ptwXY( ptwXY1, ptw 164 diff = ptwXY_binary_ptwXY( ptwXY1, ptwXY2, 1., -1., 0., status );
165 } 165 }
166 return( diff ); 166 return( diff );
167 } 167 }
168 /* 168 /*
169 ********************************************** 169 ************************************************************
170 */ 170 */
171 ptwXYPoints *ptwXY_mul_ptwXY( ptwXYPoints *ptw 171 ptwXYPoints *ptwXY_mul_ptwXY( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status ) {
172 172
173 ptwXYPoints *mul; 173 ptwXYPoints *mul;
174 174
175 if( ptwXY1->length == 0 ) { 175 if( ptwXY1->length == 0 ) {
176 mul = ptwXY_clone( ptwXY1, status ); } 176 mul = ptwXY_clone( ptwXY1, status ); }
177 else if( ptwXY2->length == 0 ) { 177 else if( ptwXY2->length == 0 ) {
178 mul = ptwXY_clone( ptwXY2, status ); } 178 mul = ptwXY_clone( ptwXY2, status ); }
179 else { 179 else {
180 mul = ptwXY_binary_ptwXY( ptwXY1, ptwX 180 mul = ptwXY_binary_ptwXY( ptwXY1, ptwXY2, 0., 0., 1., status );
181 } 181 }
182 return( mul ); 182 return( mul );
183 } 183 }
184 /* 184 /*
185 ********************************************** 185 ************************************************************
186 */ 186 */
187 ptwXYPoints *ptwXY_mul2_ptwXY( ptwXYPoints *pt 187 ptwXYPoints *ptwXY_mul2_ptwXY( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status ) {
188 188
189 int64_t i, length; 189 int64_t i, length;
190 ptwXYPoints *n = NULL; 190 ptwXYPoints *n = NULL;
191 int found; 191 int found;
192 double x1, y1, x2, y2, u1, u2, v1, v2, xz1 192 double x1, y1, x2, y2, u1, u2, v1, v2, xz1 = 0, xz2 = 0, x;
193 193
194 *status = nfu_otherInterpolation; 194 *status = nfu_otherInterpolation;
195 if( ptwXY1->interpolation == ptwXY_interpo 195 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
196 if( ptwXY2->interpolation == ptwXY_interpo 196 if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
197 if( ( n = ptwXY_mul_ptwXY( ptwXY1, ptwXY2, 197 if( ( n = ptwXY_mul_ptwXY( ptwXY1, ptwXY2, status ) ) == NULL ) return( n );
198 if( ptwXY1->interpolation == ptwXY_interpo 198 if( ptwXY1->interpolation == ptwXY_interpolationFlat ) return( n );
199 if( ptwXY2->interpolation == ptwXY_interpo 199 if( ptwXY2->interpolation == ptwXY_interpolationFlat ) return( n );
200 length = n->length - 1; 200 length = n->length - 1;
201 if( length > 0 ) { 201 if( length > 0 ) {
202 x2 = n->points[length].x; 202 x2 = n->points[length].x;
203 for( i = length - 1; i >= 0; i-- ) { 203 for( i = length - 1; i >= 0; i-- ) { /* Find and add y zeros not currently in n's. */
204 x1 = n->points[i].x; 204 x1 = n->points[i].x;
205 if( ( *status = ptwXY_getValueAtX_ 205 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) goto Err;
206 if( ( *status = ptwXY_getValueAtX_ 206 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) goto Err;
207 if( ( *status = ptwXY_getValueAtX_ 207 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x1, &v1 ) ) != nfu_Okay ) goto Err;
208 if( ( *status = ptwXY_getValueAtX_ 208 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x2, &v2 ) ) != nfu_Okay ) goto Err;
209 found = 0; 209 found = 0;
210 if( u1 * u2 < 0 ) { 210 if( u1 * u2 < 0 ) {
211 xz1 = ( u1 * x2 - u2 * x1 ) / 211 xz1 = ( u1 * x2 - u2 * x1 ) / ( u1 - u2 );
212 if( ( *status = ptwXY_setValue 212 if( ( *status = ptwXY_setValueAtX( n, xz1, 0. ) ) != nfu_Okay ) goto Err;
213 found = 1; 213 found = 1;
214 } 214 }
215 if( v1 * v2 < 0 ) { 215 if( v1 * v2 < 0 ) {
216 xz2 = ( v1 * x2 - v2 * x1 ) / 216 xz2 = ( v1 * x2 - v2 * x1 ) / ( v1 - v2 );
217 if( ( *status = ptwXY_setValue 217 if( ( *status = ptwXY_setValueAtX( n, xz2, 0. ) ) != nfu_Okay ) goto Err;
218 found += 1; 218 found += 1;
219 } 219 }
220 if( found > 1 ) { 220 if( found > 1 ) {
221 x = 0.5 * ( xz1 + xz2 ); 221 x = 0.5 * ( xz1 + xz2 );
222 if( ( *status = ptwXY_getValue 222 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x, &u1 ) ) != nfu_Okay ) goto Err;
223 if( ( *status = ptwXY_getValue 223 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x, &v1 ) ) != nfu_Okay ) goto Err;
224 if( ( *status = ptwXY_setValue 224 if( ( *status = ptwXY_setValueAtX( n, x, u1 * v1 ) ) != nfu_Okay ) goto Err;
225 } 225 }
226 x2 = x1; 226 x2 = x1;
227 } 227 }
228 228
229 if( ( *status = ptwXY_simpleCoalescePo 229 if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
230 length = n->length; 230 length = n->length;
231 x2 = n->points[n->length-1].x; 231 x2 = n->points[n->length-1].x;
232 y2 = n->points[n->length-1].y; 232 y2 = n->points[n->length-1].y;
233 for( i = n->length - 2; i >= 0; i-- ) 233 for( i = n->length - 2; i >= 0; i-- ) { /* Make interpolation fit accuracy. Work backwards so new points will not mess up loop. */
234 x1 = n->points[i].x; 234 x1 = n->points[i].x;
235 y1 = n->points[i].y; 235 y1 = n->points[i].y;
236 if( ( *status = ptwXY_mul2_s_ptwXY 236 if( ( *status = ptwXY_mul2_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x2, y2, 0 ) ) != nfu_Okay ) goto Err;
237 x2 = x1; 237 x2 = x1;
238 y2 = y1; 238 y2 = y1;
239 } 239 }
240 ptwXY_update_biSectionMax( n, (double) 240 ptwXY_update_biSectionMax( n, (double) length );
241 } 241 }
242 return( n ); 242 return( n );
243 243
244 Err: 244 Err:
245 if( n ) ptwXY_free( n ); 245 if( n ) ptwXY_free( n );
246 return( NULL ); 246 return( NULL );
247 } 247 }
248 /* 248 /*
249 ********************************************** 249 ************************************************************
250 */ 250 */
251 static nfu_status ptwXY_mul2_s_ptwXY( ptwXYPoi 251 static nfu_status ptwXY_mul2_s_ptwXY( ptwXYPoints *n, ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double x1, double y1, double x2, double y2, int level ) {
252 252
253 nfu_status status; 253 nfu_status status;
254 double u1, u2, v1, v2, x, y, yp, dx, a1, a 254 double u1, u2, v1, v2, x, y, yp, dx, a1, a2;
255 255
256 if( ( x2 - x1 ) < ClosestAllowXFactor * DB 256 if( ( x2 - x1 ) < ClosestAllowXFactor * DBL_EPSILON * ( std::fabs( x1 ) + std::fabs( x2 ) ) ) return( nfu_Okay );
257 if( level >= n->biSectionMax ) return( nfu 257 if( level >= n->biSectionMax ) return( nfu_Okay );
258 level++; 258 level++;
259 if( ( status = ptwXY_getValueAtX_ignore_XO 259 if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) return( status );
260 if( ( status = ptwXY_getValueAtX_ignore_XO 260 if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) return( status );
261 if( ( status = ptwXY_getValueAtX_ignore_XO 261 if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x1, &v1 ) ) != nfu_Okay ) return( status );
262 if( ( status = ptwXY_getValueAtX_ignore_XO 262 if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x2, &v2 ) ) != nfu_Okay ) return( status );
263 if( ( u1 == u2 ) || ( v1 == v2 ) ) return( 263 if( ( u1 == u2 ) || ( v1 == v2 ) ) return( nfu_Okay );
264 a1 = u1 * v1; 264 a1 = u1 * v1;
265 if( y1 == 0 ) a1 = 0.; 265 if( y1 == 0 ) a1 = 0.; /* Fix rounding problem. */
266 a2 = u2 * v2; 266 a2 = u2 * v2;
267 if( y2 == 0 ) a2 = 0.; 267 if( y2 == 0 ) a2 = 0.; /* Fix rounding problem. */
268 if( ( a1 == 0. ) || ( a2 == 0. ) ) { 268 if( ( a1 == 0. ) || ( a2 == 0. ) ) { /* Handle special case of 0 where accuracy can never be met. */
269 x = 0.5 * ( x1 + x2 ); } 269 x = 0.5 * ( x1 + x2 ); }
270 else { 270 else {
271 if( ( a1 * a2 < 0. ) ) return( nfu_Oka 271 if( ( a1 * a2 < 0. ) ) return( nfu_Okay ); /* Assume rounding error and no point needed as zero crossings are set in ptwXY_mul2_ptwXY. */
272 a1 = std::sqrt( std::fabs( a1 ) ); 272 a1 = std::sqrt( std::fabs( a1 ) );
273 a2 = std::sqrt( std::fabs( a2 ) ); 273 a2 = std::sqrt( std::fabs( a2 ) );
274 x = ( a2 * x1 + a1 * x2 ) / ( a2 + a1 274 x = ( a2 * x1 + a1 * x2 ) / ( a2 + a1 );
275 } 275 }
276 dx = x2 - x1; 276 dx = x2 - x1;
277 yp = ( u1 * v1 * ( x2 - x ) + u2 * v2 * ( 277 yp = ( u1 * v1 * ( x2 - x ) + u2 * v2 * ( x - x1 ) ) / dx;
278 y = ( u1 * ( x2 - x ) + u2 * ( x - x1 ) ) 278 y = ( u1 * ( x2 - x ) + u2 * ( x - x1 ) ) * ( v1 * ( x2 - x ) + v2 * ( x - x1 ) ) / ( dx * dx );
279 if( std::fabs( y - yp ) < std::fabs( y * n 279 if( std::fabs( y - yp ) < std::fabs( y * n->accuracy ) ) return( nfu_Okay );
280 if( ( status = ptwXY_setValueAtX( n, x, y 280 if( ( status = ptwXY_setValueAtX( n, x, y ) ) != nfu_Okay ) return( status );
281 if( ( status = ptwXY_mul2_s_ptwXY( n, ptwX 281 if( ( status = ptwXY_mul2_s_ptwXY( n, ptwXY1, ptwXY2, x, y, x2, y2, level ) ) != nfu_Okay ) return( status );
282 status = ptwXY_mul2_s_ptwXY( n, ptwXY1, pt 282 status = ptwXY_mul2_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x, y, level );
283 return( status ); 283 return( status );
284 } 284 }
285 /* 285 /*
286 ********************************************** 286 ************************************************************
287 */ 287 */
288 ptwXYPoints *ptwXY_div_ptwXY( ptwXYPoints *ptw 288 ptwXYPoints *ptwXY_div_ptwXY( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int safeDivide ) {
289 289
290 int isNAN1, isNAN2; 290 int isNAN1, isNAN2;
291 int64_t i, j, k, zeros = 0, length, iYs; 291 int64_t i, j, k, zeros = 0, length, iYs;
292 double x1, x2, y1, y2, u1, u2, v1, v2, y, 292 double x1, x2, y1, y2, u1, u2, v1, v2, y, xz, nan = nfu_getNAN( ), s1, s2;
293 ptwXYPoints *n = NULL; 293 ptwXYPoints *n = NULL;
294 ptwXYPoint *p; 294 ptwXYPoint *p;
295 295
296 if( ( *status = ptwXY_simpleCoalescePoints 296 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
297 if( ( *status = ptwXY_simpleCoalescePoints 297 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
298 *status = nfu_otherInterpolation; 298 *status = nfu_otherInterpolation;
299 if( ptwXY1->interpolation == ptwXY_interpo 299 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
300 if( ptwXY2->interpolation == ptwXY_interpo 300 if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
301 if( ( ptwXY1->interpolation == ptwXY_inter << 301 if( ( ptwXY1->interpolation == ptwXY_interpolationFlat ) || ( ptwXY1->interpolation == ptwXY_interpolationFlat ) )
302 return( ptwXY_div_ptwXY_forFlats( ptwX 302 return( ptwXY_div_ptwXY_forFlats( ptwXY1, ptwXY2, status, safeDivide ) );
303 303
304 if( ( *status = ptwXY_areDomainsMutual( pt 304 if( ( *status = ptwXY_areDomainsMutual( ptwXY1, ptwXY2 ) ) != nfu_Okay ) return( NULL );
305 if( ( n = ptwXY_union( ptwXY1, ptwXY2, sta 305 if( ( n = ptwXY_union( ptwXY1, ptwXY2, status, ptwXY_union_fill | ptwXY_union_mergeClosePoints ) ) != NULL ) {
306 for( i = 0, p = n->points; i < n->leng 306 for( i = 0, p = n->points; i < n->length; i++, p++ ) {
307 if( ( *status = ptwXY_getValueAtX_ 307 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, p->x, &y ) ) != nfu_Okay ) goto Err;
308 if( y == 0. ) { 308 if( y == 0. ) {
309 if( p->y == 0. ) { 309 if( p->y == 0. ) {
310 iYs = 0; 310 iYs = 0;
311 y1 = 0.; 311 y1 = 0.;
312 y2 = 0.; 312 y2 = 0.;
313 if( i > 0 ) { 313 if( i > 0 ) {
314 if( ( *status = ptwXY_ 314 if( ( *status = ptwXY_getSlopeAtX( ptwXY1, p->x, '-', &s1 ) ) != nfu_Okay ) {
315 if( *status != nfu 315 if( *status != nfu_XOutsideDomain ) goto Err;
316 s1 = 0.; 316 s1 = 0.;
317 } 317 }
318 if( ( *status = ptwXY_ 318 if( ( *status = ptwXY_getSlopeAtX( ptwXY2, p->x, '-', &s2 ) ) != nfu_Okay ) goto Err;
319 if( s2 == 0. ) { 319 if( s2 == 0. ) {
320 y1 = nan; } 320 y1 = nan; }
321 else { 321 else {
322 y1 = s1 / s2; 322 y1 = s1 / s2;
323 } 323 }
324 iYs++; 324 iYs++;
325 } 325 }
326 if( i < ( n->length - 1 ) 326 if( i < ( n->length - 1 ) ) {
327 if( ( *status = ptwXY_ 327 if( ( *status = ptwXY_getSlopeAtX( ptwXY1, p->x, '+', &s1 ) ) != nfu_Okay ) {
328 if( *status != nfu 328 if( *status != nfu_XOutsideDomain ) goto Err;
329 s1 = 0.; 329 s1 = 0.;
330 } 330 }
331 if( ( *status = ptwXY_ 331 if( ( *status = ptwXY_getSlopeAtX( ptwXY2, p->x, '+', &s2 ) ) != nfu_Okay ) goto Err;
332 if( s2 == 0. ) { 332 if( s2 == 0. ) {
333 y2 = nan; } 333 y2 = nan; }
334 else { 334 else {
335 y2 = s1 / s2; 335 y2 = s1 / s2;
336 } 336 }
337 iYs++; 337 iYs++;
338 } 338 }
339 p->y = ( y1 + y2 ) / iYs; 339 p->y = ( y1 + y2 ) / iYs;
340 if( nfu_isNAN( p->y ) ) ze 340 if( nfu_isNAN( p->y ) ) zeros++; }
341 else { 341 else {
342 if( !safeDivide ) { 342 if( !safeDivide ) {
343 *status = nfu_divByZer 343 *status = nfu_divByZero;
344 goto Err; 344 goto Err;
345 } 345 }
346 zeros++; 346 zeros++;
347 p->y = nan; 347 p->y = nan;
348 } } 348 } }
349 else { 349 else {
350 p->y /= y; 350 p->y /= y;
351 } 351 }
352 } 352 }
353 length = n->length - 1; 353 length = n->length - 1;
354 if( length > 0 ) { 354 if( length > 0 ) {
355 x2 = n->points[length].x; 355 x2 = n->points[length].x;
356 for( i = length - 1; i >= 0; i-- ) 356 for( i = length - 1; i >= 0; i-- ) { /* Find and add y zeros and NAN not currently in n's. */
357 x1 = n->points[i].x; 357 x1 = n->points[i].x;
358 if( ( *status = ptwXY_getValue 358 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) goto Err;
359 if( ( *status = ptwXY_getValue 359 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) goto Err;
360 if( ( *status = ptwXY_getValue 360 if( ( *status = ptwXY_getValueAtX( ptwXY2, x1, &v1 ) ) != nfu_Okay ) goto Err;
361 if( ( *status = ptwXY_getValue 361 if( ( *status = ptwXY_getValueAtX( ptwXY2, x2, &v2 ) ) != nfu_Okay ) goto Err;
362 if( u1 * u2 < 0 ) { 362 if( u1 * u2 < 0 ) {
363 xz = ( u1 * x2 - u2 * x1 ) 363 xz = ( u1 * x2 - u2 * x1 ) / ( u1 - u2 );
364 if( ( *status = ptwXY_setV 364 if( ( *status = ptwXY_setValueAtX( n, xz, 0. ) ) != nfu_Okay ) goto Err;
365 } 365 }
366 if( v1 * v2 < 0 ) { 366 if( v1 * v2 < 0 ) {
367 if( !safeDivide ) { 367 if( !safeDivide ) {
368 *status = nfu_divByZer 368 *status = nfu_divByZero;
369 goto Err; 369 goto Err;
370 } 370 }
371 zeros++; 371 zeros++;
372 xz = ( v1 * x2 - v2 * x1 ) 372 xz = ( v1 * x2 - v2 * x1 ) / ( v1 - v2 );
373 if( ( *status = ptwXY_setV 373 if( ( *status = ptwXY_setValueAtX( n, xz, nan ) ) != nfu_Okay ) goto Err;
374 } 374 }
375 x2 = x1; 375 x2 = x1;
376 } 376 }
377 if( ( *status = ptwXY_simpleCoales 377 if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
378 length = n->length; 378 length = n->length;
379 x2 = n->points[n->length-1].x; 379 x2 = n->points[n->length-1].x;
380 y2 = n->points[n->length-1].y; 380 y2 = n->points[n->length-1].y;
381 isNAN2 = nfu_isNAN( y2 ); 381 isNAN2 = nfu_isNAN( y2 );
382 for( i = n->length - 2; i >= 0; i- 382 for( i = n->length - 2; i >= 0; i-- ) { /* Make interpolation fit accuracy. Work backwards so new points will not mess up loop. */
383 x1 = n->points[i].x; 383 x1 = n->points[i].x;
384 y1 = n->points[i].y; 384 y1 = n->points[i].y;
385 isNAN1 = nfu_isNAN( y1 ); 385 isNAN1 = nfu_isNAN( y1 );
386 if( !isNAN1 || !isNAN2 ) { 386 if( !isNAN1 || !isNAN2 ) {
387 if( ( *status = ptwXY_div_ 387 if( ( *status = ptwXY_div_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x2, y2, 0, isNAN1, isNAN2 ) ) != nfu_Okay ) goto Err;
388 } 388 }
389 x2 = x1; 389 x2 = x1;
390 y2 = y1; 390 y2 = y1;
391 isNAN2 = isNAN1; 391 isNAN2 = isNAN1;
392 } 392 }
393 ptwXY_update_biSectionMax( n, (dou 393 ptwXY_update_biSectionMax( n, (double) length );
394 if( zeros ) { 394 if( zeros ) {
395 if( ( *status = ptwXY_simpleCo 395 if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
396 for( i = 0; i < n->length; i++ 396 for( i = 0; i < n->length; i++ ) if( !nfu_isNAN( n->points[i].y ) ) break;
397 if( nfu_isNAN( n->points[0].y 397 if( nfu_isNAN( n->points[0].y ) ) { /* Special case for first point. */
398 if( i == n->length ) { 398 if( i == n->length ) { /* They are all nan's, what now? */
399 zeros = 0; 399 zeros = 0;
400 for( i = 0; i < n->len 400 for( i = 0; i < n->length; i++ ) n->points[i].y = 0.; }
401 else { 401 else {
402 n->points[0].y = 2. * 402 n->points[0].y = 2. * n->points[i].y;
403 zeros--; 403 zeros--;
404 } 404 }
405 } 405 }
406 for( i = n->length - 1; i > 0; 406 for( i = n->length - 1; i > 0; i-- ) if( !nfu_isNAN( n->points[i].y ) ) break;
407 if( nfu_isNAN( n->points[n->le 407 if( nfu_isNAN( n->points[n->length - 1].y ) ) { /* Special case for last point. */
408 n->points[n->length - 1].y 408 n->points[n->length - 1].y = 2. * n->points[i].y;
409 zeros--; 409 zeros--;
410 } 410 }
411 if( zeros ) { 411 if( zeros ) {
412 for( i = 0; i < n->length; 412 for( i = 0; i < n->length; i++ ) if( nfu_isNAN( n->points[i].y ) ) break;
413 for( k = i + 1, j = i; k < 413 for( k = i + 1, j = i; k < n->length; k++ ) {
414 if( nfu_isNAN( n->points[k 414 if( nfu_isNAN( n->points[k].y ) ) continue;
415 n->points[j] = n->poin 415 n->points[j] = n->points[k];
416 j++; 416 j++;
417 } 417 }
418 n->length = j; 418 n->length = j;
419 } 419 }
420 } 420 }
421 } 421 }
422 } 422 }
423 return( n ); 423 return( n );
424 424
425 Err: 425 Err:
426 if( n ) ptwXY_free( n ); 426 if( n ) ptwXY_free( n );
427 return( NULL ); 427 return( NULL );
428 } 428 }
429 /* 429 /*
430 ********************************************** 430 ************************************************************
431 */ 431 */
432 static nfu_status ptwXY_div_s_ptwXY( ptwXYPoin 432 static nfu_status ptwXY_div_s_ptwXY( ptwXYPoints *n, ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double x1, double y1, double x2, double y2,
433 int level, int isNAN1, int isNAN2 ) { 433 int level, int isNAN1, int isNAN2 ) {
434 434
435 nfu_status status; 435 nfu_status status;
436 double u1, u2, v1, v2, v, x, y, yp, dx, a1 436 double u1, u2, v1, v2, v, x, y, yp, dx, a1, a2;
437 437
438 if( ( x2 - x1 ) < ClosestAllowXFactor * DB 438 if( ( x2 - x1 ) < ClosestAllowXFactor * DBL_EPSILON * ( std::fabs( x1 ) + std::fabs( x2 ) ) ) return( nfu_Okay );
439 if( level >= n->biSectionMax ) return( nfu 439 if( level >= n->biSectionMax ) return( nfu_Okay );
440 level++; 440 level++;
441 if( ( status = ptwXY_getValueAtX_ignore_XO 441 if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) return( status );
442 if( ( status = ptwXY_getValueAtX_ignore_XO 442 if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) return( status );
443 if( ( status = ptwXY_getValueAtX( ptwXY2, 443 if( ( status = ptwXY_getValueAtX( ptwXY2, x1, &v1 ) ) != nfu_Okay ) return( status );
444 if( ( status = ptwXY_getValueAtX( ptwXY2, 444 if( ( status = ptwXY_getValueAtX( ptwXY2, x2, &v2 ) ) != nfu_Okay ) return( status );
445 if( isNAN1 ) { 445 if( isNAN1 ) {
446 x = 0.5 * ( x1 + x2 ); 446 x = 0.5 * ( x1 + x2 );
447 if( ( status = ptwXY_getValueAtX_ignor 447 if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x, &u1 ) ) != nfu_Okay ) return( status );
448 if( ( status = ptwXY_getValueAtX( ptwX 448 if( ( status = ptwXY_getValueAtX( ptwXY2, x, &v1 ) ) != nfu_Okay ) return( status );
449 y = u1 / v1; } 449 y = u1 / v1; }
450 else if( isNAN2 ) { 450 else if( isNAN2 ) {
451 x = 0.5 * ( x1 + x2 ); 451 x = 0.5 * ( x1 + x2 );
452 if( ( status = ptwXY_getValueAtX_ignor 452 if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x, &u2 ) ) != nfu_Okay ) return( status );
453 if( ( status = ptwXY_getValueAtX( ptwX 453 if( ( status = ptwXY_getValueAtX( ptwXY2, x, &v2 ) ) != nfu_Okay ) return( status );
454 y = u2 / v2; } 454 y = u2 / v2; }
455 else { 455 else {
456 if( ( u1 == u2 ) || ( v1 == v2 ) ) ret 456 if( ( u1 == u2 ) || ( v1 == v2 ) ) return( nfu_Okay );
457 if( ( y1 == 0. ) || ( y2 == 0. ) ) { 457 if( ( y1 == 0. ) || ( y2 == 0. ) ) { /* Handle special case of 0 where accuracy can never be met. */
458 x = 0.5 * ( x1 + x2 ); } 458 x = 0.5 * ( x1 + x2 ); }
459 else { 459 else {
460 if( ( u1 * u2 < 0. ) ) return( nfu 460 if( ( u1 * u2 < 0. ) ) return( nfu_Okay ); /* Assume rounding error and no point needed. */
461 a1 = std::sqrt( std::fabs( u1 ) ); 461 a1 = std::sqrt( std::fabs( u1 ) );
462 a2 = std::sqrt( std::fabs( u2 ) ); 462 a2 = std::sqrt( std::fabs( u2 ) );
463 x = ( a2 * x1 + a1 * x2 ) / ( a2 + 463 x = ( a2 * x1 + a1 * x2 ) / ( a2 + a1 );
464 } 464 }
465 dx = x2 - x1; 465 dx = x2 - x1;
466 v = v1 * ( x2 - x ) + v2 * ( x - x1 ); 466 v = v1 * ( x2 - x ) + v2 * ( x - x1 );
467 if( ( v1 == 0. ) || ( v2 == 0. ) || ( 467 if( ( v1 == 0. ) || ( v2 == 0. ) || ( v == 0. ) ) return( nfu_Okay ); /* Probably not correct, but I had to do something. */
468 yp = ( u1 / v1 * ( x2 - x ) + u2 / v2 468 yp = ( u1 / v1 * ( x2 - x ) + u2 / v2 * ( x - x1 ) ) / dx;
469 y = ( u1 * ( x2 - x ) + u2 * ( x - x1 469 y = ( u1 * ( x2 - x ) + u2 * ( x - x1 ) ) / v;
470 if( std::fabs( y - yp ) < std::fabs( y 470 if( std::fabs( y - yp ) < std::fabs( y * n->accuracy ) ) return( nfu_Okay );
471 } 471 }
472 if( ( status = ptwXY_setValueAtX( n, x, y 472 if( ( status = ptwXY_setValueAtX( n, x, y ) ) != nfu_Okay ) return( status );
473 if( ( status = ptwXY_div_s_ptwXY( n, ptwXY 473 if( ( status = ptwXY_div_s_ptwXY( n, ptwXY1, ptwXY2, x, y, x2, y2, level, 0, isNAN2 ) ) != nfu_Okay ) return( status );
474 status = ptwXY_div_s_ptwXY( n, ptwXY1, ptw 474 status = ptwXY_div_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x, y, level, isNAN1, 0 );
475 return( status ); 475 return( status );
476 } 476 }
477 /* 477 /*
478 ********************************************** 478 ************************************************************
479 */ 479 */
480 static ptwXYPoints *ptwXY_div_ptwXY_forFlats( 480 static ptwXYPoints *ptwXY_div_ptwXY_forFlats( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int safeDivide ) {
481 481
482 int64_t i; 482 int64_t i;
483 ptwXYPoints *n = NULL; 483 ptwXYPoints *n = NULL;
484 ptwXYPoint *p; 484 ptwXYPoint *p;
485 double y; 485 double y;
486 486
487 *status = nfu_invalidInterpolation; 487 *status = nfu_invalidInterpolation;
488 if( ptwXY1->interpolation != ptwXY_interpo 488 if( ptwXY1->interpolation != ptwXY_interpolationFlat ) return( NULL );
489 if( ptwXY2->interpolation != ptwXY_interpo 489 if( ptwXY2->interpolation != ptwXY_interpolationFlat ) return( NULL );
490 if( ( n = ptwXY_union( ptwXY1, ptwXY2, sta 490 if( ( n = ptwXY_union( ptwXY1, ptwXY2, status, ptwXY_union_fill | ptwXY_union_mergeClosePoints ) ) != NULL ) {
491 for( i = 0, p = n->points; i < n->leng 491 for( i = 0, p = n->points; i < n->length; i++, p++ ) {
492 if( ( *status = ptwXY_getValueAtX_ 492 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, p->x, &y ) ) != nfu_Okay ) goto Err;
493 if( y == 0. ) { 493 if( y == 0. ) {
494 if( ( safeDivide ) && ( p->y = 494 if( ( safeDivide ) && ( p->y == 0 ) ) {
495 *status = nfu_divByZero; 495 *status = nfu_divByZero;
496 goto Err; 496 goto Err;
497 } } 497 } }
498 else { 498 else {
499 p->y /= y; 499 p->y /= y;
500 } 500 }
501 } 501 }
502 } 502 }
503 return( n ); 503 return( n );
504 504
505 Err: 505 Err:
506 if( n ) ptwXY_free( n ); 506 if( n ) ptwXY_free( n );
507 return( NULL ); 507 return( NULL );
508 } 508 }
509 /* 509 /*
510 ********************************************** 510 ************************************************************
511 */ 511 */
512 static nfu_status ptwXY_getValueAtX_ignore_XOu 512 static nfu_status ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXYPoints *ptwXY1, double x, double *y ) {
513 513
514 nfu_status status = ptwXY_getValueAtX( ptw 514 nfu_status status = ptwXY_getValueAtX( ptwXY1, x, y );
515 515
516 if( status == nfu_XOutsideDomain ) status 516 if( status == nfu_XOutsideDomain ) status = nfu_Okay;
517 return( status ); 517 return( status );
518 } 518 }
519 519
520 #if defined __cplusplus 520 #if defined __cplusplus
521 } 521 }
522 #endif 522 #endif
523 523