source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/detector/G4Detector/G4fresnellens/src/LensUtils.cc @ 117

#include "LensUtils.h"
#include <globals.hh>
#include <geomdefs.hh>

using namespace G4FresnelLens;
using std::swap;

#define NUM_TOLERANCE 1.e-15

//__________________________________________________________________________________
int G4FresnelLens::SolveQuadratic(double a, double b, double c, double& x1, double& x2)
{
    //
    // this function tries to avoid "catastrophic cancellation" due to subtraction
    // of almost equal numbers
    //
    if(a==0.0) {
        if(b==0.0) return 0;
        x1=x2=-c/b;
        return 1;
    }


    double b2a=b/a*0.5;
    double ca=c/a;
    double d=b2a*b2a-ca;
    if( d<0.0 ) return 0;
    if( d==0 ) {
        x1=x2=-b2a;
        return 1;
    }
    d=sqrt(d);
    if ( b2a<0 ) {
        x2=-b2a+d;
        x1=ca/x2;
    }
    else {
        x1=-b2a-d;
        x2=ca/x1;
    }
    return 2;
}

//__________________________________________________________________________________
int G4FresnelLens::BinarySearch(double* a, double x, int aleft, int aright){
//
// adopted from Ltrace Lrt_bsearch
// int Lrt_bsearch(double a[], double x, int left, int right)
  int center,left(aleft),right(aright);

  if(x<a[left])  return -1;
  if(x>a[right]) return -2;

  while(right-left>1){
      center=(left+right)/2;
      if(x<a[center]) right=center;
      else left=center;
  }

  return left;
}