// PROGRAM laplace
// 長方形領域内のラプラス方程式をヤコビの緩和法で解く

#include "TrueBASIC.h"
int main();
void assign(double V[][101], int *nx, int *ny, double *min_change);
void iterate(double V[][101], int nx, int ny, double min_change, int *iterations);
void show_output(double V[][101], int nx, int ny, int iterations);

int main()
{
    int iterations, nx, ny;
    double V[101][101], min_change;

    assign(V, &nx, &ny, &min_change);
    iterate(V, nx, ny, min_change, &iterations);

    return 0;
}

void assign(double V[][101], int *nx, int *ny, double *min_change)
// 長方形領域の1辺の長さ
{
    int x, y;
    double V0;
    char Stmp1_[_LBUFF_];

    (*nx) = 9;                              // 領域内の x 方向の点の数
    (*ny) = 9;                              // 領域内の y 方向の点の数
    V0 = 10;                                // 長方形の境界のポテンシャル
    printf("percentage change = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", min_change);
    (*min_change) = (*min_change)/100;
    // 長方形の境界の電位を指定
    for(x = 0; x <= (*nx) + 1; ++x)
    {
        V[x][0] = V0;
        V[x][(*ny) + 1] = V0;
    }
    for(y = 0; y <= (*ny) + 1; ++y)
    {
        V[0][y] = V0;
        V[(*nx) + 1][y] = V0;
    }
    // 内部の点の電位の値を予想
    for(y = 1; y <= (*ny); ++y)
        for(x = 1; x <= (*nx); ++x)
            V[x][y] = 0.9*V0;
    show_output(V, *nx, *ny, 0);
}

void iterate(double V[][101], int nx, int ny, double min_change, int *iterations)
{
    int x, y;
    double Vave[101][101], change, dV;

    (*iterations) = 0;
    do
    {
        // 反復計算における格子点の電位の値の最大差
        change = 0;
        for(y = 1; y <= ny; ++y)
        {
            for(x = 1; x <= nx; ++x)
            {
                // 隣接したセルの電位の値を平均
                Vave[x][y] = V[x + 1][y] + V[x - 1][y];
                Vave[x][y] = Vave[x][y] + V[x][y + 1] + V[x][y - 1];
                Vave[x][y] = 0.25*Vave[x][y];
                // 電位の値の変化の割合を計算
                if(Vave[x][y] != 0)
                {
                    dV = fabs((V[x][y] - Vave[x][y])/Vave[x][y]);
                    if(dV > change)
                        change = dV;
                }
            }
        }
        for(y = 1; y <= ny; ++y)   // 各格子点の電位の値を更新
        {
            for(x = 1; x <= nx; ++x)
            {
                V[x][y] = Vave[x][y];
            }
        }
        ++(*iterations);
        show_output(V, nx, ny, *iterations);
    } while(!(change <= min_change));
}

void show_output(double V[][101], int nx, int ny, int iterations)
// 電位の値を表示
{
    int x, y;

    printf("\n");
    printf("iteration = %d\n", iterations);
    for(y = 0; y <= ny + 1; ++y)
    {
        for(x = nx + 1; x >= 0; --x)
            printf("%6.2f ", V[x][y]);
        printf("\n");
    }
}