#include <winbgim.h>
#include <stdio.h>
#include <math.h>

/* Resuelve la ecuación de laplace entre dos cilindros
 de radios 10 y 50 unidades arbitrarias
 Pixelado en un cuadrado de 100 X 100 puntos
 Las condiciones de contorno son los potenciales fijados de 
 100 V para r<=10 y 0 V para r >=50
 Para comparar, este problema admite solución analitica que es 
 V(r)= 100*ln(50/r)/ln(50/10)
*/ 
main()
{
FILE *filout;
// Resuelve la ecuaci¢n de Laplace en un circulo
int i,j,k,l=0,m, Nk=20; char c;
int color;
float N[101][101];float d,p;
/* datos iniciales*/
initwindow(800,600);         // Abre una ventana gráfica de  400x300 pixeles

for(i=0;i<=100;i++)
{
 for(j=0;j<=100;j++)
 { N[i][j]=1.;
   m=(50-i)*(50-i)+(50-j)*(50-j);
   if(m<=100) N[i][j]=100.;
   if(m>=2500) N[i][j]=0.;
 }
}

printf("Comienzo del programa\n");
getchar();
do
{ // Escribe en pantalla cada 10 unidades

  for(i=0;i<=100;i=i+10)
  {
      for(j=0;j<=100;j=j+10)
      {
        printf("%3.0f ",N[i][j]);
      }
      printf("\n");
 
  }
      printf("\n\n");
// Dubuja lo calculado hasta este momento
      for(i=0;i<=100;i++)for(j=0;j<=100;j++)
      { if(N[i][j]<=0)color=0;
        if(N[i][j]>0&&N[i][j]<=20)color=1;
        if(N[i][j]>20&&N[i][j]<=40)color=2;
        if(N[i][j]>40&&N[i][j]<=60)color=3;
        if(N[i][j]>60&&N[i][j]<=80)color=4;
        if(N[i][j]>80&&N[i][j]<100)color=9;
        if(N[i][j]>=100)color=15; 
        putpixel(200+4*i,80+4*j,color);
       }      
c=getchar();

//closegraph();                       // Cierra la ventana gráfica 

  for(k=1;k<=Nk;k++) // Hace ciclos de Nk pasos
  {
    for(i=1;i<100;i++)
    {
      for(j=1;j<100;j++)
      { 
        m=(50-i)*(50-i)+(50-j)*(50-j);       
        if(m<=100||m>=2500) continue;
        N[i][j]=(N[i][j+1]+N[i][j-1]+N[i-1][j]+N[i+1][j])/4.;
     }
     }
  }
    
l++;
  
//printf("\n Ciclos: %d  para parar dar 0",l);
//scanf("%c",&c);

}while(!kbhit()&&c!='0');                    // Espera a que el usuario presione 0 


filout=fopen("laplace.dat","w");
    for(i=0;i<=100;i=i+10)
    {
      for(j=0;j<=100;j=j+10)
      {
       fprintf(filout,"%7.2f ",N[i][j]);
      }
      fprintf(filout,"\n");
    }

fclose(filout);
getchar(); // espera a cerrar el programa y el dibujo
}