00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "common.h"
00043
00044 #include <math.h>
00045 #include <string.h>
00046 #include "bio.h"
00047
00048 #include "vmath.h"
00049 #include "bu.h"
00050 #include "bn.h"
00051
00052
00053
00054
00055 void
00056 bn_table_free(struct bn_table *tabp)
00057 {
00058 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_free(%p)\n", (void *)tabp);
00059 BN_CK_TABLE(tabp);
00060
00061 tabp->nx = 0;
00062 bu_free(tabp, "struct bn_table");
00063 }
00064
00065
00066
00067
00068 void
00069 bn_tabdata_free(struct bn_tabdata *data)
00070 {
00071 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_free(%p)\n", (void *)data);
00072
00073 BN_CK_TABDATA(data);
00074 BN_CK_TABLE(data->table);
00075
00076 data->ny = 0;
00077 data->table = NULL;
00078 bu_free(data, "struct bn_tabdata");
00079 }
00080
00081
00082
00083
00084 void
00085 bn_ck_table(const struct bn_table *tabp)
00086 {
00087 register size_t i;
00088
00089 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_ck_table(%p)\n", (void *)tabp);
00090
00091 BN_CK_TABLE(tabp);
00092
00093 if (tabp->nx < 2)
00094 bu_bomb("bn_ck_table() less than 2 wavelengths\n");
00095
00096 for (i=0; i < tabp->nx; i++) {
00097 if (tabp->x[i] >= tabp->x[i+1])
00098 bu_bomb("bn_ck_table() wavelengths not in strictly ascending order\n");
00099 }
00100 }
00101
00102
00103
00104
00105
00106
00107
00108 struct bn_table *
00109 bn_table_make_uniform(size_t num, double first, double last) {
00110 struct bn_table *tabp;
00111 fastf_t *fp;
00112 fastf_t delta;
00113 int j;
00114
00115 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_make_uniform( num=%zu, %g, %g )\n", num, first, last);
00116
00117 if (first >= last) bu_bomb("bn_table_make_uniform() first >= last\n");
00118
00119 BN_GET_TABLE(tabp, num);
00120
00121 delta = (last - first) / (double)num;
00122
00123 fp = &tabp->x[0];
00124 for (j = num; j > 0; j--) {
00125 *fp++ = first;
00126 first += delta;
00127 }
00128 tabp->x[num] = last;
00129
00130 return tabp;
00131 }
00132
00133
00134
00135
00136
00137
00138 void
00139 bn_tabdata_add(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2)
00140 {
00141 register int j;
00142 register fastf_t *op;
00143 register const fastf_t *i1, *i2;
00144
00145 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_add(%p, %p, %p)\n", (void *)out, (void *)in1, (void *)in2);
00146
00147 BN_CK_TABDATA(out);
00148 BN_CK_TABDATA(in1);
00149 BN_CK_TABDATA(in2);
00150
00151 if (in1->table != in2->table || in1->table != out->table)
00152 bu_bomb("bn_tabdata_add(): samples drawn from different tables\n");
00153 if (in1->ny != in2->ny || in1->ny != out->ny)
00154 bu_bomb("bn_tabdata_add(): different tabdata lengths?\n");
00155
00156 op = out->y;
00157 i1 = in1->y;
00158 i2 = in2->y;
00159 for (j = in1->ny; j > 0; j--)
00160 *op++ = *i1++ + *i2++;
00161
00162 }
00163
00164
00165
00166
00167
00168
00169 void
00170 bn_tabdata_mul(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2)
00171 {
00172 register int j;
00173 register fastf_t *op;
00174 register const fastf_t *i1, *i2;
00175
00176 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mul(%p, %p, %p)\n", (void *)out, (void *)in1, (void *)in2);
00177
00178 BN_CK_TABDATA(out);
00179 BN_CK_TABDATA(in1);
00180 BN_CK_TABDATA(in2);
00181
00182 if (in1->table != in2->table || in1->table != out->table)
00183 bu_bomb("bn_tabdata_mul(): samples drawn from different tables\n");
00184 if (in1->ny != in2->ny || in1->ny != out->ny)
00185 bu_bomb("bn_tabdata_mul(): different tabdata lengths?\n");
00186
00187 op = out->y;
00188 i1 = in1->y;
00189 i2 = in2->y;
00190 for (j = in1->ny; j > 0; j--)
00191 *op++ = *i1++ * *i2++;
00192
00193 }
00194
00195
00196
00197
00198
00199
00200 void
00201 bn_tabdata_mul3(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, const struct bn_tabdata *in3)
00202 {
00203 register int j;
00204 register fastf_t *op;
00205 register const fastf_t *i1, *i2, *i3;
00206
00207 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mul3(%p, %p, %p, %p)\n", (void *)out, (void *)in1, (void *)in2, (void *)in3);
00208
00209 BN_CK_TABDATA(out);
00210 BN_CK_TABDATA(in1);
00211 BN_CK_TABDATA(in2);
00212 BN_CK_TABDATA(in3);
00213
00214 if (in1->table != in2->table || in1->table != out->table || in1->table != in2->table)
00215 bu_bomb("bn_tabdata_mul(): samples drawn from different tables\n");
00216 if (in1->ny != in2->ny || in1->ny != out->ny)
00217 bu_bomb("bn_tabdata_mul(): different tabdata lengths?\n");
00218
00219 op = out->y;
00220 i1 = in1->y;
00221 i2 = in2->y;
00222 i3 = in3->y;
00223 for (j = in1->ny; j > 0; j--)
00224 *op++ = *i1++ * *i2++ * *i3++;
00225
00226 }
00227
00228
00229
00230
00231
00232
00233
00234
00235 void
00236 bn_tabdata_incr_mul3_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, const struct bn_tabdata *in3, register double scale)
00237 {
00238 register int j;
00239 register fastf_t *op;
00240 register const fastf_t *i1, *i2, *i3;
00241
00242 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_incr_mul3_scale(%p, %p, %p, %p, %g)\n", (void *)out, (void *)in1, (void *)in2, (void *)in3, scale);
00243
00244 BN_CK_TABDATA(out);
00245 BN_CK_TABDATA(in1);
00246 BN_CK_TABDATA(in2);
00247 BN_CK_TABDATA(in3);
00248
00249 if (in1->table != in2->table || in1->table != out->table || in1->table != in3->table)
00250 bu_bomb("bn_tabdata_mul(): samples drawn from different tables\n");
00251 if (in1->ny != in2->ny || in1->ny != out->ny)
00252 bu_bomb("bn_tabdata_mul(): different tabdata lengths?\n");
00253
00254 op = out->y;
00255 i1 = in1->y;
00256 i2 = in2->y;
00257 i3 = in3->y;
00258 for (j = in1->ny; j > 0; j--)
00259 *op++ += *i1++ * *i2++ * *i3++ * scale;
00260 }
00261
00262
00263
00264
00265
00266
00267
00268
00269 void
00270 bn_tabdata_incr_mul2_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, register double scale)
00271 {
00272 register int j;
00273 register fastf_t *op;
00274 register const fastf_t *i1, *i2;
00275
00276 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_incr_mul2_scale(%p, %p, %p, %g)\n", (void *)out, (void *)in1, (void *)in2, scale);
00277
00278 BN_CK_TABDATA(out);
00279 BN_CK_TABDATA(in1);
00280 BN_CK_TABDATA(in2);
00281
00282 if (in1->table != in2->table || in1->table != out->table)
00283 bu_bomb("bn_tabdata_mul(): samples drawn from different tables\n");
00284 if (in1->ny != in2->ny || in1->ny != out->ny)
00285 bu_bomb("bn_tabdata_mul(): different tabdata lengths?\n");
00286
00287 op = out->y;
00288 i1 = in1->y;
00289 i2 = in2->y;
00290 for (j = in1->ny; j > 0; j--)
00291 *op++ += *i1++ * *i2++ * scale;
00292 }
00293
00294
00295
00296
00297
00298
00299 void
00300 bn_tabdata_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale)
00301 {
00302 register int j;
00303 register fastf_t *op;
00304 register const fastf_t *i1;
00305
00306 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_scale(%p, %p, %g)\n", (void *)out, (void *)in1, scale);
00307
00308 BN_CK_TABDATA(out);
00309 BN_CK_TABDATA(in1);
00310
00311 if (in1->table != out->table)
00312 bu_bomb("bn_tabdata_scale(): samples drawn from different tables\n");
00313 if (in1->ny != out->ny)
00314 bu_bomb("bn_tabdata_scale(): different tabdata lengths?\n");
00315
00316 op = out->y;
00317 i1 = in1->y;
00318 for (j = in1->ny; j > 0; j--)
00319 *op++ = *i1++ * scale;
00320
00321 }
00322
00323
00324
00325
00326
00327
00328 void
00329 bn_table_scale(struct bn_table *tabp, register double scale)
00330 {
00331 register int j;
00332 register fastf_t *op;
00333
00334 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_scale(%p, %g)\n", (void *)tabp, scale);
00335
00336 BN_CK_TABLE(tabp);
00337
00338 op = tabp->x;
00339 for (j = tabp->nx+1; j > 0; j--)
00340 *op++ *= scale;
00341
00342 }
00343
00344
00345
00346
00347
00348
00349
00350
00351 void
00352 bn_tabdata_join1(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale, const struct bn_tabdata *in2)
00353 {
00354 register int j;
00355 register fastf_t *op;
00356 register const fastf_t *i1, *i2;
00357
00358 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_join1(%p, %p, %g, %p)\n", (void *)out, (void *)in1, scale, (void *)in2);
00359
00360 BN_CK_TABDATA(out);
00361 BN_CK_TABDATA(in1);
00362 BN_CK_TABDATA(in2);
00363
00364 if (in1->table != out->table)
00365 bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
00366 if (in1->table != in2->table)
00367 bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
00368 if (in1->ny != out->ny)
00369 bu_bomb("bn_tabdata_join1(): different tabdata lengths?\n");
00370
00371 op = out->y;
00372 i1 = in1->y;
00373 i2 = in2->y;
00374 for (j = in1->ny; j > 0; j--)
00375 *op++ = *i1++ + scale * *i2++;
00376
00377 }
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387 void
00388 bn_tabdata_join2(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2, register double scale3, const struct bn_tabdata *in3)
00389 {
00390 register int j;
00391 register fastf_t *op;
00392 register const fastf_t *i1, *i2, *i3;
00393
00394 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_join2(%p, %p, %g, %p, %g, %p)\n", (void *)out, (void *)in1, scale2, (void *)in2, scale3, (void *)in3);
00395
00396 BN_CK_TABDATA(out);
00397 BN_CK_TABDATA(in1);
00398 BN_CK_TABDATA(in2);
00399
00400 if (in1->table != out->table)
00401 bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
00402 if (in1->table != in2->table)
00403 bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
00404 if (in1->table != in3->table)
00405 bu_bomb("bn_tabdata_join1(): samples drawn from different tables\n");
00406 if (in1->ny != out->ny)
00407 bu_bomb("bn_tabdata_join1(): different tabdata lengths?\n");
00408
00409 op = out->y;
00410 i1 = in1->y;
00411 i2 = in2->y;
00412 i3 = in3->y;
00413 for (j = in1->ny; j > 0; j--)
00414 *op++ = *i1++ + scale2 * *i2++ + scale3 * *i3++;
00415
00416 }
00417
00418
00419
00420
00421 void
00422 bn_tabdata_blend2(struct bn_tabdata *out, register double scale1, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2)
00423 {
00424 register int j;
00425 register fastf_t *op;
00426 register const fastf_t *i1, *i2;
00427
00428 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_blend2(%p, %g, %p, %g, %p)\n", (void *)out, scale1, (void *)in1, scale2, (void *)in2);
00429
00430 BN_CK_TABDATA(out);
00431 BN_CK_TABDATA(in1);
00432 BN_CK_TABDATA(in2);
00433
00434 if (in1->table != out->table)
00435 bu_bomb("bn_tabdata_blend2(): samples drawn from different tables\n");
00436 if (in1->table != in2->table)
00437 bu_bomb("bn_tabdata_blend2(): samples drawn from different tables\n");
00438 if (in1->ny != out->ny)
00439 bu_bomb("bn_tabdata_blend2(): different tabdata lengths?\n");
00440
00441 op = out->y;
00442 i1 = in1->y;
00443 i2 = in2->y;
00444 for (j = in1->ny; j > 0; j--)
00445 *op++ = scale1 * *i1++ + scale2 * *i2++;
00446
00447 }
00448
00449
00450
00451
00452 void
00453 bn_tabdata_blend3(struct bn_tabdata *out, register double scale1, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2, register double scale3, const struct bn_tabdata *in3)
00454 {
00455 register int j;
00456 register fastf_t *op;
00457 register const fastf_t *i1, *i2, *i3;
00458
00459 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_blend3(%p, %g, %p, %g, %p, %g, %p)\n", (void *)out, scale1, (void *)in1, scale2, (void *)in2, scale3, (void *)in3);
00460
00461 BN_CK_TABDATA(out);
00462 BN_CK_TABDATA(in1);
00463 BN_CK_TABDATA(in2);
00464 BN_CK_TABDATA(in3);
00465
00466 if (in1->table != out->table)
00467 bu_bomb("bn_tabdata_blend3(): samples drawn from different tables\n");
00468 if (in1->table != in2->table)
00469 bu_bomb("bn_tabdata_blend3(): samples drawn from different tables\n");
00470 if (in1->table != in3->table)
00471 bu_bomb("bn_tabdata_blend3(): samples drawn from different tables\n");
00472 if (in1->ny != out->ny)
00473 bu_bomb("bn_tabdata_blend3(): different tabdata lengths?\n");
00474
00475 op = out->y;
00476 i1 = in1->y;
00477 i2 = in2->y;
00478 i3 = in3->y;
00479 for (j = in1->ny; j > 0; j--)
00480 *op++ = scale1 * *i1++ + scale2 * *i2++ + scale3 * *i3++;
00481
00482 }
00483
00484
00485
00486
00487
00488
00489
00490
00491 double
00492 bn_tabdata_area1(const struct bn_tabdata *in)
00493 {
00494 register fastf_t area;
00495 register const fastf_t *ip;
00496 register int j;
00497
00498 BN_CK_TABDATA(in);
00499
00500 area = 0;
00501 ip = in->y;
00502 for (j = in->ny; j > 0; j--)
00503 area += *ip++;
00504
00505 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_area(%p) = %g\n", (void *)in, area);
00506
00507 return area;
00508 }
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518 double
00519 bn_tabdata_area2(const struct bn_tabdata *in)
00520 {
00521 const struct bn_table *tabp;
00522 register fastf_t area;
00523 fastf_t width;
00524 register int j;
00525
00526 BN_CK_TABDATA(in);
00527 tabp = in->table;
00528 BN_CK_TABLE(tabp);
00529
00530 area = 0;
00531 for (j = in->ny-1; j >= 0; j--) {
00532 width = tabp->x[j+1] - tabp->x[j];
00533 area += in->y[j] * width;
00534 }
00535
00536 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_area2(%p) = %g\n", (void *)in, area);
00537 return area;
00538 }
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549 double
00550 bn_tabdata_mul_area1(const struct bn_tabdata *in1, const struct bn_tabdata *in2)
00551 {
00552 register fastf_t area;
00553 register const fastf_t *i1, *i2;
00554 register int j;
00555
00556 BN_CK_TABDATA(in1);
00557 BN_CK_TABDATA(in2);
00558
00559 area = 0;
00560 i1 = in1->y;
00561 i2 = in2->y;
00562 for (j = in1->ny; j > 0; j--)
00563 area += *i1++ * *i2++;
00564
00565 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mul_area1(%p, %p) = %g\n", (void *)in1, (void *)in2, area);
00566 return area;
00567 }
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577 double
00578 bn_tabdata_mul_area2(const struct bn_tabdata *in1, const struct bn_tabdata *in2)
00579 {
00580 const struct bn_table *tabp;
00581 register fastf_t area;
00582 fastf_t width;
00583 register int j;
00584
00585 BN_CK_TABDATA(in1);
00586 BN_CK_TABDATA(in2);
00587 tabp = in1->table;
00588 BN_CK_TABLE(tabp);
00589
00590 area = 0;
00591 for (j = in1->ny-1; j >= 0; j--) {
00592 width = tabp->x[j+1] - tabp->x[j];
00593 area += in1->y[j] * in2->y[j] * width;
00594 }
00595
00596 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mul_area2(%p, %p) = %g\n", (void *)in1, (void *)in2, area);
00597 return area;
00598 }
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615 long
00616 bn_table_find_x(const struct bn_table *tabp, double xval)
00617 {
00618 register int i;
00619
00620 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_find_x(%p, %g)\n", (void *)tabp, xval);
00621 BN_CK_TABLE(tabp);
00622
00623 if (xval > tabp->x[tabp->nx]) return -2;
00624 if (xval >= tabp->x[tabp->nx-1]) return tabp->nx-1;
00625
00626
00627 for (i = tabp->nx-2; i >=0; i--) {
00628 if (xval >= tabp->x[i]) return i;
00629 }
00630
00631 return -1;
00632 }
00633
00634
00635
00636
00637
00638
00639
00640
00641 fastf_t
00642 bn_table_lin_interp(const struct bn_tabdata *samp, register double wl)
00643 {
00644 const struct bn_table *tabp;
00645 size_t i;
00646 int idx;
00647 fastf_t fract;
00648 fastf_t ret;
00649
00650 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_lin_interp(%p, %g)\n", (void *)samp, wl);
00651
00652 BN_CK_TABDATA(samp);
00653 tabp = samp->table;
00654 BN_CK_TABLE(tabp);
00655
00656 idx = bn_table_find_x(tabp, wl);
00657 if (idx < 0) {
00658 if (bu_debug&BU_DEBUG_TABDATA)
00659 bu_log("bn_table_lin_interp(%g) out of range %g to %g\n", wl, tabp->x[0], tabp->x[tabp->nx]);
00660 return 0;
00661 }
00662 i = (size_t)idx;
00663
00664 if (wl < tabp->x[i] || wl >= tabp->x[i+1]) {
00665 bu_log("bn_table_lin_interp(%g) assertion1 failed at %g\n", wl, tabp->x[i]);
00666 bu_bomb("bn_table_lin_interp() assertion1 failed\n");
00667 }
00668
00669 if (i >= tabp->nx-2) {
00670
00671 if (bu_debug&BU_DEBUG_TABDATA)bu_log("bn_table_lin_interp(%g)=%g off end of range %g to %g\n", wl, samp->y[tabp->nx-1], tabp->x[0], tabp->x[tabp->nx]);
00672 return samp->y[tabp->nx-1];
00673 }
00674
00675
00676 fract = (wl - tabp->x[i]) / (tabp->x[i+1] - tabp->x[i]);
00677 if (fract < 0 || fract > 1) bu_bomb("bn_table_lin_interp() assertion2 failed\n");
00678 ret = (1-fract) * samp->y[i] + fract * samp->y[i+1];
00679 if (bu_debug&BU_DEBUG_TABDATA)bu_log("bn_table_lin_interp(%g)=%g in range %g to %g\n",
00680 wl, ret, tabp->x[i], tabp->x[i+1]);
00681 return ret;
00682 }
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695 struct bn_tabdata *
00696 bn_tabdata_resample_max(const struct bn_table *newtable, const struct bn_tabdata *olddata) {
00697 const struct bn_table *oldtable;
00698 struct bn_tabdata *newsamp;
00699 size_t i;
00700 int j, k;
00701
00702 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_resample_max(%p, %p)\n", (void *)newtable, (void *)olddata);
00703
00704 BN_CK_TABLE(newtable);
00705 BN_CK_TABDATA(olddata);
00706 oldtable = olddata->table;
00707 BN_CK_TABLE(oldtable);
00708
00709 if (oldtable == newtable) bu_log("bn_tabdata_resample_max() NOTICE old and new bn_table structs are the same\n");
00710
00711 BN_GET_TABDATA(newsamp, newtable);
00712
00713 for (i = 0; i < newtable->nx; i++) {
00714
00715
00716
00717
00718 j = bn_table_find_x(oldtable, newtable->x[i]);
00719 k = bn_table_find_x(oldtable, newtable->x[i+1]);
00720 if (k == -1) {
00721
00722 newsamp->y[i] = 0;
00723 continue;
00724 }
00725 if (j == -2) {
00726
00727 newsamp->y[i] = 0;
00728 continue;
00729 }
00730
00731 if (j == k && j > 0) {
00732 register fastf_t tmp;
00733
00734
00735
00736
00737
00738
00739 newsamp->y[i] = bn_table_lin_interp(olddata, newtable->x[i]);
00740 tmp = bn_table_lin_interp(olddata, newtable->x[i+1]);
00741 if (tmp > newsamp->y[i]) newsamp->y[i] = tmp;
00742 } else {
00743 register fastf_t tmp, n;
00744 register int s;
00745
00746
00747
00748
00749
00750
00751 n = bn_table_lin_interp(olddata, newtable->x[i]);
00752 tmp = bn_table_lin_interp(olddata, newtable->x[i+1]);
00753 if (tmp > n) n = tmp;
00754 for (s = j+1; s <= k; s++) {
00755 if ((tmp = olddata->y[s]) > n)
00756 n = tmp;
00757 }
00758 newsamp->y[i] = n;
00759 }
00760 }
00761 return newsamp;
00762 }
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775 struct bn_tabdata *
00776 bn_tabdata_resample_avg(const struct bn_table *newtable, const struct bn_tabdata *olddata) {
00777 const struct bn_table *oldtable;
00778 struct bn_tabdata *newsamp;
00779 size_t i;
00780 int j, k;
00781
00782 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_resample_avg(%p, %p)\n", (void *)newtable, (void *)olddata);
00783
00784 BN_CK_TABLE(newtable);
00785 BN_CK_TABDATA(olddata);
00786 oldtable = olddata->table;
00787 BN_CK_TABLE(oldtable);
00788
00789 if (oldtable == newtable) bu_log("bn_tabdata_resample_avg() NOTICE old and new bn_table structs are the same\n");
00790
00791 BN_GET_TABDATA(newsamp, newtable);
00792
00793 for (i = 0; i < newtable->nx; i++) {
00794
00795
00796
00797
00798 j = bn_table_find_x(oldtable, newtable->x[i]);
00799 k = bn_table_find_x(oldtable, newtable->x[i+1]);
00800
00801 if (j < 0 || k < 0 || j == k) {
00802
00803
00804
00805
00806
00807
00808 newsamp->y[i] = 0.5 * (
00809 bn_table_lin_interp(olddata, newtable->x[i]) +
00810 bn_table_lin_interp(olddata, newtable->x[i+1]));
00811 } else {
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821 fastf_t wsum;
00822 fastf_t a, b;
00823 int s;
00824
00825
00826 a = bn_table_lin_interp(olddata, newtable->x[i]);
00827 b = olddata->y[j+1];
00828 wsum = 0.5 * (a+b) * (oldtable->x[j+1] - newtable->x[i]);
00829
00830
00831 for (s = j+1; s < k; s++) {
00832 a = olddata->y[s];
00833 b = olddata->y[s+1];
00834 wsum += 0.5 * (a+b) * (oldtable->x[s+1] - oldtable->x[s]);
00835 }
00836
00837
00838 a = olddata->y[k];
00839 b = bn_table_lin_interp(olddata, newtable->x[i+1]);
00840 wsum += 0.5 * (a+b) * (newtable->x[i+1] - oldtable->x[k]);
00841
00842
00843 newsamp->y[i] =
00844 wsum / (newtable->x[i+1] - newtable->x[i]);
00845 }
00846 }
00847 return newsamp;
00848 }
00849
00850
00851
00852
00853
00854
00855
00856
00857 int
00858 bn_table_write(const char *filename, const struct bn_table *tabp)
00859 {
00860 FILE *fp;
00861 size_t j;
00862
00863 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_write(%s, %p)\n", filename, (void *)tabp);
00864
00865 BN_CK_TABLE(tabp);
00866
00867 bu_semaphore_acquire(BU_SEM_SYSCALL);
00868 fp = fopen(filename, "wb");
00869 bu_semaphore_release(BU_SEM_SYSCALL);
00870
00871 if (fp == NULL) {
00872 perror(filename);
00873 bu_log("bn_table_write(%s, %p) FAILED\n", filename, (void *)tabp);
00874 return -1;
00875 }
00876
00877 bu_semaphore_acquire(BU_SEM_SYSCALL);
00878 fprintf(fp, " %lu sample starts, and one end.\n", (long unsigned int)tabp->nx);
00879 for (j=0; j <= tabp->nx; j++) {
00880 fprintf(fp, "%g\n", tabp->x[j]);
00881 }
00882 fclose(fp);
00883 bu_semaphore_release(BU_SEM_SYSCALL);
00884 return 0;
00885 }
00886
00887
00888
00889
00890
00891
00892
00893
00894 struct bn_table *
00895 bn_table_read(const char *filename) {
00896 int ret;
00897 struct bn_table *tabp;
00898 struct bu_vls line = BU_VLS_INIT_ZERO;
00899 FILE *fp;
00900 size_t nw;
00901 size_t j;
00902
00903 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_read(%s)\n", filename);
00904
00905 bu_semaphore_acquire(BU_SEM_SYSCALL);
00906 fp = fopen(filename, "rb");
00907 bu_semaphore_release(BU_SEM_SYSCALL);
00908
00909 if (fp == NULL) {
00910 perror(filename);
00911 bu_log("bn_table_read(%s) FAILED\n", filename);
00912 return NULL;
00913 }
00914
00915 if (bu_vls_gets(&line, fp) < 0) {
00916 perror(filename);
00917 bu_log("Failed to read line\n");
00918 fclose(fp);
00919 return NULL;
00920 }
00921 nw = 0;
00922
00923
00924
00925 sscanf(bu_vls_addr(&line), "%lu", (long unsigned int *)(&nw));
00926 bu_vls_free(&line);
00927
00928 if (nw <= 0) bu_bomb("bn_table_read() bad nw value\n");
00929
00930 BN_GET_TABLE(tabp, nw);
00931
00932 bu_semaphore_acquire(BU_SEM_SYSCALL);
00933 for (j=0; j <= tabp->nx; j++) {
00934
00935 ret = fscanf(fp, "%lf", &tabp->x[j]);
00936 if (ret != 1) {
00937 bu_log("bn_table_read(%s) READ FAILURE. Abort\n", filename);
00938 break;
00939 }
00940 }
00941 fclose(fp);
00942 bu_semaphore_release(BU_SEM_SYSCALL);
00943
00944 bn_ck_table(tabp);
00945
00946 return tabp;
00947 }
00948
00949
00950
00951
00952 void
00953 bn_pr_table(const char *title, const struct bn_table *tabp)
00954 {
00955 size_t j;
00956
00957 bu_log("%s\n", title);
00958 BN_CK_TABLE(tabp);
00959
00960 for (j=0; j <= tabp->nx; j++) {
00961 bu_log("%3zu: %g\n", j, tabp->x[j]);
00962 }
00963 }
00964
00965
00966
00967
00968 void
00969 bn_pr_tabdata(const char *title, const struct bn_tabdata *data)
00970 {
00971 size_t j;
00972
00973 bu_log("%s: ", title);
00974 BN_CK_TABDATA(data);
00975
00976 for (j=0; j < data->ny; j++) {
00977 bu_log("%g, ", data->y[j]);
00978 }
00979 bu_log("\n");
00980 }
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992 int
00993 bn_print_table_and_tabdata(const char *filename, const struct bn_tabdata *data)
00994 {
00995 FILE *fp;
00996 const struct bn_table *tabp;
00997 size_t j;
00998
00999 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_print_table_and_tabdata(%s, %p)\n", filename, (void *)data);
01000
01001 BN_CK_TABDATA(data);
01002 tabp = data->table;
01003 BN_CK_TABLE(tabp);
01004
01005 bu_semaphore_acquire(BU_SEM_SYSCALL);
01006 fp = fopen(filename, "wb");
01007 bu_semaphore_release(BU_SEM_SYSCALL);
01008
01009 if (fp == NULL) {
01010 perror(filename);
01011 bu_log("bn_print_table_and_tabdata(%s, %p) FAILED\n", filename, (void *)data);
01012 return -1;
01013 }
01014
01015 bu_semaphore_acquire(BU_SEM_SYSCALL);
01016 for (j=0; j < tabp->nx; j++) {
01017 fprintf(fp, "%g %g\n", tabp->x[j], data->y[j]);
01018 }
01019 fprintf(fp, "%g (novalue)\n", tabp->x[tabp->nx]);
01020 fclose(fp);
01021 bu_semaphore_release(BU_SEM_SYSCALL);
01022 return 0;
01023 }
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036 struct bn_tabdata *
01037 bn_read_table_and_tabdata(const char *filename) {
01038 struct bn_table *tabp;
01039 struct bn_tabdata *data;
01040 FILE *fp;
01041 char buf[128];
01042 size_t count = 0;
01043 size_t i;
01044
01045 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_read_table_and_tabdata(%s)\n", filename);
01046
01047 bu_semaphore_acquire(BU_SEM_SYSCALL);
01048 fp = fopen(filename, "rb");
01049 bu_semaphore_release(BU_SEM_SYSCALL);
01050
01051 if (fp == NULL) {
01052 perror(filename);
01053 bu_log("bn_read_table_and_tabdata(%s) FAILED. Couldn't open file.\n", filename);
01054 return NULL;
01055 }
01056
01057
01058 bu_semaphore_acquire(BU_SEM_SYSCALL);
01059 for (;;) {
01060 if (bu_fgets(buf, sizeof(buf), fp) == NULL) break;
01061 count++;
01062 }
01063 fclose(fp);
01064 bu_semaphore_release(BU_SEM_SYSCALL);
01065
01066 if (count < 2) {
01067 perror(filename);
01068 bu_log("bn_read_table_and_tabdata(%s) FAILED. Expected at least 2 lines in file.\n", filename);
01069 return NULL;
01070 }
01071
01072
01073 BN_GET_TABLE(tabp, count);
01074 BN_GET_TABDATA(data, tabp);
01075
01076
01077 bu_semaphore_acquire(BU_SEM_SYSCALL);
01078 fp = fopen(filename, "rb");
01079 for (i=0; i < count; i++) {
01080 buf[0] = '\0';
01081 if (bu_fgets(buf, sizeof(buf), fp) == NULL) {
01082 bu_log("bn_read_table_and_tabdata(%s) unexpected EOF on line %zu\n", filename, i);
01083 break;
01084 }
01085 sscanf(buf, "%lf %lf", &tabp->x[i], &data->y[i]);
01086 }
01087 fclose(fp);
01088 bu_semaphore_release(BU_SEM_SYSCALL);
01089
01090
01091 tabp->x[count] = 2 * tabp->x[count-1] - tabp->x[count-2];
01092
01093 bn_ck_table(tabp);
01094
01095 return data;
01096 }
01097
01098
01099
01100
01101 struct bn_tabdata *
01102 bn_tabdata_binary_read(const char *filename, size_t num, const struct bn_table *tabp) {
01103 struct bn_tabdata *data;
01104 char *cp;
01105 int nbytes;
01106 int len;
01107 int got;
01108 int fd;
01109 long i;
01110
01111 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_binary_read(%s, num=%zu, %p)\n", filename, num, (void *)tabp);
01112
01113 BN_CK_TABLE(tabp);
01114
01115 nbytes = BN_SIZEOF_TABDATA(tabp);
01116 len = num * nbytes;
01117
01118 bu_semaphore_acquire(BU_SEM_SYSCALL);
01119 fd = open(filename, 0);
01120 bu_semaphore_release(BU_SEM_SYSCALL);
01121 if (fd <= 0) {
01122 perror(filename);
01123 bu_log("bn_tabdata_binary_read open failed on \"%s\"\n", filename);
01124 return (struct bn_tabdata *)NULL;
01125 }
01126
01127
01128 data = (struct bn_tabdata *)bu_malloc(len+8, "bn_tabdata[]");
01129
01130 bu_semaphore_acquire(BU_SEM_SYSCALL);
01131 got = read(fd, (char *)data, len);
01132 bu_semaphore_release(BU_SEM_SYSCALL);
01133 if (got != len) {
01134 if (got < 0) {
01135 perror(filename);
01136 bu_log("bn_tabdata_binary_read read error on \"%s\"\n", filename);
01137 } else {
01138 bu_log("bn_tabdata_binary_read(%s) expected %d got %d\n", filename, len, got);
01139 }
01140 bu_free(data, "bn_tabdata[]");
01141 bu_semaphore_acquire(BU_SEM_SYSCALL);
01142 close(fd);
01143 bu_semaphore_release(BU_SEM_SYSCALL);
01144 return (struct bn_tabdata *)NULL;
01145 }
01146 bu_semaphore_acquire(BU_SEM_SYSCALL);
01147 close(fd);
01148 bu_semaphore_release(BU_SEM_SYSCALL);
01149
01150
01151 cp = (char *)data;
01152 for (i = num-1; i >= 0; i--, cp += nbytes) {
01153 register struct bn_tabdata *sp;
01154 sp = (struct bn_tabdata *)cp;
01155 BN_CK_TABDATA(sp);
01156 sp->table = tabp;
01157 }
01158
01159 return data;
01160 }
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170 struct bn_tabdata *
01171 bn_tabdata_malloc_array(const struct bn_table *tabp, size_t num) {
01172 struct bn_tabdata *data;
01173 char *cp;
01174 size_t i;
01175 size_t nw;
01176 size_t nbytes;
01177
01178 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_malloc_array(%p, num=%zu)\n", (void *)tabp, num);
01179
01180 BN_CK_TABLE(tabp);
01181 nw = tabp->nx;
01182 nbytes = BN_SIZEOF_TABDATA(tabp);
01183
01184 data = (struct bn_tabdata *)bu_calloc(num,
01185 nbytes, "struct bn_tabdata[]");
01186
01187 cp = (char *)data;
01188 for (i = 0; i < num; i++) {
01189 register struct bn_tabdata *sp;
01190
01191 sp = (struct bn_tabdata *)cp;
01192 sp->magic = BN_TABDATA_MAGIC;
01193 sp->ny = nw;
01194 sp->table = tabp;
01195 cp += nbytes;
01196 }
01197 return data;
01198 }
01199
01200
01201
01202
01203 void
01204 bn_tabdata_copy(struct bn_tabdata *out, const struct bn_tabdata *in)
01205 {
01206 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_copy(%p, %p)\n", (void *)out, (void *)in);
01207
01208 BN_CK_TABDATA(out);
01209 BN_CK_TABDATA(in);
01210
01211 if (in->table != out->table)
01212 bu_bomb("bn_tabdata_copy(): samples drawn from different tables\n");
01213 if (in->ny != out->ny)
01214 bu_bomb("bn_tabdata_copy(): different tabdata lengths?\n");
01215
01216 memcpy((char *)out->y, (const char *)in->y, BN_SIZEOF_TABDATA_Y(in));
01217 }
01218
01219
01220
01221
01222 struct bn_tabdata *
01223 bn_tabdata_dup(const struct bn_tabdata *in) {
01224 struct bn_tabdata *data;
01225
01226 BN_CK_TABDATA(in);
01227 BN_GET_TABDATA(data, in->table);
01228
01229 memcpy((char *)data->y, (const char *)in->y, BN_SIZEOF_TABDATA_Y(in));
01230
01231 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_dup(%p) = %p\n", (void *)in, (void *)data);
01232 return data;
01233 }
01234
01235
01236
01237
01238
01239
01240
01241 struct bn_tabdata *
01242 bn_tabdata_get_constval(double val, const struct bn_table *tabp) {
01243 struct bn_tabdata *data;
01244 int todo;
01245 register fastf_t *op;
01246
01247 BN_CK_TABLE(tabp);
01248 BN_GET_TABDATA(data, tabp);
01249
01250 op = data->y;
01251 for (todo = data->ny-1; todo >= 0; todo--)
01252 *op++ = val;
01253
01254 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_get_constval(val=%g, %p)=%p\n", val, (void *)tabp, (void *)data);
01255
01256 return data;
01257 }
01258
01259
01260
01261
01262
01263
01264 void
01265 bn_tabdata_constval(struct bn_tabdata *data, double val)
01266 {
01267 int todo;
01268 register fastf_t *op;
01269
01270 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_constval(%p, val=%g)\n", (void *)data, val);
01271
01272 BN_CK_TABDATA(data);
01273
01274 op = data->y;
01275 for (todo = data->ny-1; todo >= 0; todo--)
01276 *op++ = val;
01277 }
01278
01279
01280
01281
01282
01283
01284
01285
01286 void
01287 bn_tabdata_to_tcl(struct bu_vls *vp, const struct bn_tabdata *data)
01288 {
01289 const struct bn_table *tabp;
01290 register size_t i;
01291 register fastf_t minval = MAX_FASTF, maxval = -MAX_FASTF;
01292
01293 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_to_tcl(%p, %p)\n", (void *)vp, (void *)data);
01294
01295 BU_CK_VLS(vp);
01296 BN_CK_TABDATA(data);
01297 tabp = data->table;
01298 BN_CK_TABLE(tabp);
01299
01300 bu_vls_strcat(vp, "x {");
01301 for (i=0; i < tabp->nx; i++) {
01302 bu_vls_printf(vp, "%g ", tabp->x[i]);
01303 }
01304 bu_vls_strcat(vp, "} y {");
01305 for (i=0; i < data->ny; i++) {
01306 register fastf_t val = data->y[i];
01307 bu_vls_printf(vp, "%g ", val);
01308 if (val < minval) minval = val;
01309 if (val > maxval) maxval = val;
01310 }
01311 bu_vls_printf(vp, "} nx %zu ymin %g ymax %g",
01312 tabp->nx, minval, maxval);
01313 }
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325 struct bn_tabdata *
01326 bn_tabdata_from_array(const double *array) {
01327 register const double *dp;
01328 size_t len = 0;
01329 struct bn_table *tabp;
01330 struct bn_tabdata *data;
01331 register size_t i;
01332
01333
01334 for (dp = array; *dp > 0; dp += 2)
01335 ;
01336 len = (dp - array) >> 1;
01337
01338
01339 BN_GET_TABLE(tabp, len);
01340 for (i = 0; i < len; i++) {
01341 tabp->x[i] = array[i<<1];
01342 }
01343 tabp->x[len] = tabp->x[len-1] + 1;
01344
01345
01346 BN_GET_TABDATA(data, tabp);
01347 for (i = 0; i < len-1; i++) {
01348 data->y[i] = array[(i<<1)+1];
01349 }
01350
01351 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_from_array(%p) = %p\n", (void *)array, (void *)data);
01352 return data;
01353 }
01354
01355
01356
01357
01358
01359
01360
01361 void
01362 bn_tabdata_freq_shift(struct bn_tabdata *out, const struct bn_tabdata *in, double offset)
01363 {
01364 const struct bn_table *tabp;
01365 register size_t i;
01366
01367 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_freq_shift(%p, %p, offset=%g)\n", (void *)out, (void *)in, offset);
01368
01369 BN_CK_TABDATA(out);
01370 BN_CK_TABDATA(in);
01371 tabp = in->table;
01372
01373 if (tabp != out->table)
01374 bu_bomb("bn_tabdata_freq_shift(): samples drawn from different tables\n");
01375 if (in->ny != out->ny)
01376 bu_bomb("bn_tabdata_freq_shift(): different tabdata lengths?\n");
01377
01378 for (i=0; i < out->ny; i++) {
01379 out->y[i] = bn_table_lin_interp(in, tabp->x[i]+offset);
01380 }
01381 }
01382
01383
01384
01385
01386
01387
01388 int
01389 bn_table_interval_num_samples(const struct bn_table *tabp, double low, double hi)
01390 {
01391 register size_t i;
01392 register size_t count = 0;
01393
01394 BN_CK_TABLE(tabp);
01395
01396 for (i=0; i < tabp->nx-1; i++) {
01397 if (tabp->x[i] >= low && tabp->x[i] <= hi) count++;
01398 }
01399
01400 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_interval_num_samples(%p, low=%g, hi=%g) = %zu\n", (void *)tabp, low, hi, count);
01401 return count;
01402 }
01403
01404
01405
01406
01407
01408
01409
01410
01411 int
01412 bn_table_delete_sample_pts(struct bn_table *tabp, unsigned int i, unsigned int j)
01413 {
01414 size_t tokill;
01415 unsigned int k;
01416
01417 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_delete_samples(%p, %d, %d)\n", (void *)tabp, i, j);
01418
01419 BN_CK_TABLE(tabp);
01420
01421 if (i >= tabp->nx || j >= tabp->nx)
01422 bu_bomb("bn_table_delete_sample_pts() index out of range\n");
01423
01424 tokill = j - i + 1;
01425 if (tokill < 1) bu_bomb("bn_table_delete_sample_pts(): nothing to delete\n");
01426 if (tokill >= tabp->nx) bu_bomb("bn_table_delete_sample_pts(): you can't kill 'em all!\n");
01427
01428 tabp->nx -= tokill;
01429
01430 for (k = i; k < tabp->nx; k++) {
01431 tabp->x[k] = tabp->x[k+tokill];
01432 }
01433 return tokill;
01434 }
01435
01436
01437
01438
01439
01440
01441
01442 struct bn_table *
01443 bn_table_merge2(const struct bn_table *a, const struct bn_table *b) {
01444 struct bn_table *newtab;
01445 register size_t i, j, k;
01446
01447 BN_CK_TABLE(a);
01448 BN_CK_TABLE(b);
01449
01450 BN_GET_TABLE(newtab, a->nx + b->nx + 2);
01451
01452 i = j = 0;
01453 k = 0;
01454 while (i <= a->nx || j <= b->nx) {
01455 if (i > a->nx) {
01456 while (j <= b->nx)
01457 newtab->x[k++] = b->x[j++];
01458 break;
01459 }
01460 if (j > b->nx) {
01461 while (i <= a->nx)
01462 newtab->x[k++] = a->x[i++];
01463 break;
01464 }
01465
01466 if (ZERO(a->x[i] - b->x[j])) {
01467 newtab->x[k++] = a->x[i++];
01468 j++;
01469 continue;
01470 }
01471 if (a->x[i] <= b->x[j]) {
01472 newtab->x[k++] = a->x[i++];
01473 } else {
01474 newtab->x[k++] = b->x[j++];
01475 }
01476 }
01477 if (k > newtab->nx) bu_bomb("bn_table_merge2() assertion failed, k>nx?\n");
01478 newtab->nx = k-1;
01479
01480 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_table_merge2(%p, %p) = %p\n", (void *)a, (void *)b, (void *)newtab);
01481 return newtab;
01482 }
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496 struct bn_tabdata *
01497 bn_tabdata_mk_linear_filter(const struct bn_table *spectrum, double lower_wavelen, double upper_wavelen) {
01498 struct bn_tabdata *filt;
01499 int first;
01500 int last;
01501 fastf_t filt_range;
01502 fastf_t cell_range;
01503 fastf_t frac;
01504 size_t i;
01505
01506 BN_CK_TABLE(spectrum);
01507
01508 if (bu_debug&BU_DEBUG_TABDATA) bu_log("bn_tabdata_mk_linear_filter(%p, low=%g, up=%g)\n", (void *)spectrum, lower_wavelen, upper_wavelen);
01509
01510 if (lower_wavelen < spectrum->x[0])
01511 if (upper_wavelen > spectrum->x[spectrum->nx])
01512 bu_log("bn_tabdata_mk_linear_filter() warning, upper_wavelen %g > hightest sampled wavelen %g\n",
01513 upper_wavelen, spectrum->x[spectrum->nx]);
01514
01515
01516 first = bn_table_find_x(spectrum, lower_wavelen);
01517 if (first == -1) {
01518 first = 0;
01519 bu_log("bn_tabdata_mk_linear_filter() warning, lower_wavelen %g < lowest sampled wavelen %g\n",
01520 lower_wavelen, spectrum->x[0]);
01521 } else if (first <= -2) {
01522 bu_log("bn_tabdata_mk_linear_filter() ERROR, lower_wavelen %g > highest sampled wavelen %g\n",
01523 lower_wavelen, spectrum->x[spectrum->nx]);
01524 return NULL;
01525 }
01526
01527
01528 last = bn_table_find_x(spectrum, upper_wavelen);
01529 if (last == -1) {
01530 bu_log("bn_tabdata_mk_linear_filter() ERROR, upper_wavelen %g < lowest sampled wavelen %g\n",
01531 upper_wavelen, spectrum->x[0]);
01532 return NULL;
01533 } else if (last <= -2) {
01534 last = spectrum->nx-1;
01535 bu_log("bn_tabdata_mk_linear_filter() warning, upper_wavelen %g > highest sampled wavelen %g\n",
01536 upper_wavelen, spectrum->x[spectrum->nx]);
01537 }
01538
01539
01540 BN_GET_TABDATA(filt, spectrum);
01541
01542
01543 if (first == last) {
01544 filt_range = upper_wavelen - lower_wavelen;
01545 cell_range = spectrum->x[first+1] - spectrum->x[first];
01546 frac = filt_range / cell_range;
01547
01548
01549 BU_ASSERT((frac >= 0.0) && (frac <= 1.0));
01550
01551 filt->y[first] = frac;
01552 return filt;
01553 }
01554
01555
01556 filt_range = spectrum->x[first+1] - lower_wavelen;
01557 cell_range = spectrum->x[first+1] - spectrum->x[first];
01558 frac = filt_range / cell_range;
01559 if (frac > 1) frac = 1;
01560 BU_ASSERT((frac >= 0.0) && (frac <= 1.0));
01561 filt->y[first] = frac;
01562
01563 filt_range = upper_wavelen - spectrum->x[last];
01564 cell_range = spectrum->x[last+1] - spectrum->x[last];
01565 frac = filt_range / cell_range;
01566 if (frac > 1) frac = 1;
01567 BU_ASSERT((frac >= 0.0) && (frac <= 1.0));
01568 filt->y[last] = frac;
01569
01570
01571 for (i = first+1; i < (size_t)last; i++)
01572 filt->y[i] = 1.0;
01573
01574 return filt;
01575 }
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585