/************************************************************************/
/*  		yrms : is a part of Yao utilities and tools		*/
/************************************************************************/

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

#define BUFSIZE	80

int main (int argc, char *argv[])
{    /* . Syntaxe: yrms f1 f2 zspace [f3] */
     /* . Purpose: calculer l'erreur rms entre 2 trajectoires d'espace
		 dont les valeurs sont respectvement dans les fichiers
		 f1 et f2 passes en parametres qui doivent contenir
		 le meme nombre d'elements qui doit d'ailleurs etre
		 un multiple de zspace qui est la taille de l'espace
		 le resultat, {pour chaque pas d'espace, zspace} est
		 ecrit dans f3 s'il est mentionne sinon sur la sortie
		 standard.
		 la somme des rms est affichee a la fin
	. Warning: la detection de ligne blanche ou vide n'est pas geree
	. Compilation :	cc -Wall -o yrms yrms.c -lm			*/

	FILE    *f1, *f2, *f3;
	char	buf1[BUFSIZE+1], buf2[BUFSIZE+1];
	int	zspace, rstelt, nbt;
	double	x1, x2, errsumt, errsumallt, rms, sumrms;

	if (argc<4 || argc > 5)
	{	printf("syntaxe error !!! : yrms f1 f2 zspace [f3]\n");
		exit(0);
	}

	if ((f1 = fopen(argv[1], "r")) <= 0)
	{	printf ("file1 %s not found !!!\n", argv[1]);
		exit(0);
	}
	if ((f2 = fopen(argv[2], "r")) <= 0)
	{	printf ("file2 %s not found !!!\n", argv[2]);
		exit(0);
	}
	if (argc==5)
	{	if ((f3 = fopen(argv[4], "w")) <= 0)
		{	printf ("failed while opening file3 %s !!!\n", argv[3]);
			exit(0);
		}
	}
	else
		f3 = stdout;

	zspace = atoi(argv[3]);
	if (zspace<=0)
	{	printf ("zspace must be a integer greater than zero !!!\n");
		exit(0);
	}

	rstelt=0; nbt=0; errsumt=0.0; errsumallt=0.0, sumrms=0.0;
	while (1) /* boucle de */
	{	/* lecture paralle des 2 fichiers */
		if (fgets (buf1, BUFSIZE, f1) == NULL)
		{  /* fin du 1er fichier => fin du 2eme fichier */
		   if (fgets (buf2, BUFSIZE, f2) != NULL)
		   {  printf("erratum 1: the files have not the same size !!!\n");
		      exit(0);
		   }

		   /* fin: le nombre d'elt devait etre un multiple de zspace */
		   if (rstelt!=0)
		   {  printf("erratum 3: number of element is not o multiple of zspace\n");
		      exit(0);
		   }

		   /* fin fin */
		   break;
		}

		/* on doit pouvoir lire le 2eme fichier */
		if (fgets (buf2, BUFSIZE, f2) == NULL)
		{  printf("erratum 2: the files have not the same size !!!\n");
		   exit(0);
		}

		x1 = atof(buf1);	x2 = atof(buf2);
		//printf("x1=%f x2=%f \n", x1, x2);

		errsumt += ((x1-x2)*(x1-x2));
		errsumallt += errsumt;

		rstelt = (rstelt+1)%zspace;
		if (rstelt==0)
		{	//printf ("errsumt=%f \n", errsumt);
			rms = sqrt(errsumt/zspace);
			fprintf(f3, "%f\n",  rms);
			errsumt = 0.0;
			++nbt;
			sumrms += rms;
		}

	} /* fin du while de la boucle de lecture parallele */
	printf("SUMt(rms) = %30.20e\n", sumrms);

	/* rms sur tout la trajectoire */
	rms = sqrt(errsumallt/(zspace*nbt));
	printf("RMS(trajectoire::%i pas de temps) = %30.20e\n", nbt, rms);

	printf("yrms: process end normaly\n");
	return (0);
}