// PROGRAM maxwell
// イー-ビッシャーのアルゴリズム

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

int main();
void initial(double *escale, double *bscale, double *jscale);
void current();
void newE();
void newB();
void screen(int *IC1_, int *IC2_, int *IC3_);
void plotfields(double escale, double bscale, double jscale, int IC1_, int IC2_, 
                int IC3_);
void plotarrow(double V, int x, int y, double scale, double shiftx, 
               double shifty, double angle1, double angle2);
void rotate(float *x, float *y, float angle);

double B[22][22][22][4];
double E[22][22][22][4];
double dl;
double dl_1;
double dt;
double fpi;
double j[22][22][22][4];
int mid;
int n[4];

int main()
{
    int IC1_ = 1, IC2_ = 2, IC3_ = 3, k;
    double bscale, escale, jscale;

    GWopen(0);

    initial(&escale, &bscale, &jscale);
    screen(&IC1_, &IC2_, &IC3_);
    do
    {
        current();
        newE();
        newB();
        plotfields(escale, bscale, jscale, IC1_, IC2_, IC3_);
        do
        {
        } while(!TBkeyinput());
        k = TBgetkey();
    } while(k != 's');

    GWquit();
    return 0;
}

void initial(double *escale, double *bscale, double *jscale)
{
    int comp, x, y, z;
    double dl;

    dt = 0.03;
    n[1] = 8;
    n[2] = 8;
    n[3] = 8;
    mid = n[1]/2;
    dl = 0.1;
    dl_1 = 1/dl;
    fpi = 4*pi;
    (*escale) = dl/(4*pi*dt);
    (*bscale) = (*escale)*dl/dt;
    (*jscale) = 1;
    for(x = 1; x <= n[1]; ++x)
    {
        for(y = 1; y <= n[2]; ++y)
        {
            for(z = 1; z <= n[3]; ++z)
            {
                for(comp = 1; comp <= 3; ++comp)
                {
                    E[x][y][z][comp] = 0;
                    B[x][y][z][comp] = 0;
                }
            }
        }
    }
}

void current()
{
    // t = 0 に流れはじめた，x-y 面内の定常的なループ電流
    j[mid][mid][mid][2] = 1;
    j[mid][mid][mid][1] = -1;
    j[mid - 1][mid][mid][2] = -1;
    j[mid][mid - 1][mid][1] = 1;
}

void newE()
// 面上で定義された E
{
    int x, y, z;
    double curlBx, curlBy, curlBz;

    for(x = 1; x <= n[1]; ++x)
    {
        for(y = 1; y <= n[2]; ++y)
        {
            for(z = 1; z <= n[3]; ++z)
            {
                curlBx = B[x][y][z][2]+B[x][y + 1][z][3]-B[x][y][z + 1][2]-
                         B[x][y][z][3];
                curlBx = curlBx*dl_1;
                E[x][y][z][1] = E[x][y][z][1] + dt*(curlBx - fpi*j[x][y][z][1]);
                curlBy = B[x][y][z][3]-B[x][y][z][1]+B[x][y][z + 1][1]-
                         B[x + 1][y][z][3];
                curlBy = curlBy*dl_1;
                E[x][y][z][2] = E[x][y][z][2] + dt*(curlBy - fpi*j[x][y][z][2]);
                curlBz = B[x][y][z][1]+B[x + 1][y][z][2]-B[x][y + 1][z][1]-
                         B[x][y][z][2];
                curlBz = curlBz*dl_1;
                E[x][y][z][3] = E[x][y][z][3] + dt*(curlBz - fpi*j[x][y][z][3]);
            }
        }
    }
}

void newB()
// 辺上で定義された B
{
    int x, y, z;
    double curlEx, curlEy, curlEz;

    for(x = 1; x <= n[1]; ++x)
    {
        for(y = 1; y <= n[2]; ++y)
        {
            for(z = 1; z <= n[3]; ++z)
            {
                curlEx = E[x][y][z][3]-E[x][y][z][2]-E[x][y - 1][z][3]+
                         E[x][y][z - 1][2];
                curlEx = curlEx*dl_1;
                B[x][y][z][1] = B[x][y][z][1] - curlEx*dt;
                curlEy = E[x][y][z][1]-E[x][y][z][3]-E[x][y][z - 1][1]+
                         E[x - 1][y][z][3];
                curlEy = curlEy*dl_1;
                B[x][y][z][2] = B[x][y][z][2] - curlEy*dt;
                curlEz = E[x][y][z][2]-E[x][y][z][1]-E[x - 1][y][z][2]+
                         E[x][y - 1][z][1];
                curlEz = curlEz*dl_1;
                B[x][y][z][3] = B[x][y][z][3] - curlEz*dt;
            }
        }
    }
}

void screen(int *IC1_, int *IC2_, int *IC3_)
{
    int IC4_ = 4;
    float xwin, ywin;

    compute_aspect_ratio((float)(n[1]+1), &xwin, &ywin);
    GWcolor(BLACK, CLR_BACKGRND);
    GWcolor(WHITE, CLR_TEXT);
    GWvport(0.0f, 0.0f, 0.5f*GWaspect(1), 0.5f);
    GWindow(0.0f, 0.0f, xwin, ywin);
    GWsavevp(*IC1_);
    GWsetogn(0, *IC1_);
    GWvport(0.5f*GWaspect(1), 0.5f, GWaspect(1), 1.0f);
    GWindow(0.0f, 0.0f, xwin, ywin);
    GWsavevp(*IC2_);
    GWsetogn(0, *IC2_);
    GWvport(0.5f*GWaspect(1), 0.0f, GWaspect(1), 0.5f);
    GWindow(0.0f, 0.0f, xwin, ywin);
    GWsavevp(*IC3_);
    GWsetogn(0, *IC3_);
    GWvport(0.0f, 0.5f, 0.5f*GWaspect(1), 1.0f);
    GWindow(0.0f, 24.0f, 80.0f, 0.0f);
    GWsavevp(IC4_);
    GWsetogn(0, IC4_);
    GWputtxt(1.0f, 2.0f, "Type s to stop");
    GWputtxt(1.0f, 4.0f, "Type any other key for next time step");
}

void plotfields(double escale, double bscale, double jscale, int IC1_, int IC2_, 
                int IC3_)
{
    int x, y, z;
    float X1, Y1, X2, Y2;

    GWselvp(IC1_);
    GWsetogn(0, IC1_);
    GWerase(IC1_, -1);
    GWgetwn(&X1, &Y1, &X2, &Y2);
    GWputtxt((float)(X1+(X2-X1)/10), (float)(Y2-(Y2-Y1)/10), "E(x,y)");
    for(x = 1; x <= n[1]; ++x)
    {
        for(y = 1; y <= n[2]; ++y)
        {
            plotarrow(E[x][y][mid][1], x, y, escale, 0, 0.5, 0, pi);
            plotarrow(E[x][y][mid][2], x, y, escale, 0.5, 0, pi/2.0, 3*pi/2);
        }
    }
    GWselvp(IC2_);
    GWsetogn(0, IC2_);
    GWerase(IC2_, -1);
    GWgetwn(&X1, &Y1, &X2, &Y2);
    GWputtxt((float)(X1+(X2-X1)/10), (float)(Y2-(Y2-Y1)/10), "B(x,z)");
    for(x = 1; x <= n[1]; ++x)
    {
        for(z = 1; z <= n[3]; ++z)
        {
            plotarrow(B[x][mid][z][1], x, z, bscale, 0.5, 0, 0, pi);
            plotarrow(B[x][mid][z][3], x, z, bscale, 0, 0.5, pi/2.0, 3*pi/2);
        }
    }
    GWselvp(IC3_);
    GWsetogn(0, IC3_);
    GWerase(IC3_, -1);
    GWgetwn(&X1, &Y1, &X2, &Y2);
    GWputtxt((float)(X1+(X2-X1)/10), (float)(Y2-(Y2-Y1)/10), "j(x,y)");
    for(x = 1; x <= n[1]; ++x)
    {
        for(y = 1; y <= n[2]; ++y)
        {
            plotarrow(j[x][y][mid][1], x, y, jscale, 0, 0.5, 0, pi);
            plotarrow(j[x][y][mid][2], x, y, jscale, 0.5, 0, pi/2.0, 3*pi/2);
        }
    }
}

void plotarrow(double V, int x, int y, double scale, double shiftx, 
               double shifty, double angle1, double angle2)
{
    int i;
    float xx[] = {-0.25f, 0.25f, 0.12f, 0.25f,  0.12f}, 
          yy[] = { 0.0f,  0.0f,  0.12f, 0.0f,  -0.12f};
    if(V > 0)
    {
        for(i = 0; i < 5; ++i) {
            xx[i] *= (float)(V/scale);
            yy[i] *= (float)(V/scale);
            rotate(&(xx[i]), &(yy[i]), (float)angle1);
            xx[i] += (float)(x+shiftx);
            yy[i] += (float)(y+shifty);
        }
    }
    else if(V < 0)
    {
        for(i = 0; i < 5; ++i) {
            xx[i] *= (float)(-V/scale);
            yy[i] *= (float)(-V/scale);
            rotate(&(xx[i]), &(yy[i]), (float)angle2);
            xx[i] += (float)(x+shiftx);
            yy[i] += (float)(y+shifty);
        }
    }
    GWsetpen(YELLOW, -1, -1, -1);
    GWmove2(xx[0], yy[0]);
    GWline2(xx[1], yy[1]);
    GWline2(xx[2], yy[2]);
    GWmove2(xx[3], yy[3]);
    GWline2(xx[4], yy[4]);
    GWsetpen(WHITE, -1, -1, -1);
}

void rotate(float *x, float *y, float angle)
{
    float tx;
    tx = (*x)*cos(angle) - (*y)*sin(angle);
    *y = (*x)*sin(angle) + (*y)*cos(angle);
    *x = tx;
}