Enable monotonic negative growth of abscissa values.

This commit is contained in:
Holger Vogt 2024-06-09 17:29:57 +02:00
parent b3fb5c00b3
commit 8d854fcbfc
2 changed files with 45 additions and 14 deletions

View File

@ -1078,13 +1078,31 @@ static INPparseNode *prepare_PTF_PWL(INPparseNode *p)
fprintf(stderr, " (%lf %lf)\n", data->vals[i], data->vals[i+1]);
#endif
for (i = 2 ; i < data->n ; i += 2)
if(data->vals[i-2] >= data->vals[i]) {
fprintf(stderr,
"Error: PWL(expr, points...) the abscissa of points "
"must be ascending at line %d\nfrom file\n %s\n", Current_parse_line, Sourcefile);
controlled_exit(EXIT_BAD);
}
/* check for monotonic abscissa */
if (data->vals[0] > data->vals[2]) {
for (i = 2; i < data->n; i += 2)
if (data->vals[i - 2] < data->vals[i]) {
fprintf(stderr,
"Error: PWL(expr, points...) the abscissa of points "
"must be descending at line %d\nfrom file\n %s\n", Current_parse_line, Sourcefile);
controlled_exit(EXIT_BAD);
}
}
else if (data->vals[0] < data->vals[2]) {
for (i = 2; i < data->n; i += 2)
if (data->vals[i - 2] > data->vals[i]) {
fprintf(stderr,
"Error: PWL(expr, points...) the abscissa of points "
"must be ascending at line %d\nfrom file\n %s\n", Current_parse_line, Sourcefile);
controlled_exit(EXIT_BAD);
}
}
else {
fprintf(stderr,
"Error: PWL(expr, points...) the abscissa of points "
"must be monotonic at line %d\nfrom file\n %s\n", Current_parse_line, Sourcefile);
controlled_exit(EXIT_BAD);
}
/* strip all but the first arg,
* and attach the rest as opaque data to the INPparseNode

View File

@ -351,14 +351,27 @@ PTpwl(double arg, void *data)
int k0 = 0;
int k1 = thing->n/2 - 1;
while(k1-k0 > 1) {
int k = (k0+k1)/2;
if(thing->vals[2*k] > arg)
k1 = k;
else
k0 = k;
/* monotonically increasing abscissa */
if (thing->vals[0] < thing->vals[2]) {
while (k1 - k0 > 1) {
int k = (k0 + k1) / 2;
if (thing->vals[2 * k] > arg)
k1 = k;
else
k0 = k;
}
}
/* monotonically decreasing abscissa */
else {
while (k1 - k0 > 1) {
int k = (k0 + k1) / 2;
if (thing->vals[2 * k] < arg)
k1 = k;
else
k0 = k;
}
}
/* interpolate the ordinate */
y = thing->vals[2*k0+1] +
(thing->vals[2*k1+1] - thing->vals[2*k0+1]) *
(arg - thing->vals[2*k0]) / (thing->vals[2*k1] - thing->vals[2*k0]);