| 51 |
|
if (i) spinvector(tvec, tvec, outvec, phinc); |
| 52 |
|
if (tvec[2] > 0 ^ output_orient > 0) |
| 53 |
|
continue; |
| 54 |
< |
sum += eval_rbfrep(rbf, tvec) * output_orient*tvec[2]; |
| 54 |
> |
sum += eval_rbfrep(rbf, tvec) * COSF(tvec[2]); |
| 55 |
|
++n; |
| 56 |
|
} |
| 57 |
|
if (n < 2) /* should never happen! */ |
| 64 |
|
est_DSFrad(const RBFNODE *rbf, const FVECT outvec) |
| 65 |
|
{ |
| 66 |
|
const double rad_epsilon = 0.03; |
| 67 |
< |
const double DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec) |
| 68 |
< |
* output_orient*outvec[2]; |
| 67 |
> |
const double DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec) * |
| 68 |
> |
COSF(outvec[2]); |
| 69 |
|
double inside_rad = rad_epsilon; |
| 70 |
|
double outside_rad = 0.5; |
| 71 |
|
double DSFinside = eval_DSFsurround(rbf, outvec, inside_rad); |
| 72 |
|
double DSFoutside = eval_DSFsurround(rbf, outvec, outside_rad); |
| 73 |
|
#define interp_rad inside_rad + (outside_rad-inside_rad) * \ |
| 74 |
|
(DSFtarget-DSFinside) / (DSFoutside-DSFinside) |
| 75 |
< |
/* interpolation search */ |
| 76 |
< |
while (outside_rad-inside_rad > rad_epsilon) { |
| 75 |
> |
/* Newton's method (sort of) */ |
| 76 |
> |
do { |
| 77 |
|
double test_rad = interp_rad; |
| 78 |
< |
double DSFtest = eval_DSFsurround(rbf, outvec, test_rad); |
| 79 |
< |
if (DSFtarget < DSFtest) { |
| 78 |
> |
double DSFtest; |
| 79 |
> |
if (test_rad >= outside_rad) |
| 80 |
> |
return(test_rad); |
| 81 |
> |
if (test_rad <= inside_rad) |
| 82 |
> |
return(test_rad*(test_rad>0)); |
| 83 |
> |
DSFtest = eval_DSFsurround(rbf, outvec, test_rad); |
| 84 |
> |
if (DSFtest > DSFtarget) { |
| 85 |
|
inside_rad = test_rad; |
| 86 |
|
DSFinside = DSFtest; |
| 87 |
|
} else { |
| 88 |
|
outside_rad = test_rad; |
| 89 |
|
DSFoutside = DSFtest; |
| 90 |
|
} |
| 91 |
< |
} |
| 91 |
> |
if (DSFoutside >= DSFinside) |
| 92 |
> |
return(test_rad); |
| 93 |
> |
} while (outside_rad-inside_rad > rad_epsilon); |
| 94 |
|
return(interp_rad); |
| 95 |
|
#undef interp_rad |
| 96 |
|
} |
| 107 |
|
|
| 108 |
|
if (dsf_list == NULL) { |
| 109 |
|
bsdf_spec_peak = 0; |
| 110 |
< |
bsdf_spec_crad = 0; |
| 110 |
> |
bsdf_spec_rad = 0; |
| 111 |
|
return; |
| 112 |
|
} |
| 113 |
|
for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) { |
| 119 |
|
++n; |
| 120 |
|
} |
| 121 |
|
bsdf_spec_peak = peak_sum/(double)n; |
| 122 |
< |
bsdf_spec_crad = ANG2R( rad_sum/(double)n ); |
| 122 |
> |
bsdf_spec_rad = rad_sum/(double)n; |
| 123 |
|
} |
| 124 |
|
|
| 125 |
|
/* Create a new migration holder (sharing memory for multiprocessing) */ |
