| 393 |
|
return((rbfv[0] != NULL) + (rbfv[1] != NULL)); |
| 394 |
|
} |
| 395 |
|
|
| 396 |
+ |
/* Advect and allocate new RBF along edge (internal call) */ |
| 397 |
+ |
RBFNODE * |
| 398 |
+ |
e_advect_rbf(const MIGRATION *mig, const FVECT invec, int lobe_lim) |
| 399 |
+ |
{ |
| 400 |
+ |
double cthresh = FTINY; |
| 401 |
+ |
RBFNODE *rbf; |
| 402 |
+ |
int n, i, j; |
| 403 |
+ |
double t, full_dist; |
| 404 |
+ |
/* get relative position */ |
| 405 |
+ |
t = Acos(DOT(invec, mig->rbfv[0]->invec)); |
| 406 |
+ |
if (t < M_PI/grid_res) { /* near first DSF */ |
| 407 |
+ |
n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1); |
| 408 |
+ |
rbf = (RBFNODE *)malloc(n); |
| 409 |
+ |
if (rbf == NULL) |
| 410 |
+ |
goto memerr; |
| 411 |
+ |
memcpy(rbf, mig->rbfv[0], n); /* just duplicate */ |
| 412 |
+ |
rbf->next = NULL; rbf->ejl = NULL; |
| 413 |
+ |
return(rbf); |
| 414 |
+ |
} |
| 415 |
+ |
full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec)); |
| 416 |
+ |
if (t > full_dist-M_PI/grid_res) { /* near second DSF */ |
| 417 |
+ |
n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1); |
| 418 |
+ |
rbf = (RBFNODE *)malloc(n); |
| 419 |
+ |
if (rbf == NULL) |
| 420 |
+ |
goto memerr; |
| 421 |
+ |
memcpy(rbf, mig->rbfv[1], n); /* just duplicate */ |
| 422 |
+ |
rbf->next = NULL; rbf->ejl = NULL; |
| 423 |
+ |
return(rbf); |
| 424 |
+ |
} |
| 425 |
+ |
t /= full_dist; |
| 426 |
+ |
tryagain: |
| 427 |
+ |
n = 0; /* count migrating particles */ |
| 428 |
+ |
for (i = 0; i < mtx_nrows(mig); i++) |
| 429 |
+ |
for (j = 0; j < mtx_ncols(mig); j++) |
| 430 |
+ |
n += (mtx_coef(mig,i,j) > cthresh); |
| 431 |
+ |
/* are we over our limit? */ |
| 432 |
+ |
if ((lobe_lim > 0) & (n > lobe_lim)) { |
| 433 |
+ |
cthresh = cthresh*2. + 10.*FTINY; |
| 434 |
+ |
goto tryagain; |
| 435 |
+ |
} |
| 436 |
+ |
#ifdef DEBUG |
| 437 |
+ |
fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n", |
| 438 |
+ |
mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n); |
| 439 |
+ |
#endif |
| 440 |
+ |
rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1)); |
| 441 |
+ |
if (rbf == NULL) |
| 442 |
+ |
goto memerr; |
| 443 |
+ |
rbf->next = NULL; rbf->ejl = NULL; |
| 444 |
+ |
VCOPY(rbf->invec, invec); |
| 445 |
+ |
rbf->nrbf = n; |
| 446 |
+ |
rbf->vtotal = 1.-t + t*mig->rbfv[1]->vtotal/mig->rbfv[0]->vtotal; |
| 447 |
+ |
n = 0; /* advect RBF lobes */ |
| 448 |
+ |
for (i = 0; i < mtx_nrows(mig); i++) { |
| 449 |
+ |
const RBFVAL *rbf0i = &mig->rbfv[0]->rbfa[i]; |
| 450 |
+ |
const float peak0 = rbf0i->peak; |
| 451 |
+ |
const double rad0 = R2ANG(rbf0i->crad); |
| 452 |
+ |
FVECT v0; |
| 453 |
+ |
float mv; |
| 454 |
+ |
ovec_from_pos(v0, rbf0i->gx, rbf0i->gy); |
| 455 |
+ |
for (j = 0; j < mtx_ncols(mig); j++) |
| 456 |
+ |
if ((mv = mtx_coef(mig,i,j)) > cthresh) { |
| 457 |
+ |
const RBFVAL *rbf1j = &mig->rbfv[1]->rbfa[j]; |
| 458 |
+ |
double rad2; |
| 459 |
+ |
FVECT v; |
| 460 |
+ |
int pos[2]; |
| 461 |
+ |
rad2 = R2ANG(rbf1j->crad); |
| 462 |
+ |
rad2 = rad0*rad0*(1.-t) + rad2*rad2*t; |
| 463 |
+ |
rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal * |
| 464 |
+ |
rad0*rad0/rad2; |
| 465 |
+ |
rbf->rbfa[n].crad = ANG2R(sqrt(rad2)); |
| 466 |
+ |
ovec_from_pos(v, rbf1j->gx, rbf1j->gy); |
| 467 |
+ |
geodesic(v, v0, v, t, GEOD_REL); |
| 468 |
+ |
pos_from_vec(pos, v); |
| 469 |
+ |
rbf->rbfa[n].gx = pos[0]; |
| 470 |
+ |
rbf->rbfa[n].gy = pos[1]; |
| 471 |
+ |
++n; |
| 472 |
+ |
} |
| 473 |
+ |
} |
| 474 |
+ |
rbf->vtotal *= mig->rbfv[0]->vtotal; /* turn ratio into actual */ |
| 475 |
+ |
return(rbf); |
| 476 |
+ |
memerr: |
| 477 |
+ |
fprintf(stderr, "%s: Out of memory in e_advect_rbf()\n", progname); |
| 478 |
+ |
exit(1); |
| 479 |
+ |
return(NULL); /* pro forma return */ |
| 480 |
+ |
} |
| 481 |
+ |
|
| 482 |
|
/* Clear our BSDF representation and free memory */ |
| 483 |
|
void |
| 484 |
|
clear_bsdf_rep(void) |
