source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/detector/optics/src/Nutils.cc @ 117

/* -------------------------------------------------------------------------
 *   utils.c
 *
 *   --- miscellaneous functions for raytracing
 *
 * Copyright (c) 2000-2005 N.Sakaki, Y.Takizawa, Y.Kawasaki
 * All rights reserved.
 * $Id: utils.c,v 1.3 2005/04/04 00:36:27 sakaki Exp sakaki $
 * -------------------------------------------------------------------------
 */
#include <math.h>
#include "Nutils.hh"
using namespace NTraceLens;
/*
 * ------------------------------------------------------------------------
 *  binary search
 *    find an index "i" that satisfy a[i] <= x < a[i+1]
 * ------------------------------------------------------------------------
 */
int NTraceLens::rt_bsearch(double a[], double x, int left, int right)
{
  int center;

  if(x<a[left])
    return -1;
  if(x>a[right])
    return right+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 NTraceLens::SolveQuadratic(double aa[3], double x[2])
{
  double 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=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 NTraceLens::vdot(double a[3], double b[3])
{
  return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
}