ctools

// from 68 and 80 % containment radii (in degrees), calculate the sigma and gamma

// parameters for a normalized king function (WITH THE SECANT METHOD, fails if r80 < r68*1.1)

int calcKingParametersViaSecant( double rad68, double rad80, double & sigma, double & gamma )

{

  // copied almost verbatim from

  // $CTOOLS/scripts/cta_root2caldb.py's root2psf_king() function

  // no point in continuing if only 0's

  if ( rad68 <= 0.0 || rad80 <= 0.0 )

  {

    printf("                                           calcKingParameters(%.2f %.2f %.2f %.2f): Warning, a containment radius is 0, setting g/s to 0s...\n", rad68, rad80, sigma, gamma ) ;

    gamma = 0.0 ;

    sigma = 0.0 ;

    return 1 ;

  }

  // initial ratios

  double a = 1.0 - 0.68 ;

  double b = 1.0 - 0.80 ;

  double c = (rad68*rad68) / (rad80*rad80) ;

  // solve equation (a^x-1)/(b^x-1)=c for x using secant

  // method.  Stop when we are better than 1e-6;

  double x,f ;

  double x1 = -0.5 ;

  double x2 = -10000.0 ; // if too close to x1=-0.5, fit will fail on some combinations of rad68 and rad80 (if rad68 and rad80 are within 5% of each other)

  double f1 = ( ( pow(a,x1)-1.0 ) / ( pow(b,x1) -1.0 ) ) - c ;

  double f2 = ( ( pow(a,x2)-1.0 ) / ( pow(b,x2) -1.0 ) ) - c ;

  int iter  = 0 ;

  while ( true )

  {

    x = x1 - f1 * (x1-x2)/(f1-f2) ;

    f = ( pow(a,x)-1.0 ) / ( pow(b,x)-1.0 ) - c ;

    iter += 1 ;

    if ( iter == 1 ) printf("\n");

    char xwarn[100]="" ;

    if ( x > 0.0 ) sprintf(xwarn,"%sx:%.8f%s", KRED, x, KNRM ) ;

    else           sprintf(xwarn,"x:%.8f"    , x             ) ;

    printf("              calcKingParameters(%7.3f %7.3f %7.3f %7.3f): iter:%d  diff:%.8f  %s  f1=%f  f2=%f...\n", rad68, rad80, sigma, gamma, iter, fabs(f), xwarn, f1, f2 ) ;

    if ( fabs(f) < 1.0e-6 )

    {

      //printf("              calcKingParameters(%7.3f %7.3f %7.3f %7.3f): Warning, breaking!...\n", rad68, rad80, sigma, gamma ) ;

      break ;

    }

    else

    {

      f2 = f1 ;

      x2 = x1 ;

      f1 = f  ;

      x1 = x  ;

    }

    // compute gamma

    if ( x < 0.0 )

    {

      gamma = 1.0 - (1.0/x) ;

    }

    else

    {

      gamma = 1.0 ;

      //printf("              calcKingParameters(%7.3f %7.3f %7.3f %7.3f): Warning, x(%7.4f) >= 0, setting gamma to 1...\n", rad68, rad80, sigma, gamma, x ) ;

    }

    // compute sigma

    double denom = 2.0 * gamma * ( pow(b,x)-1.0 ) ;

    if ( denom > 0.0 )

    {

      sigma = rad68 * sqrt(1.0/denom) ;

    }

    else

    {

      //printf("              calcKingParameters(%7.3f %7.3f %7.3f %7.3f): Warning, first denom is 0, attempting 2nd method...\n", rad68, rad80, sigma, gamma ) ;

      denom = 2.0 * gamma * (pow(b,x)-1.0) ;

      if ( denom > 0.0 )

      {

        sigma = rad80 * sqrt(1.0/denom) ;

      }

      else

      {

        //printf("              calcKingParameters(%7.3f %7.3f %7.3f %7.3f): Warning, denom is still 0.0, setting g/s to 0s...\n", rad68, rad80, sigma, gamma ) ;

        gamma = 0.0 ;

        sigma = 0.0 ;

        return 0 ;

      }

    }

  }

  return 0 ;

}