/* Analyse der Rechnungen mit IBJparticle für das
 * Rechenbeispiel VDI-3782/1-04
 *
 * 2001-08-15 uj
 */

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

static char  Path[256];
static FILE  *Log;
static char  Buf[4000];
static int   nX=54, nY=39, blockNr=0;    
static float Delta=50;
                          
int calcStat() {
	 int i, j, k, nn;
	 float sigma[200], mean[200], sigma2[200], mean2[200], a, y;
	 float inty[200], inty2[200];  
	 float sum, dist, meanvar, meanstr, sigvar, sigstr, intyvar, intystr, f1, f2, f3; 
	 char *pc;
	 nn=0;
	 for (i=0; i<nX; i++) 
		  sigma[i] = mean[i] = sigma2[i] = mean2[i] = inty[i] = inty2[i] = 0.;
 	while (fgets(Buf, 4000, Log)) {
    if (Buf[4]==':') { 
    	 i = -1;
    	 k = sscanf(Buf, "%d:", &i);
      if (k != 1 || (i-1)%3 != blockNr)  continue;
      nn++;
      for (i=0; i<nX; i++) {
      	 if (i==0) pc = strtok(Buf+5, " \t\r\n"); 
      	 else      pc = strtok(Buf, " \t\r\n"); 
      	 y = -(nY-1)/2*Delta;
      	 sum = dist = 0.;     	 
     		 for (j=0; j<nY; j++) {
      		  if (!pc) {
        		  printf("value #%d in line %d not found!\n", i, nn);
        		  return -1;
        		}
      		  sscanf(pc, "%f", &a);
      		  sum += a;
      		  dist += (y*y + Delta*Delta/12.)*a;
      		  y += Delta;
      		  pc = strtok(NULL, " \t\r\n");    
      		}
      	 f1 = (sum == 0) ? 0. : sqrt(dist/sum);
      	 sigma[i]  += f1;
      	 sigma2[i] += f1*f1;
      	 f2 = sum*Delta;
      	 inty[i]   += f2;
      	 inty2[i]  += f2*f2;
        f3 =  (f1 > 0.) ? f2/(sqrt(2.*3.14159)*f1) : 0;
      	 mean[i]   += f3;
      	 mean2[i]  += f3*f3; 
      	 fgets(Buf, 4000, Log);
      }
    }
		}
	 printf("%3s %4s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s\n-\n",
		  "i", "x(m)", 
		  "Cy(10^6/m)", "Cy-eps", "Cy-low", "Cy-high",
		  "Sy(m)", "Sy-eps", "Sy-low", "Sy-high", 
		  "Cm(10^9/m3)", "Cm-eps", "Cm-low", "Cm-high");
	 for (i=0; i<nX; i++) {
		  meanvar  = sqrt( (mean2[i] - mean[i]*mean[i]/nn)/(nn-1) );
		  meanstr  = meanvar/sqrt(nn);
		  sigvar   = sqrt( (sigma2[i] - sigma[i]*sigma[i]/nn)/(nn-1) );
		  sigstr   = sigvar/sqrt(nn);
		  intyvar  = sqrt( fabs(inty2[i] - inty[i]*inty[i]/nn)/(nn-1) );
		  intystr  = intyvar/sqrt(nn);
		  mean[i]  /= nn;
		  sigma[i] /= nn;
		  inty[i]  /= nn;
		  f1 = (mean[i] > 0.)  ? mean[i]  : 1.;
		  f2 = (sigma[i] > 0.) ? sigma[i] : 1.;
		  f3 = (inty[i] > 0.)  ? inty[i]  : 1.;
		  printf("%3d %4.0f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f\n",
		 	  i, -100 + i*50., 
			   inty[i]*1.e6, 
			   intystr/f3, 
			   inty[i]*1.e6-2*intystr*1.e6, 
			   inty[i]*1.e6+2*intystr*1.e6,
			   sigma[i], 
			   sigstr/f2, 
		   	sigma[i]-2*sigstr, 
			   sigma[i]+2*sigstr,
		   	mean[i]*1.e9, 
		   	meanstr/f1, 
			   (mean[i]-2*meanstr)*1.e9, 
		   	(mean[i]+2*meanstr)*1.e9
			 );
		}
	 return 0;
}
		 
int main( int argc, char *argv[] ) {
  char fn[256];
  int r;
  if (argc < 2) {
    printf("usage: analyse <path> <0:cnc, 1: wet, 2:dry>\n");
    exit(0);
  }
  sscanf(argv[2], "%d", &blockNr);
  sprintf(fn, "%s\\%s",  argv[1], "prtlog.txt");
  Log = fopen(fn, "rb");
  if (!Log) {
    printf("can't read file %s!\n", fn);
    exit(1);
  }
  r = calcStat();
  fclose(Log);
  if (r) {
    printf("error while analysing!\n");
    exit(2);
  }
  return 0;
}
  
