src/gromacs/selection/sm_insolidangle.cpp

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