#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdlib.h>
#include "nrutil.h"
#include "nr.h"

//#define N 20
//#define a 0
//#define b 1

int N=20;
double a=0.0, b=1.0;
double h;
float *wezly;
float *U;
float *diag_diag;
float *lower_diag;
float *upper_diag;
float *vec_S;
float P_1(float t);
float Q_1(float t);
float F_1(float t);
float U_1(float t);


//int i,j,k,l,m,n,h=1.0/(N+1.0);
//h=1/(N+1);
float C1(int j,int i);
float C1_falka(int j,int i);
float C2(int j,int m,int i);
float C2_falka(int j,int m,int i);


float P_1(float t){return t*t;}
float Q_1(float t){return t*t;}
float U_1(float t){return t*t;}
float F_1(float t){return t*t;}

/* run this program using the console pauser or add your own getch, system("pause") or input loop */

void metoda_galerkina(float(*U_2)(float),float(*P)(float),float(*Q)(float),float(*F)(float)){
int k;

wezly=vector(0,N+1);
U = vector(0,N+1);


	for(k=1; k<=N;wezly[k]=a+k*(b-a)/(N+1.0),k++);



float C1(j,i){
	if(i==j){return 2*h/3.0;}
	if(j-i > 1 || j-i < -1){return 0;}
	else{return h/6.0;}
}

float C1_falka(j,i){
	if(i==j){return 2.0/h;}
	if(j-i > 1 || j-i < -1){return 0;}
	else{return (-1.0)/h;}
}

float C2(j,m,i){
	if(i==j && j==m){return h/2.0;}
	if(abs(j-i)>1 || abs(j-m)>1 || abs(i-m)>1){return 0;}
	else{return h/12.0;}
}

float C2_falka(j,m,i){
	if(i==j && j==m){return 0;}
	if(abs(j-i)>1 || abs(j-m)>1 || abs(i-m)>1){return 0;}
	if(j==m && m==i-1){return 1.0/3;}
	if(j==m && m==i+1){return (-1)/3.0;}
	if(j==i && m-1==i){return (-1)/6.0;}
	if(j==i && m+1==i){return 1.0/6;}
	if(j==m+1 && m+1==i+1){return (-1)/6.0;}
	if(j==m-1 && m-1==i-1){return 1.0/6;}
}


float *diag_diag = vector(1,N);
float *lower_diag = vector(1,N);
float *upper_diag = vector(1,N);
float *vec_S = vector(1,N);

lower_diag[1] = 0;

upper_diag[N] = 0;


diag_diag[1] = (-1)*C1_falka(1,1) + P_1(wezly[0])*C2_falka(1,0,1) + Q_1(wezly[0])*C2(1,0,1) + P_1(wezly[1])*C2_falka(1,1,1) + Q_1(wezly[1])*C2(1,1,1) +  P_1(wezly[2])*C2_falka(1,2,1) + Q_1(wezly[2])*C2(1,2,1);
upper_diag[1] = (-1)*C1_falka(1,2) + P_1(wezly[1])*C2_falka(1,1,2) + Q_1(wezly[1])*C2(1,1,2) + P_1(wezly[2])*C2_falka(1,2,2) + Q_1(wezly[2])*C2(1,2,2);
lower_diag[N] = (-1)*C1_falka(N,N-1) +  P_1(wezly[N-1])*C2_falka(N,N-1,N-1) + Q_1(wezly[N-1])*C2(N,N-1,N-1) + P_1(wezly[N])*C2_falka(N,N,N-1) + Q_1(wezly[N])*C2(N,N,N-1);
diag_diag[N] = (-1)*C1_falka(N,N) + P_1(wezly[N-1])*C2_falka(N,N-1,N) + Q_1(wezly[N-1])*C2(N,N-1,N) + P_1(wezly[N])*C2_falka(N,N,N) + Q_1(wezly[N])*C2(N,N,N) +  P_1(wezly[N+1])*C2_falka(N,N+1,N) + Q_1(wezly[N+1])*C2(N,N+1,N);

for(k=2;k=N-1;k++){
lower_diag[k] = (-1)*C1_falka(k,k-1) +  P_1(wezly[k-1])*C2_falka(k,k-1,k-1) + Q_1(wezly[k-1])*C2(k,k-1,k-1) + P_1(wezly[k])*C2_falka(k,k,k-1) + Q_1(wezly[k])*C2(k,k,k-1);
diag_diag[k] = (-1)*C1_falka(k,k) + P_1(wezly[k-1])*C2_falka(k,k-1,k) + Q_1(wezly[k-1])*C2(k,k-1,k) + P_1(wezly[k])*C2_falka(k,k,k) + Q_1(wezly[k])*C2(k,k,k) +  P_1(wezly[k+1])*C2_falka(k,k+1,k) + Q_1(wezly[k+1])*C2(k,k+1,k);
upper_diag[k] = (-1)*C1_falka(k,k+1) + P_1(wezly[k])*C2_falka(k,k,k+1) + Q_1(wezly[k])*C2(k,k,k+1) + P_1(wezly[k+1])*C2_falka(k,k+1,k+1) + Q_1(wezly[k+1])*C2(k,k+1,k+1);
}


vec_S = vector(1,N);
vec_S[1] = F_1(wezly[0])*C1(1,0) + F_1(wezly[1])*C1(1,1) + F_1(wezly[2])*C1(1,2) + U_1(a)*(C1_falka(1,0) - P_1(wezly[0])*C2_falka(1,0,0) - P_1(wezly[1])*C2_falka(1,1,0) - Q_1(wezly[0])*C2(1,0,0) - Q_1(wezly[1])*C2(1,1,0));
vec_S[N] = F_1(wezly[N-1])*C1(N,N-1) + F_1(wezly[N])*C1(N,N) + F_1(wezly[N+1])*C1(N,N+1) + U_1(b)*(C1_falka(N,N+1) - P_1(wezly[N])*C2_falka(N,N,N+1) - P_1(wezly[N+1])*C2_falka(N,N+1,N+1) - Q_1(wezly[N])*C2(N,N,N+1) - Q_1(wezly[N+1])*C2(N,N+1,N+1));

int j;
for(j=2;j=N-1;j++){
	vec_S[j] = F_1(wezly[j-1])*C1(j,j-1) + F_1(wezly[j])*C1(j,j) + F_1(wezly[j+1])*C1(j,j+1);
}


tridag(lower_diag, diag_diag, upper_diag, vec_S, U, N);
}

int main() {

	float *results = vector(1,N);
	float *check = vector(1,N);

	int f;
	metoda_galerkina(U_1,P_1,Q_1,F_1);

	for(f=1;f=N;f++){
	check[f] = U_1(wezly[f]) - U[f];}

	return 0;

}




//    // wektor rozwiazan
//        U=dvector(0,N+1);
//
//    //  Alg Thomasa
//        beta=dvector(1,N);
//        gamma=dvector(1,N);
//
//    //   Wyliczane wartosci funkcji w punktach
//
//    // Wektor wyrazow wolnych
//        D[1]=h*h*f[1]-(1-p[1]*h*0.5)*x_a;
//        D[N]=h*h*f[N+1]-(1-p[N]*h*0.5)*x_b;
//        for(int i=2; i<N; i++)
//        {
//            D[i]=h*h*f[i];
//        }
//
//    //   Tworzenie wektorow macierzy trojdiagonalnej
//        A[1]=0;
//        A[N]=(1-p[N]*h*0.5);
//        C[N]=0;
//        C[1]=(1+p[1]*h*0.5);
//        B[1]=(-2+q[1]*h*h);
//        B[N]=(-2+q[N]*h*h);
//
//        for(int i=2;i<N;i++)
//        {
//            A[i]=(1-p[i]*h*0.5);
//            C[i]=(1+p[i]*h*0.5);
//            B[i]=(-2+q[i]*h*h);
//        }
//
//    // Wyliczanie beta i gamma
//        beta[1]=-C[1]/B[1];
//        gamma[1]=D[1]/B[1];
//        for(int i=2;i<N+1;i++)
//        {
//            beta[i]=-C[i]/(A[i]*beta[i-1]+B[i]);
//            gamma[i]=(D[i]-A[i]*gamma[i-1])/(A[i]*beta[i-1]+B[i]);
//        }
//
//    // Wyliczanie rozwiazan
//        U[0]=x_a;
//        U[N+1]=x_b;
//        U[N]=gamma[N];
//
//        for(int i=N-1; i>0; i--)
//        {
//            U[i]=beta[i]*U[i+1]+gamma[i];
//        }