#include<iostream>
#include<fstream>
#include<math.h>
#include<stdlib.h>
#include<bibli_fonctions.h>
using namespace std;

/* Calcul du potentiel electrostatique dans une enceinte avec conducteurs */

const int nx = 61, ny = 106, nmax = 2000;
const double eps = 1.e-3;

int main(){
  double** V;
  int i, j, n, cont;
  double Vp;
  fstream fich;

  V = (double**)malloc(nx * sizeof(double*));
  for (i = 0; i < nx; i++){
    V[i] = (double*)malloc(ny * sizeof(double));
    for (j = 0; j < ny; j++)
      if ((i == 20 || i == 40) && (j >= 55 && j <= 75))
	V[i][j] = 1000.;
      else
	V[i][j] = 0.;
  }
  for (n = 1; n <= nmax; n++){
    cont = 0;
    for (i = 1; i < nx-1; i++)
      for (j = 1; j < ny-1; j++)
	if (!((i == 20 || i == 40)
	      && ((j >= 30 && j <= 50) || (j >= 55 && j <= 75)))){
	  Vp = V[i][j];
	  V[i][j] = (V[i-1][j] + V[i][j-1] + V[i+1][j] + V[i][j+1]) / 4.;
	  if (fabs(V[i][j] - Vp) > eps)
	    cont = 1;
	}
    if (cont == 0)
      break;
  }
  cout << "Nombre d'iterations : " << n - cont << endl;
  
  fich.open("lentille_electro_py.res", ios::out);
  for (i = 0; i < nx; i++){
    for (j = 0; j < ny; j++)
      fich << V[i][j] << " ";
    fich << endl;
  }
  fich.close();
  for (i = 0; i < nx; i++)
    free(V[i]);
  free(V);

  ostringstream pyth;
  pyth << "figure(1, figsize=(16,6))\n"
       << "subplot(1,2,1, projection='3d')\n" 
  // 1,2 veut dire : une "matrice" de figures de 1 (vertical) * 2 (horizontal),
  // et le dernier 1 veut dire que celle-ci est la premiere sous-figure. 
       << "Z = loadtxt('lentille_electro_py.res')\n"
       << "X = arange(Z.shape[0])\n"
       << "Y = arange(Z.shape[1])\n"
       << "X,Y = meshgrid(X,Y)\n" 
       << "gca().plot_surface(X,Y,Z.T, rstride=1, cstride=1, linewidth=0, cmap=cm.spectral_r, antialiased=False)\n"
  // Les commandes pour 3D sont seulement definies comme methodes d'un objet 
  // et non pas egalement comme fonctions independantes (contrairement aux 
  // commandes 2D ou meme plus general en Python). C'est pour cela qu'on a 
  // besoin du gca() un peu moche, qui donne la sous-figure actuelle. Je 
  // suppose que ca changera dans une version future de Matplotlib.
  // A cause de conventions un peu bizarres il faut ici la transposee de Z
  // (Z.T) et non pas Z. cmap specifie le choix de couleurs. Voir la
  // documentation de plot_surface pour les autres options.
       << "xlabel('x')\n"
       << "ylabel('y')\n"
       << "title('Potentiel electrostatique')\n"
       << "subplot(1,2,2)\n" // courbes de niveau dans l'autre sous-figure
       << "cont = contour(Z.T, linspace(0,1000,11), cmap=cm.spectral_r)\n"
  // On a besoin de la variable cont pour la commande suivante, sinon on
  // pourrait simplement enlever le 'cont = '.
       << "clabel(cont, fmt='%d')\n" // rajouter les valeurs sur les courbes
       << "xlabel('x')\n"
       << "ylabel('y')\n"
       << "title('Courbes de niveau')\n"
    ;
  make_plot_py(pyth);

  return 0;
}