#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <float.h>
#include <conio.h>
#define MAX 17

int spline(int n, double *x, double *y, double marg_0, double marg_n, 
		   int marg_cond, double *b, double *c, double *d);
double spval (int n, double x0, double *a, double *b, double *c, double *d,
			  double *x, double *ausg);

int main (void)
{
//========================Funktionswerte ausgeben===========================
	int i,n=16;
	double	x[MAX],f[MAX];
	double	h,sx,fx,x0,diff;
	double b[MAX], c[MAX], d[MAX],abl[3], dmax[3];

	FILE *Dat;
		
	Dat=fopen("Funktionswerte.sav","w");

	printf("\n Hoehere Mathematik II Praktikum "
		   "\n 3.Aufgabe (Spline-Interpolation)"
		   "\n Marcel Cuvelier - Gruppe MAMF"
		   "\n ================================\n\n ");

	printf("\n Zur Berechnung der Werte... ");
	system("pause");

	system("cls");	
	printf("\n Funktionswerte:\n");
	printf("\n =================");
	printf("\n\n  i\t    x[i]\t    f[i]");
	printf("\n ---------------------------------------");

	h = 2./n;

	for (i=0; i<=n; i++)
	{	
		x[i] = -1 + i*h;
		f[i] = exp(x[i]); 
		printf("\n %i\t %lf\t %e",i, x[i], f[i]);
		fprintf(Dat,"\n %i\t %lf\t %e",i, x[i], f[i]);					
	}

	printf("\n\n ");
	fclose(Dat);
	system("pause");


//==================Aufgabe a) Koeffizienten von s(x)======================
	
	system("cls");				
	Dat=fopen("aufgabe3a.sav","w");
	
	h = 2./n;

	for (i=0; i<=n; i++)
	{	
		x[i] = -1 + i*h;
		f[i] = exp(x[i]);
	}

	spline(n+1, x, f, f[0], f[n], 1,b,c,d);
	
	printf("\n\n Aufgabe a) Koeffizienten von s(x):\n");
	printf("\n ====================================");
	printf("\n\n  i\t    a[i]\t    b[i]\t     c[i]\t      d[i]");
	printf("\n ----------------------------------------------------------------------");

	for (i=0; i<n; i++)
	{
		printf("\n  %i\t %e\t %e\t %e\t %e", i, f[i], b[i], c[i], d[i]);   // a[i] = f[i]
		fprintf(Dat,"\n %i\t %e\t %e\t %e\t %e", i, f[i], b[i], c[i], d[i]);
	}

	printf("\n\n ");
	system("pause");
	fclose(Dat);

//==================Aufgabe b) Plotten der Funktionen======================
	
	system("cls");	
	Dat=fopen("aufgabe3b.sav", "w");

	h=2./1024.;
	for(i=0;i<1024;i++)
	{
		x0=-1+i*2./1024.;
		fx=exp(x0);
		sx=spval(n+1,x0,f,b,c,d,x,abl);
		fprintf (Dat,"\n%4d %9.5lf %9.5e %9.5e ",i,x0,fx,sx);
		fprintf (Dat,"%9.5e %9.5e %9.5e"	,abl[0],abl[1],abl[2]);
		fprintf (Dat," %9.5e %9.5e"	,fx-sx,fx-abl[0]);
	}

	printf("\n\n Diese Werte wurden in der Datei -aufgabe3b.sav- abgespeichert!");
	
	fclose(Dat);
	printf("\n\n ");
	system("pause");

//==================Aufgabe c) Äquidistante Stützstellen====================

	system("cls");	
	Dat=fopen("aufgabe3c.sav", "w");

	dmax[0]=0;
	dmax[1]=0;
	dmax[2]=0;

	for(i=0;i<16384;i++)
	{
		x0=-1.+i*2./16384.;
		fx=exp(x0);
		sx=spval(n+1,x0,f,b,c,d,x,abl);
		diff=fabs(fx-sx);
		if(dmax[0]<diff)dmax[0]=diff;
		diff=fabs(fx-abl[0]);
		if(dmax[1]<diff)dmax[1]=diff;
		diff=fabs(fx-abl[1]);
		if(dmax[2]<diff)dmax[2]=diff;
	}

	printf("\n Aufgabe 3c:\n\n");
	
	printf("\n (i):   max  | f(x)-s(x)  |     =%15.5e", dmax[0]);
	printf("\n (ii):  max  | f'(x)-s'(x)  |   =%15.5e", dmax[1]);
	printf("\n (iii): max  | f''(x)-s''(x)  | =%15.5e", dmax[2]);
	
	fprintf(Dat, "\n (i):   max|f(x)-s(x)|     =%15.5e", dmax[0]);
	fprintf(Dat, "\n (ii):  max|f''(x)-s'(x)|  =%15.5e", dmax[1]);
	fprintf(Dat, "\n (iii): max|f''(x)-s''(x)| =%15.5e", dmax[2]);

	fclose(Dat);
	printf("\n\n ");
	printf("Diese Werte wurden in der Datei -aufgabe3c.sav- abgespeichert!\n ");
	system("pause");

//==============Aufgabe d) Kubische Spline Funktion ? ==============

	system("cls");	
	Dat=fopen("aufgabe3d.sav", "w");
	for(i=0;i<512;i++)
	{
		x0=(-1.+i*4./1024.);
		fx=exp(-x0);
		printf("\n %4d %9.5lf %9.5e ", i, x0, fx);
		fprintf(Dat, "\n%4d %9.5lf %9.5e ", i, x0, fx);
	}
	printf("\n ");
	system("pause");
	for(i=512;i<=1024;i++)
	{
		x0=(-1.+i*4./1024.);
		fx=exp(x0-2);

		printf("\n %4d %9.5lf %9.5e ", i, x0, fx);
		fprintf(Dat, "\n%4d %9.5lf %9.5e ", i, x0, fx);
	}
	fclose(Dat);
	printf("\n\n Diese Werte wurden in der Datei -aufgabe3d.sav- abgespeichert!");
	printf("\n Das Programm wird nun beendet !!\n\n ");
	return 0;
}//Main

//Vorgegebene Funktionen---------------------------------------------------------
/* ----------------------------- Modul Tridiagonal --------------------------- */
int trdiag (int n, double *lower, double *diag, double *upper, double *b, int rep)
{
	register i;

	if(n < 2) return 1;
	if(rep == 0)
	{
		for(i=1; i<n; i++)
		{
			if(fabs(diag[i-1]) < DBL_EPSILON) return 2;
			lower[i] /= diag[i-1];
			diag[i] -= lower[i] * upper[i-1];
		}
	}
	if(fabs(diag[n-1]) < DBL_EPSILON) return 2;
	for(i=1; i<n; i++)
		b[i] -= lower[i] * b[i-1];
	b[n-1] /= diag[n-1];
	for(i=n-2; i>=0; i--)
		b[i] = (b[i] - upper[i] * b[i+1]) / diag[i];
	return 0;
}
/* ----------------------------- Ende Tridiagonal ---------------------------- */

/* ----------------------------- Modul Spline -------------------------------- */

int spline(int n, double *x, double *y, double marg_0, double marg_n, 
		   int marg_cond, double *b, double *c, double *d)
/********************************************************************************
* spline berechnet zu den vorgegebenen Wertepaaren                              *
*      (x[i], y[i], i=0(1)n-1)                                                  *
* die Koeffizienten eines kubischen Splines.                                    *
* Die Art der Randbedingung wird durch die Variable marg_cond                   *
* festgelegt.                                                                   *
*===============================================================================*
* Eingabeparameter:                                                             *
* -----------------                                                             *
* Name		Typ/Länge	Bedeutung                                               *
* ------------------------------------------------------------------------------*
* n			int/---		Anzahl der Stützstellen; n>2                            *
* x			double/[n]	x-Koordinaten der Wertepaare                            *
* y			double/[n]	y-Koordinaten der Wertepaare                            *
* marg_0	double/---	Randbedingung in x[0]                                   *
* marg_n	double/---	Randbedingung in x[n-1]                                 *
* marg_cond int/---		= 0 : not-a-knot-Bedingung                              *
*							  marg_0, marg_n ohne Bedeutung                     *
*						= 1 : marg_0, marg_n sind 1. Ableitungen                *
*						= 2 : marg_0, marg_n sind 2. Ableitungen                *
*                       = 3 : marg_0, marg_n sind 3. Ableitungen                *
*                                                                               *
* Ausgabeparameter:                                                             *
* -----------------                                                             *
* Name		Typ/Länge	Bedeutung                                               *
* ------------------------------------------------------------------------------*
* b			double/[n]	Splinekoeffizienten:                                    *
* c			double/[n]	s(x)=a[i]+b[i]*t+c[i]*t^2+d[i]*t^3                      *
* d         double/[n]	t=x-x[i]; x aus [ x[i], x[i+1] ]                        *
*						(a entspricht y)                                        *
*                       Die Splinekoeffizienten brauchten nur mit [n-1]         *
*						vereinbart werden. Es wird jedoch bei der Rechnung      *
*						[n] als Zwischenspeicher benötigt.                      *
*                                                                               *
* Rückgabewert:                                                                 *
* -------------                                                                 *
* = -i : Monotoniefehler: x[i-1] >= x[i]                                        *
* =  0 : kein Fehler                                                            *
* =  1 : falscher Wert für marg_cond                                            *
* =  2 : n < 3                                                                  *
* =  3 : es steht nicht genügend Speicher für die notwendigen Hilfsfelder       *
*        zur Verfügung                                                          *
* >  3 : Fehler in trdiag                                                       *
********************************************************************************/
{
	int i, error;
	double *h, *a;

	if(n < 3) return (2);
	/* es wird Speicherplatz für das Hilfsfeld zur Verfügung gestellt */
	a = (double *) malloc (n * sizeof (double));
	if(a == NULL) return 3;
	h = (double *) malloc ((--n) * sizeof (double));
	if(h == NULL) { free(a); return 3; }

	for(i=0; i<=n-1; i++)
	{
		h[i] = x[i+1] - x[i];
		if(h[i] <= 0.) return (-i);  /* Ueberprüfung der Monotonie */
	}
	/* Aufstellen des Gleichungssystems */
	for(i=0; i<=n-2; i++)
	{
		a[i] = 3. * ((y[i+2] - y[i+1]) / h[i+1]
			- (y[i+1] - y[i]) / h[i]);
		b[i] = h[i];
		c[i] = h[i+1];
		d[i] = 2. * (h[i] + h[i+1]);
	}
	/* In Abhängigkeit von der Randbedingung werden Werte neu bestezt */
	switch (marg_cond)
	{
		case 0: 
			{
				if(n == 2)
				{
					a[0] /= 3.;
					d[0] *= 0.5;
				}
				else
				{
					a[0] *= h[1] / (h[0] + h[1]);
					a[n-2] *= h[n-2] / (h[n-1] + h[n-2]);
					d[0] -= h[0];
					d[n-2] -= h[n-1];
					c[0] -= h[0];
					b[n-2] -= h[n-1];
				}
				break;
			}
		case 1: 
			{
				a[0] -= 1.5 * ((y[1] - y[0]) / h[0] - marg_0);
				a[n-2] -= 1.5 * (marg_n - (y[n] - y[n-1]) / h[n-1]);
				d[0] -= h[0] * 0.5;
				d[n-2] -= h[n-1] * 0.5;
				break;
			}
		case 2:
			{
				a[0] -= h[0] * marg_0 * 0.5;
				a[n-2] -= h[n-1] * marg_n * 0.5;
				break;
			}
		case 3:
			{
				a[0] += marg_0 * h[0] * h[0] * 0.5;
				a[n-2] -= marg_n * h[n-1] * h[n-1] * 0.5;
				d[0] += h[0];
				d[n-2] += h[n-1];
				break;
			}
		default : {free(a); free(h); return 1;}
	}
	/* Berechnen der Koeffizienten */
	switch (n-1)
	{
		case 1:
			{
				c[1] = a[0] / d[0];
				break;
			}
		default :
			{
				error = trdiag(n-1, b, d, c, a, 0);
				if(error != 0) { free(a); free(h); return error+3; }

				for(i=0; i<=n-2; i++)  /* Ueberschreiben des Lösungs- */
					c[i+1] = a[i];     /* vektors auf c               */
			}
	}
	/* In Abhängigkeit von der Randbedingung wird der erste und der letzte 
	   Wert von c neu besetzt. */
	switch (marg_cond)
	{
		case 0:
			{
				if(n == 2)
					c[0] = c[2] = c[1];
				else
				{
					c[0] = c[1] + h[0] * (c[1] - c[2]) / h[1];
					c[n] = c[n-1] + h[n-1] * 
						(c[n-1] - c[n-2]) / h[n-2];
				}
				break;
			}
		case 1:
			{
				c[0] = 1.5 * ((y[1] - y[0]) / h[0] - marg_0);
				c[0] = (c[0] - c[1] * h[0] * 0.5) / h[0];
				c[n] = -1.5 * ((y[n] - y[n-1]) / h[n-1] - marg_n);
				c[n] = (c[n] - c[n-1] * h[n-1] * 0.5) / h[n-1];
				break;
			}
		case 2:
			{
				c[0] = marg_0 * 0.5;
				c[n] = marg_n * 0.5;
				break;
			}
		case 3:
			{
				c[0] = c[1] - marg_0 * h[0] * 0.5;
				c[n] = c[n-1] + marg_n * h[n-1] * 0.5;
			}
	}
	for(i=0; i<=n-1; i++)
	{
		b[i] = (y[i+1] - y[i]) / h[i] - h[i] *
			(c[i+1] + 2. * c[i]) / 3.;
		d[i] = (c[i+1] - c[i]) / (3. * h[i]);
	}
	free(h);
	free(a);
	return 0;
}
/* ----------------------------- Ende Spline --------------------------------- */

/* ----------------------------- Modul SPVAL --------------------------------- */
double spval (int n, double x0, double *a, double *b, double *c, double *d,
			  double *x, double *ausg)
/********************************************************************************
* spval berechnet Funktionswerte eines kubischen Splines und                    *
* seine nicht-trivialen Ableitungen.                                            *
*===============================================================================*
* Eingabeparameter:                                                             *
* -----------------                                                             *
* Name		Typ/Länge		Bedeutung                                           *
* ------------------------------------------------------------------------------*
* n			int/---			Anzahl der Stützstellen                             *
* x0		double/---		Stelle, an der der Funktionswert berechnet wird     *
* a			double/[n-1]	)                                                   *
* b			double/[n-1]	) Splinekoeffizienten                               *
* c			double/[n-1]	)                                                   *
* d			double/[n-1]	)                                                   *
* x			double/[n-1]	Stützstellen                                        *
*                                                                               *
* Ausgabeparameter:                                                             *
* -----------------                                                             *
* Name		Typ/Länge		Bedeutung                                           *
* ------------------------------------------------------------------------------*
* ausg		double/[3]		Ableitungen an der Stelle x0                        *
*							ausg[0] enthält die 1. Ableitung                    *
*							ausg[1] enthält die 2. Ableitung                    *
*							ausg[2] enthält die 3. Ableitung                    *
*                                                                               *
* Rückgabewert:                                                                 *
* -------------                                                                 *
* Funktionswert der Splinefunktion an der Stelle                                *
*               x0,     falls x[0]<=x0<=x[n-1]                                  *
*               x[0],   falls x0<x[0]                                           *
*               x[n-1], falls x[n-1]<x0                                         *
********************************************************************************/
{
	int i, k, m;

	if(x0 < x[0]) x0 = x[0];
	if(x0 > x[n-1]) x0 = x[n-1];

	i = 0, k = n;
	while(m = (i+k) >> 1, m!= i)	/* das Intervall, in dem x0 */
	{								/* liegt, wird gesucht.     */
		if(x0 < x[m]) k = m;
		else i = m;
	}
	if (i >= n-1) i=n-2;
	x0 -= x[i];
	ausg[0] = (3. * d[i] * x0 + 2. * c[i]) * x0 + b[i];
	ausg[1] = 6. * d[i] * x0 + 2. * c[i];
	ausg[2] = 6. * d[i];
	return (((d[i] * x0 + c[i]) * x0 + b[i]) * x0 + a[i]);
}

/* ----------------------------- Ende SPVAL ---------------------------------- */