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