/* A simple code in C to 
 
 1. draw the graph of a function (in this case the square root) and write the result on a file;
 2. compute the numerical derivative of the function in a point;
 3. compute numerically the integra of the function in an interval.

 */


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

#define MAX_DELTA 10000000

double Function(double);


void PlotGraph(double (*)(double),double,double,int);
double Derivative(double (*)(double),double,double);
double DerivativeSymmetric(double (*)(double),double,double);
double Integral(double (*)(double),double,double,int);

main(){
	double xmin,xmax;
	int NPoints,NDelta;
	double epsilon=1.;
	double x,der,dersymm,DerAnalitic,integral,IntAnalitic;
	double (* Ff)();
	
	Ff=&Function;
	
	
	printf("Please insert the lower bound of the plot range\n");
	scanf("%lf",&xmin);
	printf("Please insert the upper bound of the plot range\n");
	scanf("%lf",&xmax);
	
	do{
		printf("Insert the number of points in the graph\n");
		scanf("%d",&NPoints);
	}while(NPoints<=0);
	
	PlotGraph(Function,xmin,xmax,NPoints); 
	
	
	
	
	printf("Please insert the value at which the derivative is to be evaluated \n");
	
	do{
		scanf("%lf",&x);
		if(x==0.)printf("The derivative is ill defined at this point, retry\n");
	}while(x==0.);
	
		DerAnalitic=0.5/sqrt(x);
		
	printf("epsilon       Num. der.   Exact result at x=%f. f(x) f(x+epsilon) \n\n",x);

	while (epsilon>1.e-18) {
		der=Derivative(Function,x,epsilon); 
		dersymm=DerivativeSymmetric(Function,x,epsilon); 
	
		printf("%e \t %f \t %f \t %f \t %.12f \t %.12f \n",epsilon,der,dersymm,DerAnalitic,sqrt(x),sqrt(x+epsilon));
		epsilon/=10.;
	}

	
	printf("\n\n ######################################################### \n\n");
	
	printf("Please insert lower extreme of the integral\n");
	scanf("%lf",&xmin);
	
	printf("Please insert upper extrem of the integral\n");
	scanf("%lf",&xmax);
	
	printf("Please insert initial number of intervals \n");
	scanf("%d",&NDelta);
	
	
	IntAnalitic=2./3.*(pow(xmax,1.5)-pow(xmin,1.5));
	
	printf("NDelta      Num. Int.   Exact result \n\n");

	
	while(NDelta<MAX_DELTA) {
		integral=Integral(Function,xmin,xmax,NDelta);
		
		printf("%.7d  %f  %f \n",NDelta,integral,IntAnalitic);
		NDelta =2.*NDelta;
	}
	

}


double Function(double x){
	
	return sqrt(x);
	
}


void PlotGraph(double (*f)(double),double min, double max,int NP){
	int i;
	double x,Delta,fun;
	
	FILE *Plot;
	
	Plot=fopen("Plot.dat","w");
	
	Delta=(max-min)/(double)NP;
	
	for (i=0; i<=NP; i++) {
		x=min+i*Delta;
		fun=f(x);
		fprintf(Plot,"%lf %lf \n",x,fun);
	}
	
	fclose(Plot);
}



double Derivative(double (*f)(double),double x,double epsilon){ 

	double d;
	
	d = (f(x+epsilon)-f(x))/epsilon;
	
	return d;
	
}

double DerivativeSymmetric(double (*f)(double),double x,double epsilon){ 
	
	double d;
	
	d = (f(x+epsilon)-f(x-epsilon))/(2.*epsilon);
	
	return d;
	
}

double Integral(double (*f)(double),double min, double max,int NP){
	int i;
	double x, Delta, integ=0.;
	
	Delta=(max-min)/NP;

	for(i=0;i<NP;i++){
		x=min+(i+0.5)*Delta;
		integ+=f(x);
	}
	integ=integ*Delta;
	
	return integ;
}