#include "stdio.h"
#include "math.h"
#define MAX 80

void nextstep(float A[MAX][MAX], float B[MAX][MAX], float
f[MAX][MAX], float alfa,float beta, float Uext)
{int i,j,k,N=MAX;
 float M=MAX;
 float h=4/M,t=h*h/4;
 float c=2*h*alfa;
 int m=N/3;

B[0][0]=A[0][0]+t*f[0][0]+beta*(2*A[0][1]+2*A[1][0]-4*A[0][0]);
B[0][N-1]=A[0][N-1]+t*f[0][N-1]+beta*(2*A[1][N-1]+2*A[0][N-2]-4*A[0][N-1]);
B[N-1][0]=A[N-1][0]+t*f[N-1][0]+beta*(2*A[N-1][1]+2*A[N-2][1]-4*A[N-1][0]);
B[N-1][N-1]=A[N-1][N-1]+t*f[N-1][N-1]+beta*(2*A[N-2][N-1]+2*A[N-1][N-2]-4*A[N-1][N-1]);


for(k=1;k<=N-2;k++)
{B[0][k]=A[0][k]+t*f[0][k]-beta*(4*A[0][k]-2*A[1][k]-A[0][k+1]-A[0][k-1]);

B[N-1][k]=A[N-1][k]+t*f[N-1][k]-beta*(4*A[N-1][k]-2*A[N-2][k]-A[N-1][k+1]-A[N-1][k-1]);}

for(k=1;k<=m-1;k++)
{B[k][0]=A[k][0]+t*f[k][0]+beta*(2*A[k][1]+A[k-1][0]+A[k+1][0]-4*A[k][0]);

B[k][N-1]=A[k][N-1]+t*f[k][N-1]+beta*(2*A[k][N-2]+A[k+1][N-1]+A[k-1][N-1]-4*A[k][N-1]);}


for(k=2*m+1;k<=N-2;k++)
{B[k][0]=A[k][0]+t*f[k][0]+beta*(2*A[k][1]+A[k-1][0]+A[k+1][0]-4*A[k][0]);

B[k][N-1]=A[k][N-1]+t*f[k][N-1]+beta*(2*A[k][N-2]+A[k+1][N-1]+A[k-1][N-1]-4*A[k][N-1]);}

for (k=m;k<=2*m;k++)
{B[k][0]=A[k][0]+t*f[k][0]+beta*(2*A[k][1]+c*Uext+A[k+1][0]+A[k-1][0]-(4+c)*A[k][0]);
 B[k][N-1]=A[k][N-1];}

for(i=1;i<=N-2;i++)
{for(j=1;j<=N-2;j++)
{B[i][j]=A[i][j]+t*f[i][j]+beta*(A[i-1][j]+A[i+1][j]+A[i][j+1]+A[i][j-1]-4*A[i][j]);}
}}


void copier(float A[MAX][MAX],float B[MAX][MAX])
{int i,j,N=MAX;

for(i=0;i<=N-1;i++)
 {for(j=0;j<=N-1;j++)
 {B[i][j]=A[i][j];}
 }
}


void afficher(float A[MAX][MAX],int k)
{int i,j,N=MAX; float M=MAX;
float h=4/M;
/*-------------------------------------------------------*/
/*  Affichage d'un fichier de commandes Matlab           */
/*  Le programme genere automatiquement une suite        */
/*  d'instructions en Matlab qui stocke des images dans  */
/*  pour en faire une animation. Pour executer cette     */
/*  animation, compilez:   gcc warm.c -o warm            */
/*  redirigez dans un fichier  ./warm > animation.m      */
/*  puis rentrez dans Matlab et tapez:   A=moviein(48)   */
/*                               puis:  animation        */
/*-------------------------------------------------------*/
printf("x=linspace(0,4,%d);\n",N);
printf("y=x;\n");
 printf("T%d=[",k);
for(i=0;i<=N-2;i++)
{for(j=0;j<=N-1;j++)
{printf("%3.1f ",A[i][j]);}
 printf(";\n");}
for(j=0;j<=N-1;j++)
{printf("%3.1f ",A[i][j]);}
printf ("];\n");
printf("surf(x,y,T%d)\n",k);
printf("axis([0 4 0 4 12 26])\n");
printf("zlabel('TEMPERATURE')\n");
printf("title('t=%f heure')\n",0.5*k);
/* on stocke une image toute les demi-heures*/
 printf("A(:,%d)=getframe;\n",k);}

main()
{float A[MAX][MAX];
float B[MAX][MAX];
float f[MAX][MAX];
float C[MAX][MAX];
float alfa=0.1;

int i,j,k,P,N=MAX,m=N/3;
 float M=MAX;
float h=4/M,t=h*h/4;
P=1/t;
/* 1600 iterations correspondent a une heure */
for(i=0;i<=N-1;i++)
{for(j=0;j<=N-1;j++)
 {A[i][j]=15;

 f[i][j]=0;}}

 for(k=m;k<2*m;k++)
   {A[k][N-1]=25;}

for(k=1;k<=48;k++)
     /* une iteration de cette boucle correspond a 1/2 heure */
{for(i=1;i<=800;i++)
{nextstep(A,B,f,alfa,0.25,-10);
 copier(B,A); }
 afficher(A,k);
 printf("\n");}
}