// PROGRAM eigen

#include "TrueBASIC.h"
int main();
void parameters(double *V0, double *a, double *xmax, double *dx);
void plot_potential(double V0, double a, double xmax);
void Euler(double V0, double a, double dx, double xmax);
double V(double x, double V0, double a);

int main()
{
    double V0, a, dx, xmax;

    GWopen(0);

    parameters(&V0, &a, &xmax, &dx);
    plot_potential(V0, a, xmax);
    Euler(V0, a, dx, xmax);

    GWquit();
    return 0;
}

void parameters(double *V0, double *a, double *xmax, double *dx)
{
    char Stmp1_[_LBUFF_];

    printf("magnitude of well depth = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", V0);
    printf("half width of well = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", a);
    printf("step size = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", dx);
    printf("maximum value of x to be plotted = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", xmax);
}

void plot_potential(double V0, double a, double xmax)
// draw potential well
{
    int IC1_ = 1;

    GWvport(0.0f, 0.0f, GWaspect(1), (float)(0.4));
    GWsetpen(RED, -1, -1, -1);
    GWindow((float)(-1.05*xmax), (float)(-1.1*V0), (float)(1.05*xmax), 
            (float)(1.1*V0));
    GWsavevp(IC1_);
    GWsetogn(0, IC1_);
    GWline((float)(-xmax), 0.0f, (float)(-a), 0.0f); // horizontal line
    GWline((float)(-a), 0.0f, (float)(-a), (float)(-V0)); // vertical line
    GWmove2((float)(-a), (float)(-V0));
    GWline2((float)(a), (float)(-V0));
    GWline2((float)(a), 0.0f);
    GWline2((float)(xmax), 0.0f);
}

void Euler(double V0, double a, double dx, double xmax)
{
    int IC2_ = 2, parity;
    double E, d2phi, dphi, phi, phi_old, x, x_old;
    char Stmp1_[_LBUFF_];

    GWvport(0.0f, (float)(0.4), GWaspect(1), 1.0f);
    GWindow((float)(-1.01*xmax), (float)(-4.0), (float)(1.01*xmax), 
            (float)(4.0));
    GWsavevp(IC2_);
    GWsetogn(0, IC2_);
    printf("even or odd parity (1 or -1) = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", &parity);
    GWerase(IC2_, -1);
    do
    {
        printf("E = ");
        fgets(Stmp1_, _LBUFF_, stdin);
        sscanf(Stmp1_, "%lg", &E);
        if(E == 0)
            return;
        GWsetpen(BLUE, -1, -1, -1);
        if(parity == -1)
        {
            phi = 0;                        // initial values at x = 0
            dphi = 1;                       // first derivative
        }
        else
        {
            phi = 1;
            dphi = 0;
        }
        x = 0;
        do                                  // compute wave function
        {
            x_old = x;
            phi_old = phi;
            x += dx;
            d2phi = 2*(V(x, V0, a) - E)*phi;    // dimensionless units
            dphi += d2phi*dx;         // Euler-Cromer algorithm
            phi += dphi*dx;
            // plot wave function
            GWline((float)(x_old), (float)(phi_old), (float)(x), (float)(phi));
            GWline((float)(-x_old), (float)(phi_old*parity), (float)(-x), 
                   (float)(phi*parity));
        } while(x <= xmax);
    } while(1);
}

double V(double x, double V0, double a)
// potential function V(x)
{
    double V_;

    if(fabs(x) > a)
        V_ = V0;
    else
        V_ = 0;
    return V_;
}