src/gmxlib/selection/sm_insolidangle.c

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