home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / autopsp.zip / REGRESS / CSIMP.CPP < prev    next >
C/C++ Source or Header  |  1995-07-22  |  3KB  |  113 lines
  1. // CSimpson
  2. // Method definition for CSimpson integration class
  3.  
  4. #include "csimp.h"
  5.  
  6. CSimpson::CSimpson(void)
  7. {
  8.   n=20;                           // Safe values.
  9.   result=0;                       // Ensure no uninitialization problems
  10.   cur_result=0;
  11.   old_result=0;
  12.   result=0;
  13.   width=0;
  14.   tolerance=0.005;                // Resonable tolerance.
  15. }
  16.  
  17. // I couldn't compile with out this place holder.
  18. double CSimpson::IntMe(double x)
  19. {
  20.   return(x);
  21. }
  22.  
  23. // This method handles the infinities and can add on the
  24. //   "other half" if the user passes in a non-zero total.
  25. //   Total should be the total area under the curve if known.
  26. double CSimpson::Integrate(double low, double high, double total)
  27. {
  28.   result=0;
  29.  
  30.   if (low==INFINITY && high==INFINITY){
  31.     result=total;                     // Easy case first
  32.     return(result);
  33.   }
  34.   else if (low==INFINITY && high<=0){
  35.     CalcSimpsons(high,0);
  36.     result += (total/2) - cur_result;
  37.   }
  38.   else if (low==INFINITY && high>0){
  39.     CalcSimpsons(0,high);
  40.     result += (total/2) + cur_result;
  41.   }
  42.   else if (low<=0 && high==INFINITY){
  43.     CalcSimpsons(low,0);
  44.     result += (total/2) + cur_result;
  45.   }
  46.   else if (low>0 && high==INFINITY){
  47.     CalcSimpsons(0,low);
  48.     result += (total/2) - cur_result;
  49.   }
  50.   else {
  51.     CalcSimpsons(low,high);
  52.     result=cur_result;
  53.   }
  54.   return(result);
  55. }
  56.  
  57.  
  58. // Calculates the integral using simpsons rule.
  59. void CSimpson::CalcSimpsons(double low, double high)
  60. {
  61.   double sum=0;
  62.   double slider;                  // Slides across the function.
  63.  
  64.   if(low==high){                  // Handles a special case.
  65.     cur_result=0;                 // Space between x and x is 0.
  66.     return;
  67.   }
  68.  
  69.   old_result=INFINITY;            // Avoid accidental immediate kick out.
  70.   n=10;                           // Starting at 10 makes inital run have 20 blocks.
  71.   do{
  72.     sum=0;                        // Summation starts at 0 for each attempt.
  73.     n*=2;                         // double number of blocks.
  74.     old_result=cur_result;        // Update old_result before changing cur_result.
  75.     width=(high-low)/n;           // Calc the width.
  76.     slider=0;                     // Slide starts at 0.
  77.  
  78.     sum+=IntMe(low+slider);       // Following loop implements simpsons algorythm.
  79.     slider+=width;
  80.     sum+=4*IntMe(low+slider);
  81.     while(low+slider <= high-(2*width)){
  82.       slider+=width;
  83.       sum+=2*IntMe(low+slider);
  84.       slider+=width;
  85.       sum+=4*IntMe(low+slider);
  86.     }
  87.     sum+=IntMe(high);
  88.  
  89.     cur_result=(width/3)*sum;     // Simpsons fudge factor.
  90.                                   // Loop until the value is sufficiently
  91.                                   //  cool to be accepted into the class.
  92.   } while(fabs(cur_result-old_result) >= tolerance);
  93. }
  94.  
  95. // Calculates the Gamma function based on the books description.
  96. // Will only handle integers and (integers + .5)...
  97. // (Accepts the set of all reals divisible by .5?)
  98. double CSimpson::Gamma(double x)
  99. {
  100.   if (x>1){
  101.     return((x-1)*Gamma(x-1));
  102.   }
  103.   else if (x==.5){
  104.     return(pow(M_PI,.5));
  105.   }
  106.   return(1);
  107. }
  108.  
  109. // Sets the private member: tolerance.
  110. void CSimpson::SetTolerance(double Tol)
  111. {
  112.   tolerance=Tol;
  113. }