src/gromacs/selection/sm_insolidangle.cpp

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