Line data    Source code

       1             : /*
       2             :  *  This file is part of libcxxsupport.
       3             :  *
       4             :  *  libcxxsupport is free software; you can redistribute it and/or modify
       5             :  *  it under the terms of the GNU General Public License as published by
       6             :  *  the Free Software Foundation; either version 2 of the License, or
       7             :  *  (at your option) any later version.
       8             :  *
       9             :  *  libcxxsupport is distributed in the hope that it will be useful,
      10             :  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             :  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      12             :  *  GNU General Public License for more details.
      13             :  *
      14             :  *  You should have received a copy of the GNU General Public License
      15             :  *  along with libcxxsupport; if not, write to the Free Software
      16             :  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
      17             :  */
      18             : 
      19             : /*
      20             :  *  libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
      21             :  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
      22             :  *  (DLR).
      23             :  */
      24             : 
      25             : /*
      26             :  *  Copyright (C) 2011 Max-Planck-Society
      27             :  *  \author Martin Reinecke
      28             :  */
      29             : 
      30             : #include "healpix_base/geom_utils.h"
      31             : #include "healpix_base/arr.h"
      32             : 
      33             : 
      34             : namespace healpix{
      35             : using namespace std;
      36             : 
      37             : 
      38           0 : void get_circle (const arr<vec3> &point, tsize q1, tsize q2, vec3 &center,
      39             :   double &cosrad)
      40             :   {
      41           0 :   center = (point[q1]+point[q2]).Norm();
      42           0 :   cosrad = dotprod(point[q1],center);
      43           0 :   for (tsize i=0; i<q1; ++i)
      44           0 :     if (dotprod(point[i],center)<cosrad) // point outside the current circle
      45             :       {
      46           0 :       center=crossprod(point[q1]-point[i],point[q2]-point[i]).Norm();
      47           0 :       cosrad=dotprod(point[i],center);
      48           0 :       if (cosrad<0)
      49           0 :         { center.Flip(); cosrad=-cosrad; }
      50             :       }
      51           0 :   }
      52           0 : void get_circle (const arr<vec3> &point, tsize q, vec3 &center,
      53             :   double &cosrad)
      54             :   {
      55           0 :   center = (point[0]+point[q]).Norm();
      56           0 :   cosrad = dotprod(point[0],center);
      57           0 :   for (tsize i=1; i<q; ++i)
      58           0 :     if (dotprod(point[i],center)<cosrad) // point outside the current circle
      59           0 :       get_circle(point,i,q,center,cosrad);
      60           0 :   }
      61             : 
      62             : 
      63           0 : void find_enclosing_circle (const arr<vec3> &point, vec3 &center,
      64             :   double &cosrad)
      65             :   {
      66             :   tsize np=point.size();
      67           0 :   planck_assert(np>=3,"too few points");
      68           0 :   center = (point[0]+point[1]).Norm();
      69           0 :   cosrad = dotprod(point[0],center);
      70           0 :   for (tsize i=2; i<np; ++i)
      71           0 :     if (dotprod(point[i],center)<cosrad) // point outside the current circle
      72           0 :       get_circle(point,i,center,cosrad);
      73           0 :   }
      74             : } // namespace healpix
      75             :