/* plot udrnrt.c  3 choices  uniformly distributed=1
 *                           gauss distributed=2  max 3 sigma
 *                           gausse distributed=3  max 10 sigma?
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include "udrnrt.h"

int main(int argc, char argv[])
{
  double mean, sigma, *v;
  double maxv, minv, meanv, sigmav;
  double sum, sumsq, stddev;
  int i, j, n, m, kind;
  int stat;
  FILE * outp;
  double binv[100];
  int binn[100];
  double aval, bval;

  printf("plot_udrnrt.c , enter number of samples:\n");
  scanf("%d",&n);
  printf("you have selected %d samples\n",n);
 
  printf("enter 1 for uniformly distributed:\n");
  printf("or 2 for gauss distributed:\n");
  printf("or 3 for gausse distributed:\n");
  printf("or 4 for sine wave:\n");
  scanf("%d",&kind);
  if(kind==1) printf("you have selected uniformly distributed\n",n);
  if(kind==2)
  {
    printf("you have selected gauss distributed\n",n);
    printf("enter mean sigma:\n");
    scanf("%lf %lf",&mean, &sigma);
    printf("you have selected mean=%f, sigma=%f\n", mean, sigma);
  }
  if(kind==3)
  {
    printf("you have selected gausse distributed\n",n);
    printf("enter mean sigma:\n");
    scanf("%lf %lf",&mean, &sigma);
    printf("you have selected mean=%f, sigma=%f\n", mean, sigma);
  }
  if(kind==4) printf("you have selected a sine wave\n",n);

  v=(double *)malloc(sizeof(double)*n);
  if(v==NULL){printf("malloc failed to allocate %d doubles\n",n); return 1;}

  if(kind==1)
  {
    for(i=0; i<n; i++) v[i] = udrnrt();
    meanv = 0.0;
    sumsq = 0.0;
    maxv = -99999.99;
    minv = 99999.99;
    for(i=0; i<n; i++)
    { 
      meanv = meanv + v[i];
      sumsq = sumsq + v[i]*v[i];
      if(v[i]>maxv) maxv = v[i];
      if(v[i]<minv) minv = v[i];  
    }
    sum = meanv;
    stddev = sqrt((sumsq-sum*sum/n)/(n-1));
    meanv = meanv/(double)n;
    printf("computed min=%f, max=%f, mean=%f, sigma=%f \n",
            minv, maxv, meanv, stddev);
    printf("total power=%g, avg power per sample=%f \n",
	    sumsq, sumsq/(double)n);
  } /* end kind==1 */

  if(kind==2)
  {
    for(i=0; i<n; i++) v[i] = gauss(mean, sigma);
    meanv = 0.0;
    sumsq = 0.0;
    maxv = -99999.99;
    minv = 99999.99;
    for(i=0; i<n; i++)
    { 
      meanv = meanv + v[i];
      sumsq = sumsq + v[i]*v[i];
      if(v[i]>maxv) maxv = v[i];
      if(v[i]<minv) minv = v[i];  
    }
    sum = meanv;
    stddev = sqrt((sumsq-sum*sum/(double)n)/(double)(n-1));
    meanv = meanv/(double)n;
    printf("computed min=%f, max=%f, mean=%f, sigma=%f \n",
            minv, maxv, meanv, stddev);
    printf("total power=%g, avg power per sample=%f \n",
	    sumsq, sumsq/(double)n);
  } /* end kind==2 */

  if(kind==3)
  {
    for(i=0; i<n; i++) v[i] = gausse(mean, sigma);
    meanv = 0.0;
    sumsq = 0.0;
    maxv = -99999.99;
    minv = 99999.99;
    for(i=0; i<n; i++)
    { 
      meanv = meanv + v[i];
      sumsq = sumsq + v[i]*v[i];
      if(v[i]>maxv) maxv = v[i];
      if(v[i]<minv) minv = v[i];  
    }
    sum = meanv;
    stddev = sqrt((sumsq-sum*sum/(double)n)/(double)(n-1));
    meanv = meanv/(double)n;
    printf("computed min=%f, max=%f, mean=%f, sigma=%f \n",
            minv, maxv, meanv, stddev);
    printf("total power=%g, avg power per sample=%f \n",
	    sumsq, sumsq/(double)n);
  } /* end kind==3 */

  if(kind==4)
  {
    for(i=0; i<n; i++) v[i] = sin(2.0*M_PI*(double)i/(double)n);
    meanv = 0.0;
    sumsq = 0.0;
    maxv = -99999.99;
    minv = 99999.99;
    for(i=0; i<n; i++)
    { 
      meanv = meanv + v[i];
      sumsq = sumsq + v[i]*v[i];
      if(v[i]>maxv) maxv = v[i];
      if(v[i]<minv) minv = v[i];  
    }
    sum = meanv;
    stddev = sqrt((sumsq-sum*sum/n)/(n-1));
    meanv = meanv/(double)n;
    printf("computed min=%f, max=%f, mean=%f, sigma=%f \n",
            minv, maxv, meanv, stddev);
    printf("total power=%g, avg power per sample=%f \n",
	    sumsq, sumsq/(double)n);
  } /* end kind==4 */

  m = 101;
  if(m>n) m = n;
  bval = (maxv-minv)/(double)(m-1);
  aval = minv-bval/2.0; /* center bin */
  for(i=0; i<m-1; i++)
  {
    binv[i] = aval;
    binn[i] = 0;
    for(j=0; j<n; j++) if(v[j]>=aval && v[j]<aval+bval) binn[i]++;
    aval = aval + bval; 
  }

  outp = fopen("udrnrt.dat", "w");
  if(outp==NULL) printf("could not open udrnrt.dat for write \n");
  for(i=1; i<m-1; i++) fprintf(outp, "%f %d \n", binv[i], binn[i]);
  fprintf(outp, " \n");
  fclose(outp);

  printf("plot_udrnrt.c about to run gnuplot on udrnrt.plot \n");
  fflush(stdout);
  fclose(stdout);
  stat = execl("/usr/bin/gnuplot","gnuplot","-persist","udrnrt.plot", NULL);
  /* only returns upon error, no return on sucess */
  printf("execl failed on gnuplot, stat=%d \n", stat);
  return 0;
} /* end plot_udrnrt.c */