# include <stdio.h>
# include <stdlib.h>
# include <math.h>
 
# define    Px_Y    3000           /* quantidade de píxels em y */
# define    Px_X    3000           /* quantidade de píxels em x */
# define    PlaX    1000           /* largura das placas */
# define    DDP     9.4            /* diferença de potencial entre as placas */
# define    Tol     0.00000000001  /* tolerância */
# define    Arq     "saida.dat"    /* nome do arquivo de saída */

FILE *saida; 

int main()
{
    int     i, j, h, l1, l2, l3, d, d2, k;
    double  acm, ant, aux, uni, ***matriz;
    
    /**************************************************************/
    
        /* Alocação da matriz principal */   
        matriz = malloc(Px_Y*sizeof(double**));
        for(i = 0; i < Px_Y; i++){
            matriz[i] = malloc(Px_X*sizeof(double*));
            for(j = 0; j < Px_X; j++)
                matriz[i][j] = malloc(2*sizeof(double)); 
        }
        
    /**************************************************************/
    
        /* Inicialização da matriz principal */
        for(i=0; i < Px_Y; i++)
            for(j=0; j < Px_X; j++){
                matriz[i][j][0] = 0;        /* 0 indica meio, 1 indica placa condutora */ 
                matriz[i][j][1] = DDP/2;    /* Armazena o potencial do ponto */
            }
    
    /**************************************************************/
    
        /* Contrução das barras */
        i = (int)Px_Y/2;
        j = (int)Px_X/2;
        
        d = (int)ceil(PlaX/30);
        h = (int)ceil(PlaX/10);
        l1= (int)ceil(PlaX/3);
        l2= (int)ceil(2*PlaX/3);
        l3= (int)floor(l2/h);
        d2=d+h+1;

        for(k=0; k<l1; k++){
            matriz[i+d][j-k-1-(int)(l1/2)][0] = 1; 
            matriz[i+d][j-k-1-(int)(l1/2)][1] = DDP;
            matriz[i-d][j-k-1-(int)(l1/2)][0] = 1; 
            matriz[i-d][j-k-1-(int)(l1/2)][1] = 0; 
        }   
        
        for(k=0; k<=l2; k++){
            if(k%l3==0)d2--;
            if(d2<d)d2=d;
            matriz[i+d2][j+l2-k-(int)(l1/2)][0] = 1; 
            matriz[i+d2][j+l2-k-(int)(l1/2)][1] = DDP;
            matriz[i-d2][j+l2-k-(int)(l1/2)][0] = 1; 
            matriz[i-d2][j+l2-k-(int)(l1/2)][1] = 0;            
        }
        
    /**************************************************************/        
    
        /* Aplica a Equação de Laplace */
        acm = 1;
        while(acm>Tol*DDP){
            acm = 0;
            for(i = 0; i < Px_Y; i++)
                for(j = 0; j < Px_X; j++)
                    if(matriz[i][j][0]==0){
                        ant = matriz[i][j][1];
                        aux = 0;
                        
                        if(i==0) aux = 0.25*matriz[Px_Y-1][j][1];
                        else     aux = 0.25*matriz[i-1][j][1];
                        
                        if(i==Px_Y-1) aux += 0.25*matriz[0][j][1];
                        else          aux += 0.25*matriz[i+1][j][1];
                        
                        if(j==0) aux += 0.25*matriz[i][Px_X-1][1];
                        else     aux += 0.25*matriz[i][j-1][1];
                        
                        if(j==Px_X-1) aux += 0.25*matriz[i][0][1];
                        else          aux += 0.25*matriz[i][j+1][1];
                        
                        matriz[i][j][1] = aux;
                        if(fabs(ant-aux)>acm) acm = fabs(ant-aux);
                    }
        }
        
    
    /**************************************************************/    
    
         /* Imprime matriz final */
        saida = fopen(Arq, "w"); 
        uni=15.0/PlaX;
        fprintf(saida, "Simulação do Campo Elétrico das Barras Defletoras do TRC\n\n");
        fprintf(saida, "Grade usada:         %d x %d píxels\nLargura das placas:  %d píxels \nEscala:              %.2f cm/píxel \nDDP entre as placas: %.2f V\n\n", Px_X, Px_Y, PlaX, uni, DDP);
        
        fprintf(saida, "\nComponentes próximas ao eixo x:\n\ny\\x ");
        for(j=0, aux=-uni*Px_X/2; j < Px_X; j++, aux+=uni)
            fprintf(saida, "%f ", aux);
        fprintf(saida, "\n%f: ", uni);
        for(j=(int)Px_Y/2+1, i=0; i < Px_X; i++)
            fprintf(saida, "%f ", matriz[j][i][1]);     
        fprintf(saida, "\n0: ");
        for(j=(int)Px_Y/2, i=0; i < Px_X; i++)
            fprintf(saida, "%f ", matriz[j][i][1]);     
        fprintf(saida, "\n%f: ", -uni);
        for(j=(int)Px_Y/2-1, i=0; i < Px_X; i++)
            fprintf(saida, "%f ", matriz[j][i][1]);         
        
        fprintf(saida, "\n\n\nGrade completa:\n\n");
        for(i = 0; i < Px_Y; i++){
            for(j = 0; j < Px_X; j++)
                fprintf(saida ,"%f ", matriz[i][j][1]);
            fprintf(saida ,"\n");
        }      
 
    /**************************************************************/
    
        /* Liberação das posições de memória */
        for(i = 0; i < Px_Y; i++){
            for(j = 0; j < Px_X; j++)
                free(matriz[i][j]);
            free(matriz[i]);
        }   
        free(matriz);
        
return 0;
}