src/gmxlib/selection/sm_insolidangle.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2009, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /*! \page sm_insolidangle Selection method: insolidangle
  39  *
  40  * This method selects a subset of particles that are located in a solid
  41  * angle defined by a center and a set of points.
  42  * The solid angle is constructed as a union of small cones whose axis
  43  * goes through the center and a point.
  44  * So there's such a cone for each position, and a
  45  * point is in the solid angle if it lies within any of these cones.
  46  * The width of the cones can be adjusted.
  47  *
  48  * \internal
  49  *
  50  * The method is implemented by partitioning the surface of the unit sphere
  51  * into bins using the polar coordinates \f$(\theta, \phi)\f$.
  52  * The partitioning is always uniform in the zenith angle \f$\theta\f$,
  53  * while the partitioning in the azimuthal angle \f$\phi\f$ varies.
  54  * For each reference point, the unit vector from the center to the point
  55  * is constructed, and it is stored in all the bins that overlap with the
  56  * cone defined by the point.
  57  * Bins that are completely covered by a single cone are marked as such.
  58  * Checking whether a point is in the solid angle is then straightforward
  59  * with this data structure: one finds the bin that corresponds to the point,
  60  * and checks whether the bin is completely covered. If it is not, one
  61  * additionally needs to check whether it is within the specified cutoff of
  62  * any of the stored points.
  63  *
  64  * The above construction gives quite a lot of flexibility for constructing
  65  * the bins without modifying the rest of the code.
  66  * The current (quite inefficient) implementation is discussed below, but
  67  * it should be optimized to get the most out of the code.
  68  *
  69  * The current way of constructing the bins constructs the boundaries
  70  * statically: the bin size in the zenith direction is set to approximately
  71  * half the angle cutoff, and the bins in the azimuthal direction have
  72  * sizes such that the shortest edge of the bin is approximately equal to
  73  * half the angle cutoff (for the regions close to the poles, a single bin
  74  * is used).
  75  * Each reference point is then added to the bins as follows:
  76  *  -# Find the zenith angle range that is spanned by the cone centered at the
  77  *     point (this is simple addition/subtraction).
  78  *  -# Calculate the maximal span of the cone in the azimuthal direction using
  79  *     the formula
  80  *     \f[\sin \Delta \phi_{max} = \frac{\sin \alpha}{\sin \theta}\f]
  81  *     (a sine formula in spherical coordinates),
  82  *     where \f$\alpha\f$ is the width of the cone and \f$\theta\f$ is the
  83  *     zenith angle of the cone center.
  84  *     Similarly, the zenith angle at which this extent is achieved is
  85  *     calculated using
  86  *     \f[\cos \theta_{max} = \frac{\cos \theta}{\cos \alpha}\f]
  87  *     (Pythagoras's theorem in spherical coordinates).
  88  *  -# For each zenith angle bin that is at least partially covered by the
  89  *     cone, calculate the span of the cone at the edges using
  90  *     \f[\sin^2 \frac{\Delta \phi}{2} = \frac{\sin^2 \frac{\alpha}{2} - \sin^2 \frac{\theta - \theta'}{2}}{\sin \theta \sin \theta'}\f]
  91  *     (distance in spherical geometry),
  92  *     where \f$\theta'\f$ is the zenith angle of the bin edge.
  93  *     Treat zenith angle bins that are completely covered by the cone (in the
  94  *     case that the cone is centered close to the pole) as a special case.
  95  *  -# Using the values calculated above, loop through the azimuthal bins that
  96  *     are partially or completely covered by the cone and update them.
  97  *
  98  * The total solid angle (for covered fraction calculations) is estimated by
  99  * taking the total area of completely covered bins plus
 100  * half the area of partially covered bins.
 101  * The second one is an approximation, but should give reasonable estimates
 102  * for the averages as well as in cases where the bin size is small.
 103  */
 104 /*! \internal \file
 105  * \brief Implementation of the \ref sm_insolidangle "insolidangle"
 106  *   selection method.
 107  *
 108  * \todo
 109  * The implementation could be optimized quite a bit.
 110  *
 111  * \todo Move the covered fraction stuff somewhere else and make it more
 112  * generic (along the lines it is handled in selection.h and trajana.h).
 113  */
 114 #ifdef HAVE_CONFIG_H
 115 #include <config.h>
 116 #endif
 117
 118 #include <math.h>
 119
 120 #include <macros.h>
 121 #include <maths.h>
 122 #include <pbc.h>
 123 #include <physics.h>
 124 #include <smalloc.h>
 125 #include <vec.h>
 126
 127 #include <indexutil.h>
 128 #include <position.h>
 129 #include <selection.h>
 130 #include <selmethod.h>
 131
 132 #include "selelem.h"
 133
 134 /*! \internal \brief
 135  * Internal data structure for the \p insolidangle selection method.
 136  *
 137  * \see \c t_partition
 138  */
 139 typedef struct
 140 {
 141     /** Left edge of the partition. */
 142     real                left;
 143     /** Bin index corresponding to this partition. */
 144     int                 bin;
 145 } t_partition_item;
 146
 147 /*! \internal \brief
 148  * Internal data structure for the \p insolidangle selection method.
 149  *
 150  * Describes the surface partitioning within one slice along the zenith angle.
 151  * The slice from azimuthal angle \p p[i].left to \p p[i+1].left belongs to
 152  * bin \p p[i].bin.
 153  */
 154 typedef struct
 155 {
 156     /** Number of partition items (\p p contains \p n+1 items). */
 157     int                 n;
 158     /** Array of partition edges and corresponding bins. */
 159     t_partition_item   *p;
 160 } t_partition;
 161
 162 /*! \internal \brief
 163  * Internal data structure for the \p insolidangle selection method.
 164  *
 165  * Contains the reference points that partially cover a certain region on the
 166  * surface of the unit sphere.
 167  * If \p n is -1, the whole region described by the bin is covered.
 168  */
 169 typedef struct
 170 {
 171     /** Number of points in the array \p x, -1 if whole bin covered. */
 172     int   n;
 173     /** Number of elements allocated for \p x. */
 174     int   n_alloc;
 175     /** Array of points that partially cover the bin. */
 176     rvec *x;
 177 } t_spheresurfacebin;
 178
 179 /*! \internal \brief
 180  * Data structure for the \p insolidangle selection method.
 181  *
 182  * All angle values are in the units of radians.
 183  */
 184 typedef struct
 185 {
 186     /** Center of the solid angle. */
 187     gmx_ana_pos_t       center;
 188     /** Positions that span the solid angle. */
 189     gmx_ana_pos_t       span;
 190     /** Cutoff angle. */
 191     real                angcut;
 192     /** Estimate of the covered fraction. */
 193     real                cfrac;
 194
 195     /** Cutoff for the cosine (equals cos(angcut)). */
 196     real                distccut;
 197     /** Bin size to be used as the target bin size when constructing the bins. */
 198     real                targetbinsize;
 199
 200     /** Number of bins in the \p tbin array. */
 201     int                 ntbins;
 202     /** Size of one bin in the zenith angle direction. */
 203     real                tbinsize;
 204     /** Array of zenith angle slices. */
 205     t_partition        *tbin;
 206     /** Number of elements allocated for the \p bin array. */
 207     int                 maxbins;
 208     /** Number of elements used in the \p bin array. */
 209     int                 nbins;
 210     /** Array of individual bins. */
 211     t_spheresurfacebin *bin;
 212 } t_methoddata_insolidangle;
 213
 214 /** Allocates data for the \p insolidangle selection method. */
 215 static void *
 216 init_data_insolidangle(int npar, gmx_ana_selparam_t *param);
 217 /** Initializes the \p insolidangle selection method. */
 218 static int
 219 init_insolidangle(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
 220 /** Frees the data allocated for the \p insolidangle selection method. */
 221 static void
 222 free_data_insolidangle(void *data);
 223 /** Initializes the evaluation of the \p insolidangle selection method for a frame. */
 224 static int
 225 init_frame_insolidangle(t_topology *top, t_trxframe *fr, t_pbc *pbc, void *data);
 226 /** Internal helper function for evaluate_insolidangle(). */
 227 static gmx_bool
 228 accept_insolidangle(rvec x, t_pbc *pbc, void *data);
 229 /** Evaluates the \p insolidangle selection method. */
 230 static int
 231 evaluate_insolidangle(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 232                       gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data);
 233
 234 /** Calculates the distance between unit vectors. */
 235 static real
 236 sph_distc(rvec x1, rvec x2);
 237 /** Does a binary search on a \p t_partition to find a bin for a value. */
 238 static int
 239 find_partition_bin(t_partition *p, real value);
 240 /** Finds a bin that corresponds to a location on the unit sphere surface. */
 241 static int
 242 find_surface_bin(t_methoddata_insolidangle *surf, rvec x);
 243 /** Clears/initializes the bins on the unit sphere surface. */
 244 static void
 245 clear_surface_points(t_methoddata_insolidangle *surf);
 246 /** Frees memory allocated for storing the reference points in the surface bins. */
 247 static void
 248 free_surface_points(t_methoddata_insolidangle *surf);
 249 /** Adds a reference point to a given bin. */
 250 static void
 251 add_surface_point(t_methoddata_insolidangle *surf, int tbin, int pbin, rvec x);
 252 /** Marks a bin as completely covered. */
 253 static void
 254 mark_surface_covered(t_methoddata_insolidangle *surf, int tbin, int pbin);
 255 /** Helper function for store_surface_point() to update a single zenith angle bin. */
 256 static void
 257 update_surface_bin(t_methoddata_insolidangle *surf, int tbin,
 258                    real phi, real pdelta1, real pdelta2, real pdeltamax,
 259                    rvec x);
 260 /** Adds a single reference point and updates the surface bins. */
 261 static void
 262 store_surface_point(t_methoddata_insolidangle *surf, rvec x);
 263 /** Optimizes the surface bins for faster searching. */
 264 static void
 265 optimize_surface_points(t_methoddata_insolidangle *surf);
 266 /** Estimates the area covered by the reference cones. */
 267 static real
 268 estimate_covered_fraction(t_methoddata_insolidangle *surf);
 269 /** Checks whether a point lies within a solid angle. */
 270 static gmx_bool
 271 is_surface_covered(t_methoddata_insolidangle *surf, rvec x);
 272
 273 /** Parameters for the \p insolidangle selection method. */
 274 static gmx_ana_selparam_t smparams_insolidangle[] = {
 275     {"center", {POS_VALUE,   1, {NULL}}, NULL, SPAR_DYNAMIC},
 276     {"span",   {POS_VALUE,  -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_VARNUM},
 277     {"cutoff", {REAL_VALUE,  1, {NULL}}, NULL, SPAR_OPTIONAL},
 278 };
 279
 280 /** Help text for the \p insolidangle selection method. */
 281 static const char *help_insolidangle[] = {
 282     "SELECTING ATOMS IN A SOLID ANGLE[PAR]",
 283
 284     "[TT]insolidangle center POS span POS_EXPR [cutoff REAL][tt][PAR]",
 285
 286     "This keyword selects atoms that are within [TT]REAL[tt] degrees",
 287     "(default=5) of any position in [TT]POS_EXPR[tt] as seen from [TT]POS[tt]",
 288     "a position expression that evaluates to a single position), i.e., atoms",
 289     "in the solid angle spanned by the positions in [TT]POS_EXPR[tt] and",
 290     "centered at [TT]POS[tt].[PAR]"
 291
 292     "Technically, the solid angle is constructed as a union of small cones",
 293     "whose tip is at [TT]POS[tt] and the axis goes through a point in",
 294     "[TT]POS_EXPR[tt]. There is such a cone for each position in",
 295     "[TT]POS_EXPR[tt], and point is in the solid angle if it lies within any",
 296     "of these cones. The cutoff determines the width of the cones.",
 297 };
 298
 299 /** \internal Selection method data for the \p insolidangle method. */
 300 gmx_ana_selmethod_t sm_insolidangle = {
 301     "insolidangle", GROUP_VALUE, SMETH_DYNAMIC,
 302     asize(smparams_insolidangle), smparams_insolidangle,
 303     &init_data_insolidangle,
 304     NULL,
 305     &init_insolidangle,
 306     NULL,
 307     &free_data_insolidangle,
 308     &init_frame_insolidangle,
 309     NULL,
 310     &evaluate_insolidangle,
 311     {"insolidangle center POS span POS_EXPR [cutoff REAL]",
 312      asize(help_insolidangle), help_insolidangle},
 313 };
 314
 315 /*!
 316  * \param[in]     npar  Not used (should be 3).
 317  * \param[in,out] param Method parameters (should point to
 318  *   \ref smparams_insolidangle).
 319  * \returns Pointer to the allocated data (\ref t_methoddata_insolidangle).
 320  *
 321  * Allocates memory for a \ref t_methoddata_insolidangle structure and
 322  * initializes the parameter as follows:
 323  *  - \p center defines the value for t_methoddata_insolidangle::center.
 324  *  - \p span   defines the value for t_methoddata_insolidangle::span.
 325  *  - \p cutoff defines the value for t_methoddata_insolidangle::angcut.
 326  */
 327 static void *
 328 init_data_insolidangle(int npar, gmx_ana_selparam_t *param)
 329 {
 330     t_methoddata_insolidangle *data;
 331
 332     snew(data, 1);
 333     data->angcut = 5.0;
 334     param[0].val.u.p = &data->center;
 335     param[1].val.u.p = &data->span;
 336     param[2].val.u.r = &data->angcut;
 337     return data;
 338 }
 339
 340 /*!
 341  * \param   top  Not used.
 342  * \param   npar Not used.
 343  * \param   param Not used.
 344  * \param   data Pointer to \ref t_methoddata_insolidangle to initialize.
 345  * \returns 0 on success, -1 on failure.
 346  *
 347  * Converts t_methoddata_insolidangle::angcut to radians and allocates
 348  * and allocates memory for the bins used during the evaluation.
 349  */
 350 static int
 351 init_insolidangle(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
 352 {
 353     t_methoddata_insolidangle *surf = (t_methoddata_insolidangle *)data;
 354     int                        i, c;
 355
 356     if (surf->angcut <= 0)
 357     {
 358         fprintf(stderr, "error: angle cutoff should be > 0");
 359         return -1;
 360     }
 361
 362     surf->angcut *= DEG2RAD;
 363
 364     surf->distccut = -cos(surf->angcut);
 365     surf->targetbinsize = surf->angcut / 2;
 366     surf->ntbins = (int) (M_PI / surf->targetbinsize);
 367     surf->tbinsize = (180.0 / surf->ntbins)*DEG2RAD;
 368
 369     snew(surf->tbin, (int)(M_PI/surf->tbinsize) + 1);
 370     surf->maxbins = 0;
 371     for (i = 0; i < surf->ntbins; ++i)
 372     {
 373         c = max(sin(surf->tbinsize*i), sin(surf->tbinsize*(i+1)))
 374               * M_2PI / surf->targetbinsize + 1;
 375         snew(surf->tbin[i].p, c+1);
 376         surf->maxbins += c;
 377     }
 378     surf->nbins = 0;
 379     snew(surf->bin, surf->maxbins);
 380
 381     return 0;
 382 }
 383
 384 /*!
 385  * \param data Data to free (should point to a \ref t_methoddata_insolidangle).
 386  *
 387  * Frees the memory allocated for \c t_methoddata_insolidangle::center and
 388  * \c t_methoddata_insolidangle::span, as well as the memory for the internal
 389  * bin structure.
 390  */
 391 static void
 392 free_data_insolidangle(void *data)
 393 {
 394     t_methoddata_insolidangle *d = (t_methoddata_insolidangle *)data;
 395     int                        i;
 396
 397     if (d->tbin)
 398     {
 399         for (i = 0; i < d->ntbins; ++i)
 400         {
 401             sfree(d->tbin[i].p);
 402         }
 403         sfree(d->tbin);
 404     }
 405     free_surface_points(d);
 406     sfree(d->bin);
 407 }
 408
 409 /*!
 410  * \param[in]  top  Not used.
 411  * \param[in]  fr   Current frame.
 412  * \param[in]  pbc  PBC structure.
 413  * \param      data Should point to a \ref t_methoddata_insolidangle.
 414  * \returns    0 on success, a non-zero error code on error.
 415  *
 416  * Creates a lookup structure that enables fast queries of whether a point
 417  * is within the solid angle or not.
 418  */
 419 static int
 420 init_frame_insolidangle(t_topology *top, t_trxframe *fr, t_pbc *pbc, void *data)
 421 {
 422     t_methoddata_insolidangle *d = (t_methoddata_insolidangle *)data;
 423     rvec                       dx;
 424     int                        i;
 425
 426     free_surface_points(d);
 427     clear_surface_points(d);
 428     for (i = 0; i < d->span.nr; ++i)
 429     {
 430         if (pbc)
 431         {
 432             pbc_dx(pbc, d->span.x[i], d->center.x[0], dx);
 433         }
 434         else
 435         {
 436             rvec_sub(d->span.x[i], d->center.x[0], dx);
 437         }
 438         unitv(dx, dx);
 439         store_surface_point(d, dx);
 440     }
 441     optimize_surface_points(d);
 442     d->cfrac = -1;
 443     return 0;
 444 }
 445
 446 /*!
 447  * \param[in] x    Test point.
 448  * \param[in] pbc  PBC data (if NULL, no PBC are used).
 449  * \param[in] data Pointer to a \c t_methoddata_insolidangle data structure.
 450  * \returns   TRUE if \p x is within the solid angle, FALSE otherwise.
 451  */
 452 static gmx_bool
 453 accept_insolidangle(rvec x, t_pbc *pbc, void *data)
 454 {
 455     t_methoddata_insolidangle *d = (t_methoddata_insolidangle *)data;
 456     rvec                       dx;
 457
 458     if (pbc)
 459     {
 460         pbc_dx(pbc, x, d->center.x[0], dx);
 461     }
 462     else
 463     {
 464         rvec_sub(x, d->center.x[0], dx);
 465     }
 466     unitv(dx, dx);
 467     return is_surface_covered(d, dx);
 468 }
 469
 470 /*!
 471  * See sel_updatefunc() for description of the parameters.
 472  * \p data should point to a \c t_methoddata_insolidangle.
 473  *
 474  * Calculates which atoms in \p g are within the solid angle spanned by
 475  * \c t_methoddata_insolidangle::span and centered at
 476  * \c t_methoddata_insolidangle::center, and stores the result in \p out->u.g.
 477  */
 478 static int
 479 evaluate_insolidangle(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 480                       gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data)
 481 {
 482     t_methoddata_insolidangle *d = (t_methoddata_insolidangle *)data;
 483     int                        b;
 484
 485     out->u.g->isize = 0;
 486     for (b = 0; b < pos->nr; ++b)
 487     {
 488         if (accept_insolidangle(pos->x[b], pbc, data))
 489         {
 490             gmx_ana_pos_append(NULL, out->u.g, pos, b, 0);
 491         }
 492     }
 493     return 0;
 494 }
 495
 496 /*!
 497  * \param[in] sel Selection element to query.
 498  * \returns   TRUE if the covered fraction can be estimated for \p sel with
 499  *   _gmx_selelem_estimate_coverfrac(), FALSE otherwise.
 500  */
 501 gmx_bool
 502 _gmx_selelem_can_estimate_cover(t_selelem *sel)
 503 {
 504     t_selelem   *child;
 505     gmx_bool         bFound;
 506     gmx_bool         bDynFound;
 507
 508     if (sel->type == SEL_BOOLEAN && sel->u.boolt == BOOL_OR)
 509     {
 510         return FALSE;
 511     }
 512     bFound    = FALSE;
 513     bDynFound = FALSE;
 514     child     = sel->child;
 515     while (child)
 516     {
 517         if (child->type == SEL_EXPRESSION)
 518         {
 519             if (child->u.expr.method->name == sm_insolidangle.name)
 520             {
 521                 if (bFound || bDynFound)
 522                 {
 523                     return FALSE;
 524                 }
 525                 bFound = TRUE;
 526             }
 527             else if (child->u.expr.method
 528                      && (child->u.expr.method->flags & SMETH_DYNAMIC))
 529             {
 530                 if (bFound)
 531                 {
 532                     return FALSE;
 533                 }
 534                 bDynFound = TRUE;
 535             }
 536         }
 537         else if (!_gmx_selelem_can_estimate_cover(child))
 538         {
 539             return FALSE;
 540         }
 541         child = child->next;
 542     }
 543     return TRUE;
 544 }
 545
 546 /*!
 547  * \param[in] sel Selection for which the fraction should be calculated.
 548  * \returns Fraction of angles covered by the selection (between zero and one).
 549  *
 550  * The return value is undefined if _gmx_selelem_can_estimate_cover() returns
 551  * FALSE.
 552  * Should be called after gmx_ana_evaluate_selections() has been called for the
 553  * frame.
 554  */
 555 real
 556 _gmx_selelem_estimate_coverfrac(t_selelem *sel)
 557 {
 558     t_selelem   *child;
 559     real         cfrac;
 560
 561     if (sel->type == SEL_EXPRESSION && sel->u.expr.method->name == sm_insolidangle.name)
 562     {
 563         t_methoddata_insolidangle *d = (t_methoddata_insolidangle *)sel->u.expr.mdata;
 564         if (d->cfrac < 0)
 565         {
 566             d->cfrac = estimate_covered_fraction(d);
 567         }
 568         return d->cfrac;
 569     }
 570     if (sel->type == SEL_BOOLEAN && sel->u.boolt == BOOL_NOT)
 571     {
 572         cfrac = _gmx_selelem_estimate_coverfrac(sel->child);
 573         if (cfrac < 1.0)
 574         {
 575             return 1 - cfrac;
 576         }
 577         return 1;
 578     }
 579
 580     /* Here, we assume that the selection is simple enough */
 581     child = sel->child;
 582     while (child)
 583     {
 584         cfrac = _gmx_selelem_estimate_coverfrac(child);
 585         if (cfrac < 1.0)
 586         {
 587             return cfrac;
 588         }
 589         child = child->next;
 590     }
 591     return 1.0;
 592 }
 593
 594 /*!
 595  * \param[in] x1  Unit vector 1.
 596  * \param[in] x2  Unit vector 2.
 597  * \returns   Minus the dot product of \p x1 and \p x2.
 598  *
 599  * This function is used internally to calculate the distance between the
 600  * unit vectors \p x1 and \p x2 to find out whether \p x2 is within the
 601  * cone centered at \p x1. Currently, the cosine of the angle is used
 602  * for efficiency, and the minus is there to make it behave like a normal
 603  * distance (larger values mean longer distances).
 604  */
 605 static real
 606 sph_distc(rvec x1, rvec x2)
 607 {
 608     return -iprod(x1, x2);
 609 }
 610
 611 /*!
 612  * \param[in] p     Partition to search.
 613  * \param[in] value Value to search for.
 614  * \returns   The partition index in \p p that contains \p value.
 615  *
 616  * If \p value is outside the range of \p p, the first/last index is returned.
 617  * Otherwise, the return value \c i satisfies \c p->p[i].left<=value and
 618  * \c p->p[i+1].left>value
 619  */
 620 static int
 621 find_partition_bin(t_partition *p, real value)
 622 {
 623     int pmin, pmax, pbin;
 624
 625     /* Binary search the partition */
 626     pmin = 0; pmax = p->n;
 627     while (pmax > pmin + 1)
 628     {
 629         pbin = pmin + (pmax - pmin) / 2;
 630         if (p->p[pbin].left <= value)
 631         {
 632             pmin = pbin;
 633         }
 634         else
 635         {
 636             pmax = pbin;
 637         }
 638     }
 639     pbin = pmin;
 640     return pbin;
 641 }
 642
 643 /*!
 644  * \param[in] surf  Surface data structure to search.
 645  * \param[in] x     Unit vector to find.
 646  * \returns   The bin index that contains \p x.
 647  *
 648  * The return value is an index to the \p surf->bin array.
 649  */
 650 static int
 651 find_surface_bin(t_methoddata_insolidangle *surf, rvec x)
 652 {
 653     real theta, phi;
 654     int  tbin, pbin;
 655
 656     theta = acos(x[ZZ]);
 657     phi = atan2(x[YY], x[XX]);
 658     tbin = floor(theta / surf->tbinsize);
 659     if (tbin >= surf->ntbins)
 660     {
 661         tbin = surf->ntbins - 1;
 662     }
 663     pbin = find_partition_bin(&surf->tbin[tbin], phi);
 664     return surf->tbin[tbin].p[pbin].bin;
 665 }
 666
 667 /*!
 668  * \param[in,out] surf Surface data structure.
 669  *
 670  * Clears the reference points from the bins and (re)initializes the edges
 671  * of the azimuthal bins.
 672  */
 673 static void
 674 clear_surface_points(t_methoddata_insolidangle *surf)
 675 {
 676     int i, j, c;
 677
 678     surf->nbins = 0;
 679     for (i = 0; i < surf->ntbins; ++i)
 680     {
 681         c = min(sin(surf->tbinsize*i), sin(surf->tbinsize*(i+1)))
 682               * M_2PI / surf->targetbinsize + 1;
 683         if (c <= 0)
 684         {
 685             c = 1;
 686         }
 687         surf->tbin[i].n = c;
 688         for (j = 0; j < c; ++j)
 689         {
 690             surf->tbin[i].p[j].left = -M_PI + j*M_2PI/c - 0.0001;
 691             surf->tbin[i].p[j].bin = surf->nbins;
 692             surf->bin[surf->nbins].n = 0;
 693             surf->nbins++;
 694         }
 695         surf->tbin[i].p[c].left = M_PI + 0.0001;
 696         surf->tbin[i].p[c].bin = -1;
 697     }
 698 }
 699
 700 /*!
 701  * \param[in,out] surf Surface data structure.
 702  */
 703 static void
 704 free_surface_points(t_methoddata_insolidangle *surf)
 705 {
 706     int i;
 707
 708     for (i = 0; i < surf->nbins; ++i)
 709     {
 710         if (surf->bin[i].x)
 711         {
 712             sfree(surf->bin[i].x);
 713         }
 714         surf->bin[i].n_alloc = 0;
 715         surf->bin[i].x = NULL;
 716     }
 717 }
 718
 719 /*!
 720  * \param[in,out] surf Surface data structure.
 721  * \param[in]     tbin Bin number in the zenith angle direction.
 722  * \param[in]     pbin Bin number in the azimuthal angle direction.
 723  * \param[in]     x    Point to store.
 724  */
 725 static void
 726 add_surface_point(t_methoddata_insolidangle *surf, int tbin, int pbin, rvec x)
 727 {
 728     int bin;
 729
 730     bin = surf->tbin[tbin].p[pbin].bin;
 731     /* Return if bin is already completely covered */
 732     if (surf->bin[bin].n == -1)
 733         return;
 734     /* Allocate more space if necessary */
 735     if (surf->bin[bin].n == surf->bin[bin].n_alloc) {
 736         surf->bin[bin].n_alloc += 10;
 737         srenew(surf->bin[bin].x, surf->bin[bin].n_alloc);
 738     }
 739     /* Add the point to the bin */
 740     copy_rvec(x, surf->bin[bin].x[surf->bin[bin].n]);
 741     ++surf->bin[bin].n;
 742 }
 743
 744 /*!
 745  * \param[in,out] surf Surface data structure.
 746  * \param[in]     tbin Bin number in the zenith angle direction.
 747  * \param[in]     pbin Bin number in the azimuthal angle direction.
 748  */
 749 static void
 750 mark_surface_covered(t_methoddata_insolidangle *surf, int tbin, int pbin)
 751 {
 752     int bin;
 753
 754     bin = surf->tbin[tbin].p[pbin].bin;
 755     surf->bin[bin].n = -1;
 756 }
 757
 758 /*!
 759  * \param[in,out] surf      Surface data structure.
 760  * \param[in]     tbin      Bin number in the zenith angle direction.
 761  * \param[in]     phi       Azimuthal angle of \p x.
 762  * \param[in]     pdelta1   Width of the cone at the lower edge of \p tbin.
 763  * \param[in]     pdelta2   Width of the cone at the uppper edge of \p tbin.
 764  * \param[in]     pdeltamax Max. width of the cone inside \p tbin.
 765  * \param[in]     x         Point to store (should have unit length).
 766  */
 767 static void
 768 update_surface_bin(t_methoddata_insolidangle *surf, int tbin,
 769                    real phi, real pdelta1, real pdelta2, real pdeltamax,
 770                    rvec x)
 771 {
 772     real pdelta, phi1, phi2;
 773     int  pbin1, pbin2, pbiniter, pbin;
 774
 775     /* Find the edges of the bins affected */
 776     pdelta = max(max(pdelta1, pdelta2), pdeltamax);
 777     phi1 = phi - pdelta;
 778     if (phi1 >= -M_PI)
 779     {
 780         pbin = find_partition_bin(&surf->tbin[tbin], phi1);
 781         pbin1 = pbin;
 782     }
 783     else
 784     {
 785         pbin = find_partition_bin(&surf->tbin[tbin], phi1 + M_2PI);
 786         pbin1 = pbin - surf->tbin[tbin].n;
 787     }
 788     phi2 = phi + pdelta;
 789     if (phi2 <= M_PI)
 790     {
 791         pbin2 = find_partition_bin(&surf->tbin[tbin], phi2);
 792     }
 793     else
 794     {
 795         pbin2 = find_partition_bin(&surf->tbin[tbin], phi2 - M_2PI);
 796         pbin2 += surf->tbin[tbin].n;
 797     }
 798     ++pbin2;
 799     if (pbin2 - pbin1 > surf->tbin[tbin].n)
 800     {
 801         pbin2 = pbin1 + surf->tbin[tbin].n;
 802     }
 803     /* Find the edges of completely covered region */
 804     pdelta = min(pdelta1, pdelta2);
 805     phi1 = phi - pdelta;
 806     if (phi1 < -M_PI)
 807     {
 808         phi1 += M_2PI;
 809     }
 810     phi2 = phi + pdelta;
 811     /* Loop over all affected bins */
 812     for (pbiniter = pbin1; pbiniter != pbin2; ++pbiniter, ++pbin)
 813     {
 814         /* Wrap bin around if end reached */
 815         if (pbin == surf->tbin[tbin].n)
 816         {
 817             pbin = 0;
 818             phi1 -= M_2PI;
 819             phi2 -= M_2PI;
 820         }
 821         /* Check if bin is completely covered and update */
 822         if (surf->tbin[tbin].p[pbin].left >= phi1
 823             && surf->tbin[tbin].p[pbin+1].left <= phi2)
 824         {
 825             mark_surface_covered(surf, tbin, pbin);
 826         }
 827         else
 828         {
 829             add_surface_point(surf, tbin, pbin, x);
 830         }
 831     }
 832 }
 833
 834 /*!
 835  * \param[in,out] surf Surface data structure.
 836  * \param[in]     x    Point to store (should have unit length).
 837  *
 838  * Finds all the bins covered by the cone centered at \p x and calls
 839  * update_surface_bin() to update them.
 840  */
 841 static void
 842 store_surface_point(t_methoddata_insolidangle *surf, rvec x)
 843 {
 844     real theta, phi;
 845     real pdeltamax, tmax;
 846     real theta1, theta2, pdelta1, pdelta2;
 847     int  tbin, pbin, bin;
 848
 849     theta = acos(x[ZZ]);
 850     phi = atan2(x[YY], x[XX]);
 851     /* Find the maximum extent in the phi direction */
 852     if (theta <= surf->angcut)
 853     {
 854         pdeltamax = M_PI;
 855         tmax = 0;
 856     }
 857     else if (theta >= M_PI - surf->angcut)
 858     {
 859         pdeltamax = M_PI;
 860         tmax = M_PI;
 861     }
 862     else
 863     {
 864         pdeltamax = asin(sin(surf->angcut) / sin(theta));
 865         tmax = acos(cos(theta) / cos(surf->angcut));
 866     }
 867     /* Find the first affected bin */
 868     tbin = max(floor((theta - surf->angcut) / surf->tbinsize), 0);
 869     theta1 = tbin * surf->tbinsize;
 870     if (theta1 < theta - surf->angcut)
 871     {
 872         pdelta1 = 0;
 873     }
 874     else
 875     {
 876         pdelta1 = M_PI;
 877     }
 878     /* Loop through all affected bins */
 879     while (tbin < ceil((theta + surf->angcut) / surf->tbinsize)
 880            && tbin < surf->ntbins)
 881     {
 882         /* Calculate the next boundaries */
 883         theta2 = (tbin+1) * surf->tbinsize;
 884         if (theta2 > theta + surf->angcut)
 885         {
 886             /* The circle is completely outside the cone */
 887             pdelta2 = 0;
 888         }
 889         else if (theta2 <= -(theta - surf->angcut)
 890                  || theta2 >= M_2PI - (theta + surf->angcut)
 891                  || tbin == surf->ntbins - 1)
 892         {
 893             /* The circle is completely inside the cone, or we are in the
 894              * 360 degree bin covering the pole. */
 895             pdelta2 = M_PI;
 896         }
 897         else
 898         {
 899             /* TODO: This formula is numerically unstable if theta is very
 900              * close to the pole.  In practice, it probably does not matter
 901              * much, but it would be nicer to adjust the theta bin boundaries
 902              * such that the case above catches this instead of falling through
 903              * here. */
 904             pdelta2 = 2*asin(sqrt(
 905                     (sqr(sin(surf->angcut/2)) - sqr(sin((theta2-theta)/2))) /
 906                     (sin(theta) * sin(theta2))));
 907         }
 908         /* Update the bin */
 909         if (tmax >= theta1 && tmax <= theta2)
 910         {
 911             update_surface_bin(surf, tbin, phi, pdelta1, pdelta2, pdeltamax, x);
 912         }
 913         else
 914         {
 915             update_surface_bin(surf, tbin, phi, pdelta1, pdelta2, 0, x);
 916         }
 917         /* Next bin */
 918         theta1 = theta2;
 919         pdelta1 = pdelta2;
 920         ++tbin;
 921     }
 922 }
 923
 924 /*!
 925  * \param[in,out] surf Surface data structure.
 926  *
 927  * Currently, this function does nothing.
 928  */
 929 static void
 930 optimize_surface_points(t_methoddata_insolidangle *surf)
 931 {
 932     /* TODO: Implement */
 933 }
 934
 935 /*!
 936  * \param[in] surf Surface data structure.
 937  * \returns   An estimate for the area covered by the reference points.
 938  */
 939 static real
 940 estimate_covered_fraction(t_methoddata_insolidangle *surf)
 941 {
 942     int  t, p, n;
 943     real cfrac, tfrac, pfrac;
 944
 945     cfrac = 0.0;
 946     for (t = 0; t < surf->ntbins; ++t)
 947     {
 948         tfrac = cos(t * surf->tbinsize) - cos((t+1) * surf->tbinsize);
 949         for (p = 0; p < surf->tbin[t].n; ++p)
 950         {
 951             pfrac = surf->tbin[t].p[p+1].left - surf->tbin[t].p[p].left;
 952             n = surf->bin[surf->tbin[t].p[p].bin].n;
 953             if (n == -1) /* Bin completely covered */
 954             {
 955                 cfrac += tfrac * pfrac;
 956             }
 957             else if (n > 0) /* Bin partially covered */
 958             {
 959                 cfrac += tfrac * pfrac / 2; /* A rough estimate */
 960             }
 961         }
 962     }
 963     return cfrac / (4*M_PI);
 964 }
 965
 966 /*!
 967  * \param[in] surf  Surface data structure to search.
 968  * \param[in] x     Unit vector to check.
 969  * \returns   TRUE if \p x is within the solid angle, FALSE otherwise.
 970  */
 971 static gmx_bool
 972 is_surface_covered(t_methoddata_insolidangle *surf, rvec x)
 973 {
 974     int  bin, i;
 975
 976     bin = find_surface_bin(surf, x);
 977     /* Check for completely covered bin */
 978     if (surf->bin[bin].n == -1)
 979     {
 980         return TRUE;
 981     }
 982     /* Check each point that partially covers the bin */
 983     for (i = 0; i < surf->bin[bin].n; ++i)
 984     {
 985         if (sph_distc(x, surf->bin[bin].x[i]) < surf->distccut)
 986         {
 987             return TRUE;
 988         }
 989     }
 990     return FALSE;
 991 }