/*
** cc cluster_prob.c -o cluster_prob -lm
*/
#include <stdio.h>
#include <math.h>

double calc_pa(double K, double sld, double area, double t_val);

main( int argc, char **argv)
{
   double K, fwhm, sld, area, t_val, p_val;
   double fourLogTwo, fwhm_true;

   if (argc < 6 )
   {      
      fprintf( stderr, "Usage: %s %s\n", argv[0],
         "clusterSize area t_val -[sld | fwhm] val" );

      exit(1);
   }
   K = atof(argv[1]);
   area = atof(argv[2]);
   t_val = atof(argv[3]);

   fourLogTwo = (double)4.0 * log((double)(2.0));

   if ( 'f' == argv[4][1] )
   {
      double denom;
      fwhm = atof(argv[5]);

      denom = 1.0 + ((fwhm * fwhm) / fourLogTwo);
      sld = (double)1.0 /denom;
      printf( "sld( fwhm_gauss = %lf ) = %lf\n", fwhm, sld );
   }
   else
      sld = atof(argv[5]);

   fwhm_true = sqrt( fourLogTwo / sld );
   printf( "fwhm_true( sld = %lf ) = %lf\n", sld, fwhm_true );

   p_val = calc_pa( K, sld, area, t_val);

   printf( "P value = %lf = %le\n", p_val, p_val );
   exit(0);
}

double calc_pa(double K, double sld, double area, double t_val)
{
   double n_bar, beta, gamma, prob_k, prob_a;

   n_bar = pow( (double)(2.0 * M_PI), (double)(-3.0/2.0));
   n_bar *= sld * area * t_val;
   n_bar *= exp( (t_val * t_val) / (double)(-2.0));
   
   printf( "n_bar = %lf\n", n_bar);

   gamma = (double)(2 * M_PI) / (t_val * t_val * sld );
   printf( "gamma = %lf\n", gamma);

   beta = (double)(1.0 - (0.47 * sld)) * gamma;
   printf( "beta = %lf\n", beta);

   prob_k = (1.0 + (K / (beta * t_val * t_val)));
   prob_k *= exp(- ((K / beta) + ((K*K)/(2.0 * beta * beta * t_val * t_val))));
   printf( "prob_k = %lf\n", prob_k);

   prob_a = 1.0 - exp( -n_bar * prob_k );

   return( prob_a );
}