// PROGRAM verify
// 離散型のラプラス方程式を正方格子上の
// 点電荷について検証する

#include "TrueBASIC.h"
#include "csgraphics.h"

int main();
void initial(double *L, double *dx, double *dy, int *IC1_, int *IC2_);
void draw_grid(double L, double dx, double dy, int *IC1_);
void potential(double L, double dx, double dy, char* stop_plot, int IC1_, 
               int IC2_, int IC3_);
void showpotential(double x, double y, int xp, int yp, double *V);

int main()
{
    int IC1_ = 1, IC2_ = 2, IC3_ = 3;
    double L, dx, dy;
    char stop_plot[_LBUFF_];

    GWopen(0);

    initial(&L, &dx, &dy, &IC1_, &IC2_);
    draw_grid(L, dx, dy, &IC3_);
    do
    {
        potential(L, dx, dy, stop_plot, IC1_, IC2_, IC3_);
    } while(!(strcmp(stop_plot, "yes") == 0));

    GWquit();
    return 0;
}

void initial(double *L, double *dx, double *dy, int *IC1_, int *IC2_)
{
    char Stmp1_[_LBUFF_];

    printf("lattice dimension = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", L);
    printf("grid spacing in x direction = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", dx);
    printf("grid spacing in y direction = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", dy);
    GWvport(0.0f, 0.8f, 0.25f*GWaspect(1), 1.0f);
    GWindow(-1.2f, -1.0f, 2.0f, 2.0f);
    GWsavevp(*IC1_);
    GWsetogn(0, (*IC1_));
    GWvport(0.0f, 0.0f, 0.25f*GWaspect(1), 0.6f);
    GWindow(0.0f, 24.0f, 80.0f, 0.0f);
    GWsavevp(*IC2_);
    GWsetogn(0, (*IC2_));
}

void draw_grid(double L, double dx, double dy, int *IC1_)
{
    double x, y;
    float a, xwin, ywin;

    GWvport(0.25f*GWaspect(1), 0.0f, GWaspect(1), 1.0f);
    compute_aspect_ratio((float)L, &xwin, &ywin);
    GWindow(-xwin, -ywin, xwin, ywin);
    GWsavevp(*IC1_);
    GWsetogn(0, (*IC1_));
    a = 0.2*dx;                             // 画面上の電荷の"半径"
    GWsetmrk(6, 2*a, BLUE, -1, -1);
    GWputmrk(0.0f, 0.0f);
    GWsetmrk(6, a/2, BLACK, -1, -1);
    for(y = -L; y <= L; y += dy)
        for(x = -L; x <= L; x += dx)
            GWputmrk((float)(x), (float)(y));
    GWrect((float)(-L), (float)(-L), (float)(L), (float)(L));
}

void potential(double L, double dx, double dy, char* stop_plot, int IC1_, 
               int IC2_, int IC3_)
{
    int Itmp1_;
    float xs, ys;
    double V0, V1, V2, V3, V4;
    char Stmp1_[_LBUFF_];

    GWselvp(IC3_);
    GWsetogn(0, IC3_);
    GWsetmsg("click on mouse to choose site\n");
    do
    {
    } while(!(GWmouse(&Itmp1_, &xs, &ys) || Itmp1_));
    GWselvp(IC1_);
    GWsetogn(0, IC1_);
    GWerase(IC1_, -1);
    strcpy(stop_plot, "");
    if(fabs(xs) <= L && fabs(ys) <= L)
    {
        showpotential((double)xs, (double)ys, 0, 0, &V0);
        showpotential(xs+dx, (double)ys, 1, 0, &V1);
        showpotential(xs-dx, (double)ys, -1, 0, &V2);
        showpotential((double)xs, ys+dy, 0, 1, &V3);
        showpotential((double)xs, ys-dy, 0, -1, &V4);
        GWselvp(IC2_);
        GWsetogn(0, IC2_);
        GWerase(IC2_, -1);
        GWputtxt(1.0f, 1.0f, " average potential");
        GWputtxt(1.0f, 2.0f, " of four neighbors:");
        sprintf(Stmp1_, "%.4f", 0.25*(V1+V2+V3+V4));
        GWputtxt(1.0f, 3.0f, Stmp1_);
    }
    else
        strcpy(stop_plot, "yes");
}

void showpotential(double x, double y, int xp, int yp, double *V)
{
    char Stmp1_[_LBUFF_];
    (*V) = 1/sqrt(x*x + y*y);               // 点電荷の電位
    sprintf(Stmp1_, "%.4f", (*V));
    GWputtxt((float)xp, (float)yp, Stmp1_);
}