#include<math.h>

#include "genlimit.h"
static void gameansigma(double *mean,double *sigma,
			int nchan,int nens,const int nobs[],const EB* ens);
static double arcfreq(double y);
#define freq(x) (0.5*erfc(-0.707106781186547524*(x)))



/*         Joel Heinrich  24 March 2005

   returns crude Gaussian approximation to upper limit.
   For use as a starting point.
*/

double galim(double beta,int nchan,int nens,const int nobs[],const EB* ens) {
  double mean=0,sigma=0;
  gameansigma(&mean,&sigma,nchan,nens,nobs,ens);
  return mean-sigma*arcfreq( (1-freq(-mean/sigma))*(1-beta) );
}



static void gameansigma(double *mean,double *sigma,
		 int nchan,int nens,const int nobs[],const EB* ens) {

  double sum=0,sum2=0,vsum=0;
  const EB* p = ens;
  int i,j;
  for(i=0;i<nens;++i) {
    double s=0,s2=0;
    for(j=0;j<nchan;++j) {
      const int n = nobs[j];
      const double eps = p->e;
      s += (n-p->b)*eps/(n+1);
      s2 += eps*eps/(n+1);
      ++p;
    }
    s /= s2;
    vsum += 1/s2;
    sum += s;
    sum2 += s*s;
  }
  
  *mean = sum/nens;
  *sigma = sqrt(vsum/nens + sum2/nens - (*mean)*(*mean));
  return;
}

#define rdfreq(x) (exp(0.5*(x)*(x))*2.50662827463100050242)

#define C0 2.515517
#define C1 0.802853
#define C2 0.010328
#define D0 1.0
#define D1 1.432788
#define D2 0.189269
#define D3 0.001308

static double arcfreq(double y) {
  const double yy = (y>0.5) ? 1-y : y, t = sqrt(-2*log(yy));
  double x = (C0+t*(C1+t*C2))/(D0+t*(D1+t*(D2+t*D3))) - t;
  x -= (freq(x) - yy)*rdfreq(x);
  x -= (freq(x) - yy)*rdfreq(x);
  return (y>0.5) ? -x : x;
}