source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/optics/src/Kutils.cc @ 114

/* -------------------------------------------------------------------------
 *   Kutils.c
 *
 *   --- miscellaneous functions for raytracing
 *
 * $Id: Kutils.cc 119 2003-09-23 13:19:04Z pallas $
 * $Id: Kutils.c,v 1.0 2003/09/21 J. W. Watts
 * -------------------------------------------------------------------------
 */
#include "TMath.h"
#include "Kutils.hh"

/*
 * ------------------------------------------------------------------------
 *  binary search
 *    find an index "i" that satisfy a[i] <= x < a[i+1]
 * ------------------------------------------------------------------------
 */
Int_t Krt_bsearch(Double_t a[], Double_t x, Int_t left, Int_t right)
{
  Int_t center;

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

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

  return left;
}

/* ------------------------------------------------------------------------
 * This function solves the quadratic equation
 *   a[0]*x^2 + a[1]*x + a[2] = 0
 * and returns the number of real roots and the roots themselves
 * in x[0] and x[1].
 * ------------------------------------------------------------------------
 */
Int_t KSolveQuadratic(Double_t aa[3], Double_t x[2])
{
  Double_t d, t, a, b, c;

  a= aa[0];
  b=-aa[1];
  c= aa[2];

  if(a==0.0)
    {
      if(b==0.0)
        return 0;

      x[0]=c/b;
      x[1]=x[0];

      return 1;
    }

  /* Discriminant */
  d = b*b-4.0*a*c;
  if(d<0.0)
    return 0;

  d=TMath::Sqrt(d);

  t=2.0*a;

  if(b>0.0)
    {
      x[1]=(b+d)/t;
      x[0]=(c/a)/x[1];
    }
  else
    {
      x[0]=(b-d)/t;
      x[1]=(c/a)/x[0];
    }

  if(d==0.0)
    return 1;
  else
    return 2;
}

/*
 * ------------------------------------------------------------------------
 *  calculate inner product of two vectors
 * ------------------------------------------------------------------------
 */
Double_t Kvdot(Double_t a[3], Double_t b[3])
{
  return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
}