#ifndef __RUNGEKUTTA_H__
#define __RUNGEKUTTA_H__


/* Ho usato questo typedef per poter facilmente trasformare questi 
   integrtatori da double a float
 */
typedef double rkfloat;


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

/*
 * Queste sono due implementazioni personali delle funzioni rkqs ed rkck 
 * descritte nella sezione 16.2 di 
 * "NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING",
 * www.nr.com
 *
 * rkck integra con un metodo di Runge-Kutta una funzione dal 
 * punto di partenza (x,y[]) al punto (x+h,ynext[]), dove gli y sono
 * dei vettori di dimensione n. yerr[] contiene per ogni variabile del
 * problema una stima dell'errore di integrazione.
 * function è un puntatore ad una funzione che, dato il punto (x,y[]) 
 * calcola le derivate dy[]/dx.
 *
 * rkqs è la funzione che, dato il punto inizial (x,y[]) calcola qual'è il
 * passo migliore per per effettuare un'integrazione con la precisione 
 * richiesta. All'uscita da rkqs, x risulta incrementato della quantità
 * stepdone.
 * La precisione richiesta è definita dal'argomento precision che rappresenta
 * la precisione assoluta se y>1. o relativa se y<1.
 */

void rkck( rkfloat x, int n, rkfloat y[], 
           rkfloat h,        rkfloat ynext[], rkfloat yerr[],
           void (*function)(rkfloat, rkfloat*, rkfloat*) );

void rkqs( rkfloat *x, int n, rkfloat y[],
           rkfloat errscal[], rkfloat precision,
           rkfloat steptry, rkfloat *stepdone, rkfloat *stepnext,
           void (*function)(rkfloat, rkfloat*, rkfloat*) );

/*
 * Per le nostre applicazioni lo step variable di rkqs potrebbe non essere 
 * adatto al passo di campionamente che vogliamo per fare il grafico di una 
 * funzione. rkad effettua l'integrazione su un intervallo fissato a priori,
 * da (xstart,ystart[]), fino a (xend,yend[]), nell'intervallo (xstart,xend)
 * il passo di integrazione viene aggiustato usanto rkqs per mantenere 
 * l'errore entro il limite fissato da precision. 
 */

void rkad( rkfloat xstart, rkfloat xend,
	   int n, rkfloat ystart[], rkfloat yend[],
	   rkfloat precision,
           void (*function)(rkfloat, rkfloat*, rkfloat*) );


#endif
