// PROGRAM fermat
// 最短の光学的経路を求めるためのモンテカルロ法

#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;                              // 異なる媒質の数(偶数)
    // 真空中の光速を 1 にとる
    v1 = 1.0;
    index = 1.5;                            // 第2の媒質の屈折率
    v2 = 1/index;
    for(i = 0; i < (*N)/2; ++i)
        v[i] = v1;                          // 左半分における光の速度
    for(i = (*N)/2; i < (*N); ++i)
        v[i] = v2;                          // 右半分における光の速度
    y[0] = 2;                               // 固定された光源
    y[(*N)-1] = 8;                          // 固定された検出器
    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]);
                                            // 境界
    // 最初に境界で端点を結ぶ経路をとる
    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);                      // 最初の経路
}

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

    delta = 0.5;                           // y の変化の最大値
    do
    {
        // x = 1 と N 以外の乱数 x を選ぶ
        x = (int)((N-2)*rnd(0)) + 1;
        ytrial = y[x] + (2*rnd(0) - 1)*delta;  // 新しい y の位置
        dy2 = pow(y[x] - y[x + 1], 2.0);   // 鉛直方向の距離の2乗
        dist = sqrt(1 + dy2);              // 水平方向の距離は 1
        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)           // 新しい位置では時間が短くなる
        {
            y[x] = ytrial;
            GWplot1(-1, N, 1.0f, (float)N, 0.0f,0.0f,1.0f,0.0f, y);
            GWflush(-10);
        }
    } while(!TBkeyinput());
}