#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>
#include"mt_mogy.h"

/* a tarolashoz hasznalt tombok meretei */
#define IMAX 20000  /* oldalankent (bal es jobb) legfeljebb ennyi reszecsket
		  tudunk kovetni. Legyen legalabb 20*T, hogy a kezdetben
		  10*T-nel kozelebb levo reszecskek (varhatoan 10*T darab)
		  biztosan beleferjenek. */
#define KMAX 10001  /* legfeljebb ennyi diszkret idopontot tudunk szamontartani.
                  Legyen legalabb T/dt . */


/* bemeneti es kimeneti fileok */
FILE *input,*output,*output2;
char outfile[80];   /* kimeneti file neve */

/* a szimulaciot meghatarozo globalis "konstansok" - minden szimulacioban (a bemeneti file minden sorara)
   csak egyszer kapnak erteket */
int uzemmod;  /* A "c" 'rugoallando' megvalasztasanak modja. Ha 1, akkor "c"-t a termeszetes modon
                   valasztjuk. Ha 2, akkor ennek "A"-szorosának. Ha 3, akkor "c" elore megadott
                   konstans. Ha 4, akkor c=W*d, ahol "W" elore adott konstans */
int teszt;    /* ha 1, akkor teszt gyanant kiiratjuk az osszes reszecske helyet az ido fv-eben */
              /* ez esetben celszeru 1-elemu mintat venni, kulonben sokaig tart */
              /* ha 2, akkor kiiratjuk a ket reszecske tavolsagat minden szimulacio elejen es vegen */
long int N1;
long int N2;   /* a mintavetel elemszama a ketto szorzata*/
double T;     /* a szimulacio idotartama */
double dt;      /* ennyi idoegysegenkent szamoljuk a szorasokat */

/* a modositott uzemmodokhoz szukseges bemeno adatok */
double parameter; /* ezt fogjuk beolvasni, es ertkeul adjuk c-nek, A-nak vagy W-nek az uzammodtol fuggoen */
double c;     /* a ket kituntetett reszecsket egymastol eltaszito rugot jellemzo szam - */
                  /* ez a szimulacio soran konstans, de a kezdeti ertektol fugg, az "uzemmod"-tol */
                  /* fuggo modon */
double A;     /* ennyiszeresre noveljuk "c"-t a 2-es uzemmodban */
double W;     /* a 4-es uzemmodhoz kell */

/* a fazister koordinatai */
double x;           /* a kituntetett reszecskek tavolsaganak fele (helyuk a tkp-i rendszerben) */
double v;           /* a kituntetett reszecskek relativ sebessegenek fele (sebesseguk tkp-tol) */
double nbal,njobb;  /* ennyi egyeb reszecske sorsat kovetjuk */
double xbal[IMAX],xjobb[IMAX],vbal[IMAX],vjobb[IMAX]; /* az egyeb reszecskek helye es sebessege */

/* a szimulacioval kovetendo mennyisegek */
double t;           /* ennyi ido telt el az eleje ota */
double X;           /* a ket kozepso reszecske helyenek atlaga (tomegkp helye) */
double V;           /* az imenti tomegkp sebessege */

/* az eredmenyt tarolo tombok */
double sX_N[KMAX],sX2_N[KMAX],sV_N[KMAX],sV2_N[KMAX];
       /* a mintabol szamolt osszegek az ido fuggvenyeben dt egysegenkent */
double sX[KMAX],sX2[KMAX],sV[KMAX],sV2[KMAX];
       /* az elozoek reszletosszegei a mintavetel belso ciklusa soran */

double akk[KMAX],mkk[KMAX],hibak[KMAX],rkk[KMAX];
  /* tesztelesi cellal minden ponthoz megjegyezzuk az illesztett egyenest */

/* tovabbi globalis konstansok */
double TT,kmaxdt;   /* T-hez kozeli ertekek, hogy a diszkretizalas ne okozzon gondot */
double c1;          /* 12.0/(dt*dt) , az egyenesilleszteshez kell */

/* az ido meresehez */
time_t ido_volt,ido_most;

/* ciklusvaltozok */
long int i,j11,j2,k;

/* logikai */
int jobbrol;

/* segedvaltozok: */
double xmax, tmax;    /* globalis konstansok, ket irrealisan nagy ertek (10*T)  */
double E,ts,vlim;     /* a kituntetett reszecske palyajanak jellemzoi */
long int hanyadik;
double tv,xv,vv,dV,Dt;   /* a kovetkezo utkozes adatai */
double kdt,Xk,xk;        /* az eredmenyek diszkret ideju feljegyzesehez */
long int kmax;           /* ennyi diszkret ideju lepes fel bele */
long int n;           /* mintevetel elemszama (n+1) az egyenes-illesztesnel */
double yi,sy,sy2,sxy,ck,rkorr,rkorr_max,hiba,hiba_min; /* egyenes-illeszteshez */
double ak,mk,tk,rk;    /* az illesztett egyenes parameterei */
double sak,sak2,smk,smk2,stk,stk2,srk,srk2; /* statisztikak az illesztett egyenesekrol*/
double v1,v2,xi,vi,Einv,Vdt,c2,c2per4,c2per2;  /* ideiglenes tarolasra */
double d,ds,dmax,fmax,ccp2; /* d-nek adott c melletti felteteles eloszlasa generalasahoz */
int elfogadom;			/* ugyszinten */



/* a veletlengenerator inicializalasahoz */
int randmode; /* ez adja meg az inicializalas modjat */
unsigned long int randseed; /* ez lesz az inicializalo mag */




main(int argc, char *argv[]){

/* parameter-file megnyitasa */
        if (argc!=3) {
                printf("A program hasznalata:\n a.out <parameter file> <eredmeny file>\n");
                return -1;
        	}
                                                                                
        input=fopen(argv[1],"rt");      /* Bemeneti file megnyitasa */
        output2=fopen(argv[2],"w");      /* Kimeneti file megnyitasa */

                                                                                
        if (input==NULL) {
                printf ("Hiba a bemeneti file beolvasasa kozben.\n");
                return -1;
        	}
        if (output2==NULL) {
                printf ("Hiba a fo kimeneti file megnyitasa kozben.\n");
                return -1;
        	}
        /* A parameterek beolvasasa a bemeneti filebol */

/*printf("A bemeneti filet sikerult megnyitni.\n");*/

/* a fo kimeneti file fejlecenek kiiratasa */
fprintf(output2,"A feldolgozott bemeneti file neve: %s\n\n",argv[1]);
fprintf(output2,"Az alabbi tablazat minden soraban a bemeneti file megfelelo soraban ");
fprintf(output2,"megadott szimulacio eredmenyei vannak.\n\n");
fprintf(output2,"Jelmagyarazat:\n\n");
fprintf(output2,"Futasi parameterek:\n");
fprintf(output2,"outfile: a szimulacio reszletes eredmenyeit tartalmazo file neve\n");
fprintf(output2,"teszt: nemnulla, ha teszt-futtatas a trajektoriak reszletes kiirasaval\n");
fprintf(output2,"randmode: a veletlengenerator inicializalasanak modja\n");
fprintf(output2,"randseed: a veletlengenerator inicializalasahoz hasznalt mag\n");
fprintf(output2,"N1: a reszmintak szama\n");
fprintf(output2,"N2: a reszmintak elemszama - igy az egyesitett minta N1*N2 elemu\n");
fprintf(output2,"T: a szimulacio idotartama\n");
fprintf(output2,"dt: a meresi idopontok tavolsaga\n");
fprintf(output2,"uzemmod: a potencial konstansa megvalasztasanak modja (a parameter jelentese)\n");
fprintf(output2,"parameter: a potencial megvalasztasahoz hasznalt parameter\n");
fprintf(output2,"m: az egyesitett mintabol szamolt meredekseg\n\n");
fprintf(output2,"Szimulacios eredmenyek:\n");
fprintf(output2,"ezek mind a kituntetett reszecske helyenek szorasnegyzetet az ido fuggvenyeben\n");
fprintf(output2,"megado fuggveny aszimptotajanak (mint egyenesnek) az adatai\n");
fprintf(output2,"m: meredekseg az egyesitett minta alapjan szamolva\n");
fprintf(output2,"a: tangelymetszet az egyesitett minta alapjan szamolva\n");
fprintf(output2,"t_min: az aszimptota illesztesehez hasznalt legkisebb idopont - innentol jo az illeszkedes\n");
fprintf(output2,"r^2: az egyenes-illesztes soran a tapasztalati korrelacios egyutthato negyzete\n");
fprintf(output2,"m_atl, a_atl, t_min_atl, r^2_atl: ugyanezen mennyisegekre a reszmintakbol szamolt");
fprintf(output2,"becslesek atlaga\n");
fprintf(output2,"m_hiba, a_hiba, t_min_hiba, r^2_hiba: becsult hiba a reszmintakbol keszult");
fprintf(output2,"becslesekbol keszult N1 elemu statisztika alapjan\n");
fprintf(output2,"szim_ido: a szimulacio idotartama masodpercben\n");
fprintf(output2,"ido_0: a szimulacio kezdetenek idopontja\n\n");

fprintf(output2,"outfile\tteszt\trandmode\trandseed\tN1\tN2\tT\tdt\t");
fprintf(output2,"uzemmod\tparameter\tm\ta\tt_min\tr^2\t");
fprintf(output2,"m_atl\tm_hiba\ta_atl\ta_hiba\t");
fprintf(output2,"t_min_atl\tt_min_hiba\tr^2_atl\tr^2_hiba\tszim_ido\tido_0\t\n");

/*printf("Tul vagyok a fejlec kiirasan.\n");*/



/* egy legeslegkulso ciklus a bemeneti file sorainak feldolgozasara */
while (!feof(input)){

        fscanf(input,"outfile=%s\tteszt=%d\trandmode=%d\trandseed=%ld\t",outfile,&teszt,&randmode,&randseed);
        fscanf(input,"N1=%ld\tN2=%ld\tT=%le\tdt=%le\tuzemmod=%d\tparameter=%le\n",&N1,&N2,&T,&dt,&uzemmod,&parameter);

/*
printf("Egy sort, ugy tunik, beolvastam. Ugy tudom, hogy\n");
printf("outfile=%s\n",outfile);
printf("randmode=%lu\nrandseed=%lu\nteszt=%d\n",randmode,randseed,teszt);
printf("T=%f\ndt=%f\nuzemmod=%d\nparameter=%f\n\n",T,dt,uzemmod,parameter);
*/

if((T>KMAX*dt)||(2*T>IMAX)){
	fprintf(stderr,"Hiba: T vagy T/dt tul nagy.\n");
	continue;
	}
/* printf("randmode=%d\trandseed=%ld\n",randmode,randseed); */

/* printf("Nekilatok a generatort inicializalni.\n"); */

/* kezdo idopont felirasa */
ido_volt=time(NULL);

/* veletlengenerator inicializalasa */
if((randmode==0)||(randseed==0)){
/* printf("ido alapjan fogom csinalni.\n"); */
	randseed=time(NULL); 
/* printf("idot megkerdeztem, es most randseed=%lu\n",randseed); */
	randseed*=630360016;
	init_genrand(randseed);
	}
else if(randmode==1) init_genrand(randseed);

/* ha randmode==0, a rendszeridot hasznaljuk; ha randmode=1, akkor a megadott seed-et */
/* barmely mas esetben nem inicializalja ujra a generatort az egyes szimulaciok kozott. */
/* Ha egyszer se inicializaljuk, az elejen az alapertelmezes */ 
/* (randseed=5489, a mostani mt.h-ban) lep eletbe.   */

/*printf("inicializaltam a randomgeneratort.\n");*/



/* ertekadas globalis konstansoknak */
A=c=W=parameter; /* parameter ugyis csak egy van, majd elvalik, melyiket hasznaljuk */
TT=T; /* TT-be mentjuk el a T erteket, hogy a vegen ki tudjuk iratni, de a valosagban
         egy kicsivel mast (tobbnyire nagyobbat) hasznalunk helyette, kerekitesi
         hibak elkerulese vegett */
xmax=10*T; tmax=10*T; kmax=(T+0.001*dt)/dt; kmaxdt=kmax*dt; T=kmaxdt+0.5*dt;
if(uzemmod==3){ds=1.5*pow(c,2.0/3.0);dmax=ds+1.0;fmax=exp(-ds);ccp2=c*c*0.5;}
	/* ezek csak a d-nek adott c melletti felteteles eloszlasahoz kellenek */
c1=12.0/(dt*dt); /* csak az egyenes-illeszteshez kell majd */


/* a statisztikat tarolo valtozok nullazasa */
sak=sak2=smk=smk2=stk=stk2=srk=srk2=0.0;

/* a statisztikakat tarolo tombok nullazasa */
for(k=0;k<=kmax;k++) { sX_N[k]=0.0; sX2_N[k]=0.0; sV_N[k]=0.0; sV2_N[k]=0.0; }
/*printf("kmax=%d\n",kmax);*/

/* N=N1*N2 elemu mintat veszunk */
for(j11=1;j11<=N1;j11++){        /* mintavetel kulso ciklusanak eleje */
	/*printf("%d\n",j11);*/

	/* a statisztikakat tarolo tombok nullazasa */
	for(k=0;k<=kmax;k++) { sX[k]=0.0; sX2[k]=0.0; sV[k]=0.0; sV2[k]=0.0; }



   for(j2=1;j2<=N2;j2++){        /* mintavetel belso ciklusanak eleje */

	/*************** kezdofeltetel generalasa ****************/
	/* ha uzemmod=3, akkor a kezdeti tavolsagot az adott c melletti felteteles eoszlassal veszem */
	if(uzemmod==3){
		elfogadom=0;
		while(!elfogadom){
			d=dmax*rrand();
			if(d<=ds) {if(rrand()*fmax<=exp(-d-ccp2/(d*d))) elfogadom=1;} /*elso eset*/
			else{   /* masodik eset */
				d=ds-log(d-ds); /* x-xs egyenletes (0,1]-en, igy -log(d-ds) exponencialis */
				if(rrand()<=exp(-ccp2/(d*d)))	elfogadom=1;
				}
			}
		x=0.5*d;
		}
	else x=-0.5*(log(rrand())+log(rrand()));  /* a kit.ek tavolsaga 2x, ami ket fuggetlen exp. osszege */
	v1=nrand(); v2=nrand();  /* a kituntetettek sebessege fuggetlen normalis */
	V=0.5*(v1+v2);  /* a tomegkozeppont sebessege */
	v=0.5*(v2-v1);  /* tomegkozepponthoz viszonyitott sebesseg */
	i=1; xi=x-log(rrand());
	while(xi<xmax){xjobb[i]=xi; vjobb[i]=nrand()-V; xi-=log(rrand()); i++;} /* jobboldali reszecskek */
	njobb=i-1;
	i=1; xi=x-log(rrand());
	while(xi<xmax){xbal[i]=xi; vbal[i]=nrand()+V; xi-=log(rrand()); i++;} /* baloldali reszecskek */
	nbal=i-1;
	switch(uzemmod){       /* a rugoallandot az uzemmodtol fuggoen valasztjuk */
		case 1: c=2.0*x*nrand(); break;      /* a termeszetes valasztas */
		case 2: c=A*2.0*x*nrand(); break;    /* a termeszetes valasztas A-szorosa */
		case 3: break;                       /* c az elore adott konstans marad */
		case 4: c=2.0*x*W; break;            /* c=d*W, ahol W az elore adott konstans */
		}
	t=0.0; X=0.0;
	k=1; kdt=dt; /* k indexeli az sX, sX2, sV es sV2 tomboket, kdt a diszkret ido */

	/* statisztikak elso elemei */
	/* sX2[0]+=0.0; */
	sV[0]+=V;
	sV2[0]+=V*V;

	/* segedvaltozok ertekadasa */
	c2=c*c; c2per4=0.25*c*c; c2per2=0.5*c2;

/* tesztelesi cellal kiirhatjuk a kezdoallapotot: */
if(teszt==1){
  printf("%f\t",0.0);
  printf("%f\t",x);printf("%f\t",-x);
  for(i=1;i<=njobb;i++) printf("%f\t",xjobb[i]);
  for(i=1;i<=nbal;i++) printf("%f\t",-xbal[i]);
  printf("\n");
  }
if(teszt==2) printf("%f\t",2.0*x);
/* eddig a tesztelesi celu kiiratas */




	/* a szimulacio addig megy, amig lesz, aki idon belul utkozik */
	jobbrol=1;
	while(jobbrol!=0)
		{                        /* szimulacio ciklusanak eleje */

		/*************** szimulacios lepes *****************/

		/* kiszamoljuk a kituntetett reszecske hiperbolapalyajanak parametereit */
		E=2*v*v+c2per4/(x*x); Einv=1.0/E; vlim=sqrt(0.5*E);
		ts=sqrt(2*x*x*E-c2per2)*Einv; if(v>0) ts=-ts;  /* ez a hiperbola t^* parametere */

		/* eloszor megkeressuk, hol es mikor lesz a legkozelebbi utkozes: xv-ben Dt ido mulva */
		/* illetve melyik oldalrol (jobbrol, ha jobbrol=1 es balrol, ha jobbrol=-1) es hanyadik*/
		/* fog eloszor utkozni */
		xv=xmax; Dt=0.0;   /* eloszor irrealisan nagyot tippelunk, ennel kell jobbat talalni */
		jobbrol=0; /* ha ennyi marad, egy reszecske se fog utkozni, vege is lesz a ciklusnak */

		/* megkeressuk a legjobb jobboldali reszecsket */
		for(i=1;i<=njobb;i++) if((vjobb[i]<vlim)&&(xjobb[i]+vjobb[i]*Dt<xv))
									/* akkor jobb az eddiginel, ha ez igaz */
			{
			vi=vjobb[i]; xi=xjobb[i];
			Dt=(2*xi*vi+E*ts+sqrt(2*E*(vi*ts+xi)*(vi*ts+xi)+vi*vi*c2*Einv-c2per2))/(E-2*vi*vi);
			xv=xi+Dt*vi;    /* es ez esetben ekkor es itt fog utkozni */
			hanyadik=i; jobbrol=1;
			}

		/* megnezzuk, bal oldalon van-e jobb */
		for(i=1;i<=nbal;i++) if((vbal[i]<vlim)&&(xbal[i]+vbal[i]*Dt<xv))
									/* akkor jobb az eddiginel, ha ez igaz */
			{
			vi=vbal[i]; xi=xbal[i];
			Dt=(2*xi*vi+E*ts+sqrt(2*E*(vi*ts+xi)*(vi*ts+xi)+vi*vi*c2*Einv-c2per2))/(E-2*vi*vi);
			xv=xi+Dt*vi;    /* es ez esetben ekkor es itt fog utkozni */
			hanyadik=i; jobbrol=-1;
			}

		/* ha nincs tobb utkozes, egy hosszu idot repulunk */
		if(jobbrol==0) Dt=tmax;

		/* feljegyezzuk a tomegkozeppont helyet es sebesseget a kovetkezo utkozesig */
		tv=t+Dt;
		if(tv>T) {tv=T; jobbrol=0;}  /* ha T-ig nincs utkozes, mindegy is, merrol lenne, ciklus vege */
		/* k ill. kdt tarolja az elso olyan diszkret idot, amit meg korabban nem ertunk el */
		/* az elso diszkret pontot kulon kell kezelni, mert addig a dt-nek csak tort resze telik el */
		if(kdt<=tv){
			Xk=X+V*(kdt-t);
			sX[k]+=Xk;
			sV[k]+=V;
			sX2[k]+=Xk*Xk;
			sV2[k]+=V*V;
/* tesztelesi cellal kiirhatjuk a trajektoriat: */
if(teszt==1){
  printf("%f\t",kdt);
  xk=sqrt(c2per4*Einv+0.5*E*(kdt-t-ts)*(kdt-t-ts));
  printf("%f\t",Xk+xk);printf("%f\t",Xk-xk);
  for(i=1;i<=njobb;i++) printf("%f\t",Xk+(xjobb[i]+(kdt-t)*vjobb[i]));
  for(i=1;i<=nbal;i++) printf("%f\t",Xk-(xbal[i]+(kdt-t)*vbal[i]));
  printf("\n");
  }
/* eddig a tesztelesi celu kiiratas */
			k++; kdt+=dt;
			Vdt=V*dt;
			}
		/* a tobbi diszkret lepes hossza mar dt */
		while(kdt<=tv){
			Xk+=Vdt;
			sX[k]+=Xk;
			sV[k]+=V;
			sX2[k]+=Xk*Xk;
			sV2[k]+=V*V;
/* tesztelesi cellal kiirhatjuk a trajektoriat: */
if(teszt==1){
  printf("%f\t",kdt);
  xk=sqrt(c2per4*Einv+0.5*E*(kdt-t-ts)*(kdt-t-ts));
  printf("%f\t",Xk+xk);printf("%f\t",Xk-xk);
  for(i=1;i<=njobb;i++) printf("%f\t",Xk+(xjobb[i]+(kdt-t)*vjobb[i]));
  for(i=1;i<=nbal;i++) printf("%f\t",Xk-(xbal[i]+(kdt-t)*vbal[i]));
  printf("\n");
  }
/* eddig a tesztelesi celu kiiratas */
			k++; kdt+=dt;
			}

		/* megfeleloen modositjuk a rendszer allapotat */
		if(jobbrol==1){
			x=xv;
			vv=sqrt(0.5*E-0.5*c2per4/(xv*xv)); if(Dt<ts) vv=-vv;
								/* a kit. reszecskek sebessege az utkozeskor */
			dV=0.5*(vi-vv); /* a tomegkozeppont sebessegenek megvaltozasa az utkozeskor */
			for(i=1;i<=njobb;i++){
				xjobb[i]+=Dt*vjobb[i];
				vjobb[i]-=dV; /* atteres az uj tomegkozepponti rendszerre */
				}
			for(i=1;i<=nbal;i++){
				xbal[i]+=Dt*vbal[i];
				vbal[i]+=dV; /* atteres az uj tomegkozepponti rendszerre */
				}
			v=vv+dV; vjobb[hanyadik]=vv-dV;
			t=tv;
			X+=Dt*V;
			V+=dV;
			}
		if(jobbrol==-1){
			x=xv;
			vv=sqrt(0.5*E-0.5*c2per4/(xv*xv)); if(Dt<ts) vv=-vv;
								/* a kit. reszecskek sebessege az utkozeskor */
			dV=0.5*(vv-vi); /* a tomegkozeppont sebessegenek megvaltozasa az utkozeskor */
			for(i=1;i<=njobb;i++){
				xjobb[i]+=Dt*vjobb[i];
				vjobb[i]-=dV; /* atteres az uj tomegkozepponti rendszerre */
				}
			for(i=1;i<=nbal;i++){
				xbal[i]+=Dt*vbal[i];
				vbal[i]+=dV; /* atteres az uj tomegkozepponti rendszerre */
				}
			v=vv-dV; vbal[hanyadik]=vv+dV;
			t=tv;
			X+=Dt*V;
			V+=dV;
			}

		/* szimulacios lepes vege */

		}  /* szimulacio ciklusanak vege */

if(teszt==2) printf("%f\n", 2.0*sqrt(c2per4*Einv+0.5*E*(kmaxdt-t-ts)*(kmaxdt-t-ts)) );

	}  /* mintavetel belso ciklusanak vege */

   /* a reszletosszegeket hozzaadjuk a statisztikat tartalmazo tombokhoz */
   for(k=0;k<=kmax;k++) {
	sX_N[k]+=sX[k];
	sX2_N[k]+=sX2[k];
	sV_N[k]+=sV[k];
	sV2_N[k]+=sV2[k];
	}

   /* A kituntetett reszecskek tomegkozeppontja helyenek a belso ciklusban szamolt
      szorasnegyzetere - vagyis sX2[i]/N2 -re, mint az ido fuggvenyere, egyenest
      illesztunk. Ezt ugy csinaljuk, hogy megkeressuk, hany pontot kell venni a
      kapott gorbe vegerol, hogy a korrelacios egyutthato korrigalt becslese
      maximalis legyen, es ennyi pontra illesztunk egyenest a legkisebb
      negyzetek modszerevel. */
   yi=sX2[kmax]/N2; sy=yi; sy2=yi*yi; sxy=kmaxdt*yi; /* az osszeg-valtozok kezdoerteke */
   hiba=hiba_min=1.0;  /* az ediggi legkisebb (1-rk)/(n-1) */
   for(k=kmax-1,kdt=k*dt,n=1;k>=0;k--,kdt-=dt,n++) {    /* a minta (n+1) elemu!! */
	yi=sX2[k]/N2; sy+=yi; sy2+=yi*yi; sxy+=kdt*yi;
	ck=(sxy-0.5*(kdt+kmaxdt)*sy)/(n+1);   /* kovariancia becslese */
	rk=c1*ck*ck/((n+2)*(sy2-sy*sy/(n+1)))*(n+1)/n;  /* korrigalatlan korrelacios egyutthato */
	if(n>1) hiba=(1.0-rk)/(n-1);
	if((n>10)&&(hiba<hiba_min)) {
		hiba_min=hiba;   /* az eddigi legkisebb hiba-jellemzo */
		mk=c1*ck/(n*(n+2));   /* meredekseg becslese */
		ak=sy/(n+1)-0.5*(kdt+kmaxdt)*mk;   /* a tengelymetszet becslese */
		tk=kdt;   /* az eddigi legjobb illeszkedes kezdopontja */
		}
	}

   /* az illesztett egyenes parametereit bevesszuk az errol keszult statisztikaba */
   sak+=ak; sak2+=ak*ak;
   smk+=mk; smk2+=mk*mk;
   stk+=tk; stk2+=tk*tk;
   srk+=rk; srk2+=rk*rk;


   }  /* mintavetel kulso ciklusanak vege */


/* Most illesztunk egy egyenest az X egyesitett minta alapjan becsult szorasnegyzetere,
   vagyis sX2_N[i]/(N1*N2) -re   - ugyanugy, mint a resz-mintak alapjan tettuk */

   yi=sX2_N[kmax]/N2/N1; sy=yi; sy2=yi*yi; sxy=kmaxdt*yi; /* az osszeg-valtozok kezdoerteke */
   hiba=hiba_min=1.0;  /* az ediggi legkisebb (1-rk)/(n-1) */
   for(k=kmax-1,kdt=k*dt,n=1;k>=0;k--,kdt-=dt,n++) {    /* a minta (n+1) elemu!! */
	yi=sX2_N[k]/N2/N1; sy+=yi; sy2+=yi*yi; sxy+=kdt*yi;
	ck=(sxy-0.5*(kdt+kmaxdt)*sy)/(n+1);   /* kovariancia becslese */
	rk=c1*ck*ck/((n+2)*(sy2-sy*sy/(n+1)))*(n+1)/n;  /* korrigalatlan korrelacios egyutthato */
	if(n>1) hiba=(1.0-rk)/(n-1);
      /* a kovetkezo 5 sor tesztelesi cellal megjegyzi az osszes becslest */
	mkk[k]=c1*ck/(n*(n+2));   /* meredekseg becslese */
	akk[k]=sy/(n+1)-0.5*(kdt+kmaxdt)*mkk[k];   /* a tengelymetszet becslese */
	rkk[k]=rk; /* korrigalatlan becsles a korrelacios egyutthato
		                  negyzetere. Ez lesz 1-hez kozel, ha jo az illeszkedes */
	hibak[k]=hiba;
	if((n>10)&&(hiba<hiba_min)) {
		hiba_min=hiba;   /* az eddigi legkisebb hiba-jellemzo */
		mk=c1*ck/(n*(n+2));   /* meredekseg becslese */
		ak=sy/(n+1)-0.5*(kdt+kmaxdt)*mk;   /* a tengelymetszet becslese */
		tk=kdt;   /* az eddigi legjobb illeszkedes kezdopontja */
		}
	}

/* a vegen az ido feljegyzese */
ido_most=time(NULL);

/* a vegeredmenyek kiszamolasa, es egyuttal kiirasa */

output=fopen(outfile,"w");

fprintf(output,"futas idotartama masodpercben: %f\n",difftime(ido_most,ido_volt));
fprintf(output,"futas kezdete: %s\n",ctime(&ido_volt));

fprintf(output,"futasi parameterek:\n");
fprintf(output,"uzemmod=%d\n",uzemmod);
fprintf(output,"teszt=%d\n",teszt);
fprintf(output,"N1=%d\n",N1);
fprintf(output,"N2=%d\n",N2);
fprintf(output,"T=%f\n",TT);
fprintf(output,"dt=%f\n",dt);
fprintf(output,"parameter=%f\n\n",parameter);


fprintf(output,"A reszmintak alapjan illesztett egyenes adatai:\n");
fprintf(output,"m=%f     a meredekseg\n",smk/N1);
fprintf(output,"Sm=%f     a hibaja\n",sqrt((N1*smk2-smk*smk)/(N1*N1*(N1-1))));
fprintf(output,"a=%f     a tengelymetszet\n",sak/N1);
fprintf(output,"Sa=%f     a hibaja\n",sqrt((N1*sak2-sak*sak)/(N1*N1*(N1-1))));
fprintf(output,"t1=%f     az illeszkedes kezdo ideje\n",stk/N1);
fprintf(output,"St1=%f     a hibaja\n",sqrt((N1*stk2-stk*stk)/(N1*N1*(N1-1))));
fprintf(output,"r2=%f     az illeszkedes josaga\n",srk/N1);
fprintf(output,"Sr2=%f     a hibaja\n\n",sqrt((N1*srk2-srk*srk)/(N1*N1*(N1-1))));

fprintf(output,"Az egyesitett minta alapjan illesztett egyenes adatai:\n");
fprintf(output,"m=%f     a meredekseg\n",mk);
fprintf(output,"a=%f     a tengelymetszet\n",ak);
fprintf(output,"t1=%f     az illeszkedes kezdo ideje\n",tk);
fprintf(output,"r2=%f     az illeszkedes josaga\n",rk);
fprintf(output,"hiba=%e     az illeszkedes hibaja\n\n",hiba_min);


fprintf(output,"Az egyesitett minta es az utolso reszminta reszletes adatai:\n\n");

fprintf(output,"k\t k*dt\t ");
fprintf(output,"sX_N[k]\t sX2_N[k]\t sV_N[k]\t sV2_N[k]\t ");
fprintf(output,"sX[k]\t sX2[k]\t sV[k]\t sV2[k]\t");
/* a kovetkezo sor tesztadatokhoz gyart fejlecet */
fprintf(output,"mkk[k]\t akk[k]\t rkk[k]\t hiba[k]\t");
/* a kovetkezo sor a konnyu abrazolas kedveert van */
fprintf(output,"\t k*dt\t D2X[k]\t egyenes\t");
fprintf(output,"\n");


for(k=0,kdt=0.0;k<=kmax;k++,kdt+=dt){
	fprintf(output,"%d\t",k);
	fprintf(output,"%f\t",kdt);
	fprintf(output,"%f\t",sX_N[k]);
	fprintf(output,"%f\t",sX2_N[k]);
	fprintf(output,"%f\t",sV_N[k]);
	fprintf(output,"%f\t",sV2_N[k]);
	fprintf(output,"%f\t",sX[k]);
	fprintf(output,"%f\t",sX2[k]);
	fprintf(output,"%f\t",sV[k]);
	fprintf(output,"%f\t",sV2[k]);
    /* a kovetkezo 5 sor teszt-adatokat ir ki */
	fprintf(output,"%f\t",mkk[k]);
	fprintf(output,"%f\t",akk[k]);
	fprintf(output,"%f\t",rkk[k]);
	fprintf(output,"%e\t",hibak[k]);
    /* a kovetkezo 3 sor a konnyu abrazolas kedveert van */
	fprintf(output,"\t%f\t",kdt);
	fprintf(output,"%f\t",sX2_N[k]/N2/N1);
	fprintf(output,"%f\t",mk*kdt+ak);
	fprintf(output,"\n");
	}
fflush(output);
fclose(output);

/* kiirunk egy sort a fo kimeneti fileba a bemeneti file minden sorahoz */
fprintf(output2,"%s\t",outfile);
fprintf(output2,"%d\t",teszt);
fprintf(output2,"%d\t",randmode);
fprintf(output2,"%lu\t",randseed);
fprintf(output2,"%d\t",N1);
fprintf(output2,"%d\t",N2);
fprintf(output2,"%f\t",TT);
fprintf(output2,"%f\t",dt);
fprintf(output2,"%d\t",uzemmod);
fprintf(output2,"%f\t",parameter);
fprintf(output2,"%f\t",mk);
fprintf(output2,"%f\t",ak);
fprintf(output2,"%f\t",tk);
fprintf(output2,"%f\t",rk);
fprintf(output2,"%f\t",smk/N1);
fprintf(output2,"%f\t",sqrt((N1*smk2-smk*smk)/(N1*N1*(N1-1))));
fprintf(output2,"%f\t",sak/N1);
fprintf(output2,"%f\t",sqrt((N1*sak2-sak*sak)/(N1*N1*(N1-1))));
fprintf(output2,"%f\t",stk/N1);
fprintf(output2,"%f\t",sqrt((N1*stk2-stk*stk)/(N1*N1*(N1-1))));
fprintf(output2,"%f\t",srk/N1);
fprintf(output2,"%f\t",sqrt((N1*srk2-srk*srk)/(N1*N1*(N1-1))));
fprintf(output2,"%f\t",difftime(ido_most,ido_volt));
fprintf(output2,"%s",ctime(&ido_volt));


} /* a bemeneti file sorait feldolgozo ciklus vege */


fclose(input);
fclose(output2);

}  /* main() vege */