/* Equation de Helmholtz en 2D:

f_xx(x,y)+f_yy(x,y)+M2*f(x,y)=0, f(0,y)=Vleft, f(Lx,y)=Vright, f(x,0)=Vdown, f(x,Ly)=Vup


Solution par les méthodes Jacobi/Gauss-Seidel/Successive Over relaxation 

La bonne forme de la récurrence (e.g. pour SOR Jacobi) est:  

f_new[i][j]=(1-w)*f_old[i][j]+w*( cx*(f_old[i+1][j]+f_old[i-1][j])+
                                cy*(f_old[i][j+1]+f_old[i][j-1])+
                                M2Mcr*f_old[i][j] )



*/

//======================== Fichiers d'en-tête ===========================
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//========================================================================

//======================= Déclarations ===================================
#define sqr(x) (x)*(x)
#define Vright 0.0 
#define Vleft  0.0
#define Vup    0.0
#define Vdown  1.0

// Les valeurs Vright, Vleft, Vup et Vdown sont fournies à titre d'exemple
// Les modifier, le cas échéant, pour les rendre conformes au texte du TP!

#define LX 2.0   //ces valeurs sont données 
#define LY 10.0  //à titre d'exemple!
#define NX 400   //Les modifier, 
#define NY 2000  //le cas échéant!
#define ax (LX)/(NX) // pas de réseau en x
#define ay (LY)/(NY) // pas de réseau en y
#define M  10000     // nombre d'itérations
#define size (NX)*(NY) // taille du réseau
#define inv_ax 1.0/(ax) 
#define inv_ay 1.0/(ay)
#define cx (0.5*sqr(inv_ax)/( sqr(inv_ax)+sqr(inv_ay) ))
#define cy (0.5*sqr(inv_ay)/( sqr(inv_ax)+sqr(inv_ay) ))
#define M2 0.0 //la valeur 0 donne Laplace
#define M2Mcr (M2)*0.5/( sqr(inv_ax)+sqr(inv_ay) ) 




double f_old[NX+1][NY+1]; // On déclare ces tableaux 
double f_new[NX+1][NY+1]; // comme des variables globales

//============================================================================

//============================ Fonctions ===================================

void bc_up(){ // Conditions au bord supérieur, y=LY=NY*ay
// A remplir!
}

void bc_down(){// Conditions au bord inférieur, y=0
//A remplir
}


void bc_right(){//Conditions au bord à droite, x=LX=NX*ax 
//A remplir
}

void bc_left(){//Conditions au bord à gauche, x=0
//A remplir
}

void init(){// Valeurs initiales pour la solution 
    int i,j;
    for(i=1;i<NX;i++){
	for(j=1;j<NY;j++){
	    f_old[i][j]=0.0; f_new[i][j]=0.0;}//initial guess
    }
}




void sor_update(){
//A remplir
}


void output(){// Cette fonction affiche la configuration finale  
              // sous la forme x y f_old[i][j]
}



double error(){// Cette fonction affiche l'écart entre les configurations
               // nouvelle et ancienne d'une itération
               // ATTENTION pour Gauss-Seidel et SOR à bien définir f_new!!
// A remplir
}


double energy_bord(){// Cette fonction calcule la contribution des bords 
                     // à la fonctionnelle d'énergie
// A remplir
}


double energy_bulk(){// Cette fonction calcule la contribution de l'intérieur
                     // du domaine à la fonctionnelle d'énergie
// A remplir
}


main(){
// Déclarer des variables locales nécessaires pour cette fonction 
// et décommenter les commandes appropriées


    init(); bc_up(); bc_down(); bc_right(); bc_left(); 

    for(k=1;k<M;k++){//Itérations 

//	sor_update();

//      Afficher l'erreur 
/*      printf("%d\t%g\n",k,sqrt(error())/(double)size);*/

//      Afficher la valeur de l'énergie

/*	printf("%g %g %g\n",energy_bulk()*ax*ay,ax*energy_bord(),energy_bulk()*ax*ay+energy_bord()*ax);*/
	} 
 
/* afficher la configuration finale */


//    output();



}