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

#define NPART_MAX 10000

double YDisplacement(double, double);
double MeanValue(double YData[], int);
double StdDev(double YData[], double, int);

int NInStdDev(double YData[], double, int);

main(){

	int i,NStep,seed;
	int n,Np,NInS;
	double dx,dh,L,test;
	double YFinal[NPART_MAX], OldY,NewY;
	double MeanEndY,StdDevEndY;
	
	printf("Inserire la profondita' del mezzo, L\n");
	scanf("%lf",&L);
	do{
	printf("\n Inserire il passo, dx\n");
	scanf("%lf",&dx);
		test=1.;
		if(dx!=0) test=L/dx;
		test=test/(int)test;
		if(test!=1.)printf("L deve essere divisibile per dx. Inserire un nuovo valore per dx");
	}while(dx==0.||test!=1.);
	
	NStep=L/dx;
	
	printf("Inserire il modulo della deviazione, dh\n");
	scanf("%lf",&dh);
	
	printf("\nInserire il numero di particelle, Np (<100)\n");
	do{
		scanf("%d",&Np);
		if(Np>NPART_MAX)printf("Numero di particelle troppo grande!\n");
	}while(Np>NPART_MAX);
	
	
	
	for(n=1;n<=Np;n++){
		OldY=0.;
		for(i=0;i<NStep;i++){
			NewY=YDisplacement(dh,OldY);
	//		printf("NewY %lf\n",NewY);
			OldY=NewY;
	}
		YFinal[n-1]=NewY;
	}
	
/*	for(n=0;n<Np;n++){
		printf("Particella %d, Posizione %lf \n",n,YFinal[n]);
	} */
	
	MeanEndY=MeanValue(YFinal,Np);
	
	StdDevEndY=StdDev(YFinal,MeanEndY,Np);
	
	printf("Il valore medio della coordinata y alla fine della traiettoria e' %lf\n",MeanEndY);
	
	printf("La deviazione standard dei valori alla fine della traiettoria e' %lf\n",StdDevEndY);
	
	
	NInS=NInStdDev(YFinal,StdDevEndY,Np);
	
	printf("Il numero di particelle che cadono in [-s,s] e' %d\n",NInS);
	


}

double YDisplacement(double yshift, double Ypos){

	double p, NewYpos; 
	
	p=(double)lrand48()/RAND_MAX;
	
	if(p<0.5){
		NewYpos=Ypos-yshift;
	}else{
		NewYpos=Ypos+yshift;
	}

	return NewYpos;

}

double MeanValue(double Data[],int NData){
	int i;
	double mu=0.;
	
	for(i=0;i<NData;i++){
		mu+=*(Data+i);
	}
	
	mu=mu/NData;
	
	return mu;
}

double StdDev(double Data[],double mu,int NData){
	int i;
	double s=0.;
	
	for(i=0;i<NData;i++){
		s+=(*(Data+i)-mu)*(*(Data+i)-mu);
	}
	
	if(NData!=1){
		s=s/(NData-1.);
	}else{
		printf("Numero di Dati insufficiente a calcolare la deviazione standard");
	}
	
	return sqrt(s);

}

int NInStdDev(double Data[], double s, int NData){
	int n, Ntot=0;
	double ease;
	
	for(n=0;n<NData;n++){
		ease=*(Data+n);
		if(fabs(ease)<=s) Ntot+=1;
	}
	
	return Ntot;
}