/* Questo programma svolge (mediante apposite funzioni) diverse "operazioni" utili 
 sulle funzioni:
 
 1. costruisce il grafico della funzione in un certo intervallo e lo stampa su file;
 2. calcola numericamente la derivata della funzione in un punto;
 3. calcola numericamente l'integrale definito della funzione in un intervallo.
 
 La funzione scelta in questo caso e' la radice quadrata. Cambiarla comporta, 
 essenzialmente, modificare il solo contenuto della funzione Function.
 
 Il punto "tecnico" nuovo sono i puntatori a funzione */


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

#define MAX_DELTA 10000000

double Function(double);

/* Notare il primo argomento nelle funzioni qui sotto: 
 e' un puntatore ad una	funzione -- per ora senza nome -- 
 che restituisce un double ed ha un argomento che e' un double */

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

main(){
	double xmin,xmax;
	int NPoints,NDelta;
	double epsilon=1.;
	double x,der,DerAnalitic,integral,IntAnalitic;
	double (* Ff)();
	
	Ff=&Function;
	
	printf("Ff=%p   Function=%p\n\n",Ff, Function);
	
	printf("Inserire l'estremo inferiore dell'intervallo in cui si vuol disegnare la funzione\n");
	scanf("%lf",&xmin);
	printf("Inserire l'estremo superiore dell'intervallo in cui si vuol disegnare la funzione\n");
	scanf("%lf",&xmax);
	
	do{
		printf("Inserire il numero di punti (>0) nell'intevallo in cui si vuol disegnare la funzione\n");
		scanf("%d",&NPoints);
	}while(NPoints<=0);
	
	PlotGraph(Function,xmin,xmax,NPoints); /* la "funzione" viene passata con il primo argomento nella 
										chiamata, passando in realta'il puntatore alla prima locazione 
										di memoria che essa occupa */
	
	
	
	
	printf("Inserire il valore del punto in cui calcolare la derivata \n");
	
	do{
		scanf("%lf",&x);
		if(x==0.)printf("La derivata e' singolare nel punto scelto, ripetere la scelta\n");
	}while(x==0.);
	
		DerAnalitic=0.5/sqrt(x);
		
	printf("epsilon       Der. num.   risultato esatto per x=%f.\n\n",x);

	while (epsilon>1.e-18) {
		der=Derivative(Function,x,epsilon); /* la "funzione" viene passata con il primo argomento nella 
											chiamata, passando in realta'il puntatore alla prima locazione 
											di memoria che essa occupa */
	
		printf("%e  %f  %f \n",epsilon,der,DerAnalitic);
		epsilon/=10.;
	}

	
	printf("\n\n ######################################################### \n\n");
	
	printf("Inserire l'estremo inferiore dell'integrale\n");
	scanf("%lf",&xmin);
	
	printf("Inserire l'estremo superiore dell'integrale\n");
	scanf("%lf",&xmax);
	
	printf("Inserire il numero di intervalli di partenza \n");
	scanf("%d",&NDelta);
	
	
	IntAnalitic=2./3.*(pow(xmax,1.5)-pow(xmin,1.5));
	
	printf("NDelta      Int. num.   risultato esatto \n\n");

	
	while(NDelta<MAX_DELTA) {
		integral=Integral(Function,xmin,xmax,NDelta); /* la "funzione" viene passata con il primo argomento nella 
													chiamata, passando in realta'il puntatore alla prima locazione 
													di memoria che essa occupa */
		
		printf("%.7d  %f  %f \n",NDelta,integral,IntAnalitic);
		NDelta =2.*NDelta;
	}
	

}


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

/*Il nome della funzione all'interno di PlotGraph
 -- f -- e'fissato nell'argomento */
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 Integral(double (*f)(double),double min, double max,int NP){
	int i;
	double x, Delta, integ=0.;
	
	Delta=(max-min)/NP;

	/* Il ciclo qui sotto, quando vi ho fatto vedere il programma in classe conteneva
	 un errore. Quale? Osservate l'effetto della correzione sulla precisione
	 dell'integrale numerico e provate a capire l'origine della differenza. */
	for(i=0;i<NP;i++){
		x=min+(i+0.5)*Delta;
		integ+=f(x);
	}
	integ=integ*Delta;
	
	return integ;
}