// PROGRAM fermat
// Monte Carlo method for finding minimum optical path

#include "TrueBASIC.h"
int main();
void initial(float y[], float v[], int *N);
void change_path(float y[], float v[], int N);

int main()
{
    int N;
    float v[101], y[101];

    GWopen(0);

    initial(y, v, &N);
    change_path(y, v, N);

    GWquit();
    return 0;
}

void initial(float y[], float v[], int *N)
{
    int i;
    float index, v1, v2;

    rnd(-1);
    (*N) = 10;                              // number of different media (even)
    // speed of light in vacuum equal to unity
    v1 = 1.0;
    index = 1.5;                            // index of refraction of medium #2
    v2 = 1/index;
    for(i = 0; i < (*N)/2; ++i)
        v[i] = v1;                          // speed of light in left half
    for(i = (*N)/2; i < (*N); ++i)
        v[i] = v2;                          // speed of light in right half
    y[0] = 2;                               // fixed source
    y[(*N)-1] = 8;                          // fixed detector
    GWindow(-0.5f, y[0]-0.5f, *N+1.0f, y[(*N)-1]+0.5f);
    GWrect(1.0f, y[0], (float)(*N), y[(*N)-1]);
    GWline((float)(*N)/2, y[0], (float)(*N)/2, y[(*N)-1]);
                                            // boundary
    // choose initial path at boundary between endpoints
    for(i = 1; i < (*N)-1; ++i) {
        y[i] = (y[(*N)-1] - y[0])*rnd(0) + y[0];
    }
    GWsetpen(RED, -1, -1, -1);
    GWanchor(1);
    GWsetogn(0, 10);
    GWplot1(-1, *N, 1.0f, (float)(*N), 0.0f,0.0f,1.0f,0.0f, y);
    GWsetogn(0, -10);                      // initial path
}

void change_path(float y[], float v[], int N)
{
    int x;
    float delta, dist, dy2, t_original, t_trial, ytrial;

    delta = 0.5;                           // maximum change in y
    do
    {
        // choose random x not including x = 1 or N
        x = (int)((N-2)*rnd(0)) + 1;
        ytrial = y[x] + (2*rnd(0) - 1)*delta;  // new y position
        dy2 = pow(y[x] - y[x + 1], 2.0);   // vertical distance squared
        dist = sqrt(1 + dy2);              // horizontal distance is unity
        t_original = dist/v[x + 1];
        dy2 = pow(y[x] - y[x - 1], 2.0);
        dist = sqrt(1 + dy2);
        t_original += dist/v[x];
        dy2 = pow(ytrial - y[x + 1], 2.0);
        dist = sqrt(1 + dy2);
        t_trial = dist/v[x + 1];
        dy2 = pow(ytrial - y[x - 1], 2.0);
        dist = sqrt(1 + dy2);
        t_trial += dist/v[x];
        if(t_trial < t_original)           // new position reduces time
        {
            y[x] = ytrial;
            GWplot1(-1, N, 1.0f, (float)N, 0.0f,0.0f,1.0f,0.0f, y);
            GWflush(-10);
        }
    } while(!TBkeyinput());
}