/*********************** Compute the invariant density for Logistic Map ************************/



//Include libraries

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



main()
{

	int i, j, n, NOB, count;

	double x_n, x_n1, mu, x_min, x_max, dx, x[5005], y;

	FILE*fa = fopen( "analytical_invariant_density.dat", "w" );

	FILE*fw = fopen( "invariant_density.dat", "w" );

	printf( "Enter the value of mu: \t");

	scanf( "%lf", &mu );

	printf( "Enter the value of n: \t");

	scanf( "%d", &n );

	printf( "Enter the value of x0 between 0 and 1: \t");

	scanf( "%lf", &x_n );

	printf( "Enter the number of boxes for histogram: \t");

	scanf( "%d", &NOB );

//~~~~~~~~~~~~~~~~~~~~~~~~~~Logistic Map~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


	for( i = 1; i <= n; i++ )
	{

		x_n1 = mu * x_n * ( 1 - x_n );

		x_n = x_n1;

		x[i] = x_n;

	}


//~~~~~~~~~~~~~~~~~Initiate x_min and x_max~~~~~~~~~~~~~~~~~~~~~~~~~~~~


	x_min = 0.5;

	x_max = 0.5;


//~~~~~~~~~~~~~~~~~Find minimum and maximum value of x~~~~~~~~~~~~~~~~~

	for( i = 1; i <= n; i++ )
	{

		if( x[i] <= x_min )	x_min = x[i];

		if( x[i] >= x_max )	x_max = x[i];

	}

	dx = ( x_max - x_min ) / ( NOB * 1.0 );


//~~~~~~~~~~~~~~~~~~~~~~~~~~~Histogram~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	
	for( i = 1; i <= NOB; i++ )
	{

		count = 0;

		for( j = 1; j <= n; j++ )
		{

			if( x[j] > ( i - 1 ) * dx && x[j] <= i * dx )	count += 1;

		}

		y = ( i - 0.5 ) * dx;

		fprintf( fw, "%0.6lf\t%lf\n", y , NOB * count / ( n * 1.0 ) );

		fprintf( fa, "%0.6lf\t%0.6lf\n", y, 1.0 / ( M_PI * sqrt ( y * ( 1 - y ) ) ) );


	}



	fclose(fa);

	fclose(fw);


}