diff options
Diffstat (limited to 'src/GeneralizedPolynomial-Uniform/molecule_posn.c')
| -rw-r--r-- | src/GeneralizedPolynomial-Uniform/molecule_posn.c | 177 |
1 files changed, 94 insertions, 83 deletions
diff --git a/src/GeneralizedPolynomial-Uniform/molecule_posn.c b/src/GeneralizedPolynomial-Uniform/molecule_posn.c index e007396..6d42ef3 100644 --- a/src/GeneralizedPolynomial-Uniform/molecule_posn.c +++ b/src/GeneralizedPolynomial-Uniform/molecule_posn.c @@ -8,7 +8,6 @@ @@*/ #include <math.h> -#include <limits.h> #include <stdlib.h> #include <stdio.h> @@ -73,6 +72,18 @@ static const char *rcsid = "$Header$"; center, i.e. the above diagram has columns labelled with [i+m]. @enddesc + @hdate 27.Jan.2003 + @hauthor Jonathan Thornburg <jthorn@aei.mpg.de> + @hdesc Complete rewrite: now supports + @var boundary_off_centering_tolerance_min, + @var boundary_off_centering_tolerance_max, + @var boundary_extrapolation_tolerance_min, + @var boundary_extrapolation_tolerance_max, + also change to returning status code, + also drop returning @var min_m and @var max_m + since they were never used. + @endhdesc + @var grid_origin @vdesc The floating-point coordinate x of the grid point i=0. @vtype fp grid_origin @@ -98,50 +109,53 @@ static const char *rcsid = "$Header$"; @vtype int molecule_size @endvar - @var out_of_range_tolerance_min, out_of_range_tolerance_max - @vdesc Specifies how out-of-range interpolation points should - be handled for the {minimum,maximum} ends of the grid - respectively. - If out_of_range_tolerance >= 0.0, - then an interpolation point is considered to be - "out of range" if and only if the interpolation - point is > out_of_range_tolerance * grid_delta - outside the grid in any coordinate. - If out_of_range_tolerance == -1.0, - then an interpolation point is considered to be - "out of range" if and only if a centered molecule - would require data from outside the grid. - Other values of out_of_range_tolerance are illegal. - @vtype fp out_of_range_tolerance_min, out_of_range_tolerance_max + @var boundary_off_centering_tolerance_min, + boundary_off_centering_tolerance_max + @vdesc Specifies how many grid spacings the interpolation point + may be beyond the default-centering region before being + declared "out of range" on the {minimum,maximum} boundaries + of the grid respectively. + @vtype fp boundary_off_centering_tolerance_min, + boundary_off_centering_tolerance_max + @endvar + + @var boundary_extrapolation_tolerance_min, + boundary_extrapolation_tolerance_max + @vdesc Specifies how many grid spacings the interpolation point + may be beyond the grid boundary before being declared + "out of range" on the {minimum,maximum} boundaries + of the grid respectively. + @vtype fp boundary_extrapolation_tolerance_min, + boundary_extrapolation_tolerance_max @endvar @var x - @vdesc The floating-point coordinate x at which we would like to center - the molecule. + @vdesc The floating-point coordinate of the interpolation point. @vtype fp x @endvar - @var x_rel - @vdesc A pointer to a floating-point value where this function should - store the x coordinate relative to the molecule center, measured - in units of the grid spacing; or NULL to skip storing this. - @vtype fp *x_rel + @var i_center + @vdesc A pointer to an value where this function should + store the integer coordinate of the molecule center, + or NULL to skip storing this + @vtype int *i_center @vio pointer to out @endvar - @var molecule_m_min, molecule_m_max - @vdesc A pointer to an int where this function should store the - {minimum,maximum} molecule coordinate m of the molecule; + @var x_rel + @vdesc A pointer to where this function should store the + interpolation point's x coordinate relative to the + molecule center, measured in units of the grid spacing; or NULL to skip storing this. - @vtype int *molecule_m_min, *molecule_m_max + @vtype fp *x_rel @vio pointer to out @endvar @returntype int @returndesc - This function returns the integer coordinate i_center at which - the molecule should be centered, or one of the error codes defined - in "InterpLocalUniform.h": + This function returns 0 if the interpolation point is "in range", + or one of the (negative) error codes defined in "InterpLocalUniform.h" + if the interpolation point is "out of range": MOLECULE_POSN_ERROR_X_LT_MIN if x is out-of-range on the min end of the grid (i.e. x < the minimum allowable x) @@ -155,73 +169,70 @@ static const char *rcsid = "$Header$"; int LocalInterp_molecule_posn(fp grid_origin, fp grid_delta, int grid_i_min, int grid_i_max, int molecule_size, - fp out_of_range_tolerance_min, - fp out_of_range_tolerance_max, + fp boundary_off_centering_tolerance_min, + fp boundary_off_centering_tolerance_max, + fp boundary_extrapolation_tolerance_min, + fp boundary_extrapolation_tolerance_max, fp x, - fp *x_rel, - int *molecule_m_min, int *molecule_m_max) + int* i_center, fp* x_rel) { -/* min/max molecule coordinates */ -const int m_max = (molecule_size >> 1); /* a positive number */ -const int m_min = m_max - molecule_size + 1; /* a negative number */ +/* molecule radia (inherently positive) in +/- directions */ +const int mr_plus = (molecule_size >> 1); +const int mr_minus = molecule_size - mr_plus - 1; -/* integer coordinate i of point x, as a floating-point number */ -const fp fp_i = (x - grid_origin) / grid_delta; +/* range of x_rel for which this molecule is centered, cf. diagram in header comment */ +const fp centered_min_x_rel = IS_EVEN(molecule_size) ? 0.0 : -0.5; +const fp centered_max_x_rel = IS_EVEN(molecule_size) ? 1.0 : +0.5; -/* is point x out-of-range? */ -if (out_of_range_tolerance_min >= 0.0) - then { - const fp fp_effective_grid_i_min - = ((fp) grid_i_min) - out_of_range_tolerance_min; - if (fp_i < fp_effective_grid_i_min) - then return INT_MIN; /*** ERROR RETURN ***/ - } -if (out_of_range_tolerance_max >= 0.0) - then { - const fp fp_effective_grid_i_max - = ((fp) grid_i_max) + out_of_range_tolerance_max; - if (fp_i > fp_effective_grid_i_max) - then return INT_MAX; /*** ERROR RETURN ***/ - } +/* range of i where we *could* center the molecule, as floating-point numbers */ +const fp fp_centered_min_possible_i = grid_i_min + mr_minus + centered_min_x_rel; +const fp fp_centered_max_possible_i = grid_i_max - mr_plus + centered_max_x_rel; -/* integer coordinate i where we will center the molecule */ -/* (see diagram in header comment for explanation) */ -/* ... as a floating-point number */ - { -const fp fp_i_center = IS_EVEN(molecule_size) ? floor(fp_i) : JT_ROUND(fp_i); -/* ... as an integer */ -int i_center = (int) fp_i_center; +/* integer coordinate i of interpolation point, as a floating-point number */ +const fp fp_i = (x - grid_origin) / grid_delta; -/* min/max integer coordinates in molecule assuming it's fully centered */ -const int centered_i_min = i_center + m_min; -const int centered_i_max = i_center + m_max; +/* is the molecule larger than the grid? */ +if (molecule_size > HOW_MANY_IN_RANGE(grid_i_min,grid_i_max)) + then return MOLECULE_POSN_ERROR_GRID_TINY; /*** ERROR RETURN ***/ -/* check if off-centered molecules are forbidden? */ -if ((out_of_range_tolerance_min == -1.0) && (centered_i_min < grid_i_min)) +/* is interpolation point x beyond the extrapolation tolerance? */ +if (fp_i < (fp)grid_i_min - boundary_extrapolation_tolerance_min) then return MOLECULE_POSN_ERROR_X_LT_MIN; /*** ERROR RETURN ***/ -if ((out_of_range_tolerance_max == -1.0) && (centered_i_max > grid_i_max)) +if (fp_i > (fp)grid_i_max + boundary_extrapolation_tolerance_max) then return MOLECULE_POSN_ERROR_X_GT_MAX; /*** ERROR RETURN ***/ -/* off-center as needed if we're close to the edge of the grid */ +/* is interpolation point x beyond the off-centering tolerance? */ +if (fp_i < fp_centered_min_possible_i - boundary_off_centering_tolerance_min) + then return MOLECULE_POSN_ERROR_X_LT_MIN; /*** ERROR RETURN ***/ +if (fp_i > fp_centered_max_possible_i + boundary_off_centering_tolerance_max) + then return MOLECULE_POSN_ERROR_X_GT_MAX; /*** ERROR RETURN ***/ + +/* now choose the actual molecule position/centering: */ +/* first set up a centered molecule */ { -int shift = 0; -if (molecule_size > HOW_MANY_IN_RANGE(grid_i_min,grid_i_max)) - then return MOLECULE_POSN_ERROR_GRID_TINY; /*** ERROR RETURN ***/ -if (centered_i_min < grid_i_min) - then shift = grid_i_min - centered_i_min; /* a positive number */ -if (centered_i_max > grid_i_max) - then shift = grid_i_max - centered_i_max; /* a negative number */ -i_center += shift; +fp fp_i_center = IS_EVEN(molecule_size) /* ... as a floating-point number */ + ? floor(fp_i) + : JT_ROUND(fp_i); +int int_i_center = (int) fp_i_center; /* ... as an integer */ + +/* then clamp the molecule at the grid boundaries */ +if (int_i_center - mr_minus < grid_i_min) + then { + int_i_center = grid_i_min + mr_minus; + fp_i_center = (fp) int_i_center; + } +if (int_i_center + mr_plus > grid_i_max) + then { + int_i_center = grid_i_max - mr_plus; + fp_i_center = (fp) int_i_center; + } /* store the results */ +if (i_center != NULL) + then *i_center = int_i_center; if (x_rel != NULL) - then *x_rel = fp_i - ((fp) i_center); -if (molecule_m_min != NULL) - then *molecule_m_min = m_min; -if (molecule_m_max != NULL) - then *molecule_m_max = m_max; + then *x_rel = fp_i - fp_i_center; -return i_center; - } +return 0; /*** NORMAL RETURN ***/ } } |