src/snark/cin.cpp

   1 /*
   2  ***********************************************************
   3  $SNARK_Header: S N A R K  1 4 - A PICTURE RECONSTRUCTION PROGRAM $
   4  $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/cin.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  THIS FUNCTION CALCULATE THE INTEGRAL OF (1.-COS(U))/U W.R.T. U
  11  FROM 0 TO X USING RATIONAL APPROXIMATIONS
  12  REFERENCES :
  13  1) POLYNOMIAL APPROXIMATIONS TO INTEGRAL TRANSFORMS,
  14  J. WIMP
  15  MATHEMATICS OF COMPUTATION, 15, 1961, P.P.177
  16  2) APPROXIMATIONS FOR DIGITAL COMPUTERS,
  17  C. HASTINGS
  18  PRINCETON UNIVERSITY PRESS, 1955, P.P.198
  19  */
  20
  21 #include <cstdio>
  22 #include <cmath>
  23
  24 #include "blkdta.h"
  25
  26 #include "cin.h"
  27
  28 REAL cin(REAL x)
  29 {
  30
  31         static const REAL gamma = (REAL) 0.5772156649;
  32         static const REAL a2n0 = (REAL) 0.1352962627;
  33         REAL ang;
  34         REAL anginc;
  35         static const REAL a2n[6] = { (REAL) -0.4232751922, (REAL) 0.0182227219, (REAL) 0.000415765, (REAL) 56716.e-10, (REAL) -511.e-10, (REAL) 3.e-10 };
  36
  37         REAL xsq;
  38         REAL f;
  39         REAL g;
  40
  41         if (x == 0)
  42                 return 0;
  43
  44         x = (REAL) fabs(x);
  45
  46         if (x < 1)
  47         {
  48
  49                 REAL cin_tmp = gamma - a2n0;
  50                 ang = 0.;
  51                 anginc = (REAL) 2.0 * (REAL) atan(sqrt(4. - x * x) / x);
  52                 for (int i = 0; i < 6; i++)
  53                 {
  54                         ang += anginc;
  55                         cin_tmp -= a2n[i] * (REAL) cos(ang);
  56                 }
  57                 return cin_tmp;
  58         }
  59
  60         xsq = x * x;
  61
  62         f = ((xsq + (REAL) 7.241163) * xsq + (REAL) 2.463936) / ((xsq + (REAL) 9.068580) * xsq + (REAL) 7.157433) / x;
  63
  64         g = ((xsq + (REAL) 7.547478) * xsq + (REAL) 1.564072) / ((xsq + (REAL) 12.723684) * xsq + (REAL) 15.723606) / xsq;
  65
  66         return gamma + (REAL) log(x) - f * (REAL) sin(x) + g * (REAL) cos(x);
  67 }