source/specialunits.cc

   1 /* This file is part of Shapes.
   2  *
   3  * Shapes is free software: you can redistribute it and/or modify
   4  * it under the terms of the GNU General Public License as published by
   5  * the Free Software Foundation, either version 3 of the License, or
   6  * any later version.
   7  *
   8  * Shapes is distributed in the hope that it will be useful,
   9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11  * GNU General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU General Public License
  14  * along with Shapes.  If not, see <http://www.gnu.org/licenses/>.
  15  *
  16  * Copyright 2008 Henrik Tidefelt
  17  */
  18
  19 #include <cmath>
  20
  21 #include "specialunits.h"
  22 #include "consts.h"
  23 #include "isnan.h"
  24
  25 #include <iostream> /* I get size_t from here... */
  26 #include <limits>
  27
  28 using namespace Shapes;
  29
  30 double
  31 Computation::specialUnitCircleHandle( double a0 )
  32 {
  33         /* a0 is half the angle of the circular arc segment.
  34          */
  35         double t = a0 / Computation::RREL_TH_STEP;
  36         size_t i = static_cast< size_t >( t );
  37         t -= static_cast< double >( i );
  38         return ( 1 - t ) * Computation::RREL_TABLE[ i ] + t * Computation::RREL_TABLE[ i + 1 ];
  39 }
  40
  41 double
  42 Computation::specialUnitCorrection( double a0, double a1 )
  43 {
  44         const double k2 = -2.0729490168751576;
  45         const double k3 = -1.381966011250105;
  46         const double wInv = 0.7821424388422419;
  47
  48         if( a0 == a1 )
  49                 {
  50                         return 1;
  51                 }
  52
  53         double aRel = wInv * ( a1 - a0 ) / a0 ;
  54         double x = 2 * exp( aRel ) / ( 1 + exp( aRel ) ) - 1;
  55         if( IS_NAN( x ) )
  56                 {
  57                         if( a1 > 0 )
  58                                 {
  59                                         x = 1;
  60                                 }
  61                         else
  62                                 {
  63                                         x = -1;
  64                                 }
  65                 }
  66
  67         double xSquare = x * x;
  68         if( aRel < 0 )
  69                 {
  70                         return 1 + k2 * xSquare + k3 * x * xSquare;
  71                 }
  72         return 1 + 3 * xSquare - 2 * x * xSquare;
  73 }
  74
  75 double
  76 Computation::specialUnitNoInflexion( double a0, double a1 )
  77 {
  78         {
  79                 /* tmp is the opposite handles angle.   It will be in the range [ -\pi, 2\pi ], so the
  80                  * handle rays will intersect iff tmp \in [ 0, \pi - a0 ]
  81                  */
  82                 double tmp = a0 + a1;
  83                 if( tmp < 0 ||
  84                                 tmp > M_PI - a0 )
  85                         {
  86                                 return HUGE_VAL;
  87                         }
  88         }
  89         return sin( a0 + a1 ) / sin( a0 + a0 + a1 );
  90 }
  91