/**************************************/
/* Estimate the CDF of Non-cen Chi^2  */ 
/* Poisson(Lambda) sum of Expo(theta) */
/**************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>
#include <math.h>

#define MIN(x,y) (((x) < (y)) ? (x): (y))

double nonchi(double z, double theta, double lamb);
double nonchicdf(double, double, double);
double mass(double theta, double lamb,double lower, double upper, int Nstep);
double tailmass(double theta, double lamb,double lower, double upper, int Nstep);

/*************
double nonchi(double z, double lamb)
{
  double old, mynew, mysum;
  int k;
  old = -99;
  mynew = lamb * z/2 * exp(-lamb-z/2);
  mysum = exp(-lamb-z/2) + mynew;
  k = 2;
  while(mynew > old*0.005){
    old = mynew;
    mynew *= lamb * z/2 /k /(k-1);
    mysum += mynew;
    k++;
  }
  return(mysum/ z);
}
*************/

double nonchi(double z, double theta, double lamb)
{
  double old, mynew, mysum;
  int k;
  if(z==0) return(0);
  old = -99;
  mynew = lamb * z*theta * exp(-lamb-z*theta);
  mysum = exp(-lamb-z*theta) + mynew;
  k = 2;
  while(mynew > old*0.005){
    old = mynew;
    mynew *= lamb * z*theta /k /(k-1);
    mysum += mynew;
    k++;
  }
  return(mysum/ z);
}

/************
double nonchicdf(double ga, double lamb)
{
  double x, z, mysum, N=30.0;
  int i;
@@  x = ga/theta /2; taoyun pointed out this was wrong! 
@@    x = ga *theta*2  I let theta=1
  x = 2 * ga;
  mysum = 0;
  for (z=x/N; z<=x; z += x/N) mysum += nonchi(z, lamb);
  return(mysum*(x/N));
}
*************/

double nonchicdf(double ga, double theta, double lamb)
{
  double x, z, mysum, *intvl;
  int i;
  double mean, std, lower,  upper, small, big, left, middle, right, temp;
  /*  x = ga/theta /2; taoyun pointed out this was wrong! 
      x = ga *theta*2  I let theta=1 */
  mysum = 0;

  /*** now the intervals need to be adjusted according to where the
    mass is ***
  for (z=ga/N; z<=ga; z += ga/N) mysum += nonchi(z, theta, lamb);
  ***** ABANDONED ****/

  mean = lamb/theta;
  std = sqrt(2*lamb)/theta;

  lower = mean - 8/3 * std;
  upper = mean + 8/3 * std;
  small = std/3;

  if(upper < ga) {	
    /*    temp = nonchi(ga + 15*std, theta, lamb);*/
    return(1 - tailmass(theta, lamb, ga, ga + 15*std, 5)
	   /* - temp*std/(nonchi(ga+14*std, theta,
					   lamb)-temp)*/);
  }

  if(lower > 0){
    if (lower >= ga)
      return( tailmass(theta,lamb, 0, ga, 5 ) ); /* leftmass(theta,lamb,lo, hi, #steps)*/
    left = tailmass(theta,lamb,0, lower, 5 );
    /*    if(upper > ga) {upper = ga}  no need */
    middle = mass(theta,lamb,lower, ga, (int)( (ga-lower)/small ) );
    /*right = tailmass(theta,lamb,upper, ga, MIN(5, (int)((ga-upper)/small)));*/
    return(left + middle);
  }
  else{ /* where lower <0 */
    /*   if(upper > ga) no need 
      return(mass(theta,lamb,0, ga, (int)((ga-0)/small)));*/
    /*      if (upper < 0) return(tailmass(theta,lamb,0, ga,
	    MIN(5,(int)(ga/small)))); no need either */
      return(mass(theta,lamb,0,ga,(int)(upper/small) ));
  }
}

double mass(double theta, double lamb,double lower, double upper, int Nstep)
{
  double small, mysum=0;
  int i;

  if(upper < lower) {nrerror2("in mass(theta,lamb,), upper < lower!\n"); exit(1);}
  else if(upper == lower) return(0);
  else { /* lower < upper normal case */
    if(Nstep < 0) {nrerror2("in mass(theta,lamb,), Nstep < 0! \n "); exit(1);}
    else if(0==Nstep) return((upper-lower)/2 * (nonchi(upper, theta, lamb) + nonchi(lower, theta, lamb)));
    else {
      small = (upper - lower)/ Nstep;
      for (i = 1; i<Nstep; i++) {
	mysum += nonchi(lower + i*small, theta, lamb);}
      mysum += 0.5 * (nonchi(lower, theta, lamb)+nonchi(upper, theta, lamb));
      return(mysum * small);
    }
  }
}

double tailmass(double theta, double lamb,double lower, double upper, int Nstep)
{
  /* in tail area, i use exponential fit to calculate the integral
     because the actual tail is approx exponential, slightly slower
     the area = (f(t2)-f(t1))/lam, 
     where lam = (log f(t2)-log f(t1) )/(t2 - t1)
     */

  double small, mysum=0, oldf, newf, lam;
  int i, start;

  start = 1;
  if(upper < lower) {nrerror2("in mass(theta,lamb,), upper < lower!\n"); exit(1);}
  else if(upper == lower) return(0);
  else { /* lower < upper normal case */
    if(Nstep < 0) {nrerror2("in mass(theta,lamb,), Nstep < 0! \n "); exit(1);}
    else {
      if(0==Nstep) Nstep = 1;
      small = (upper - lower)/ Nstep;
      if(0==lower) {
	oldf = nonchi(lower+small, theta, lamb);
	start = 2;
	mysum = oldf/2 * small;
      }
      else  oldf = nonchi(lower, theta, lamb);
      for (i = start; i<Nstep+1; i++) {
	newf = nonchi(lower + i*small, theta, lamb);
	lam = (log(newf) - log(oldf))/small;
	if(0==lam) mysum += newf * small;
	else  mysum += (newf - oldf)/lam;
      }
      return(mysum);
    }
  }
}


#undef MIN