// PROGRAM interfere

#include "TrueBASIC.h"
int main();
void initial(int *nslit, double *a, double *L, double *lambda, double *ymax, 
             int *nbar);
void pattern(int nslit, double a, double L, double lambda, double ymax, int nbar);
int main()
{
    int nbar, nslit;
    double L, a, lambda, ymax;

    GWopen(0);

    initial(&nslit, &a, &L, &lambda, &ymax, &nbar);
    pattern(nslit, a, L, lambda, ymax, nbar);

    GWquit();
    return 0;
}

void initial(int *nslit, double *a, double *L, double *lambda, double *ymax, 
             int *nbar)
{
    int itick, ntick;
    double Imax, Lt, dy;
    char Stmp1_[_LBUFF_];

    printf("number of slits = " );
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", nslit);
    printf("distance between slits (mm) = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", a);
    printf("distance to the screen (mm) = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", L);
    printf("wavelength (angstroms)= ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", lambda);
    printf("max. position on photographic plate (mm) = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", ymax);
    printf("# times to average intensity = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", nbar);
    (*lambda) = (*lambda)*1e-7;             // オングストロームを mm に変換
    Imax = pow((*nslit)/(*L), 2.0);         // 強度の最大値
    GWindow((float)(-1.1*(*ymax)), (float)(-0.1*Imax), (float)(1.1*(*ymax)), 
            (float)(1.1*Imax));
    GWline((float)(-(*ymax)), 0.0f, (float)((*ymax)), 0.0f);
                                            // 横軸を描く
    GWline(0.0f, (float)(Imax), 0.0f, 0.0f);
                                            // 強度に対する軸を描く
    printf("intensity versus vertical distance on screen\n");
    dy = (*lambda)*(*L)/(*a);               // L >> aならば，極大値間の距離
    ntick = (int)((*ymax)/dy);
    Lt = 0.01*Imax;                         // 縦軸の刻み
    for(itick = -ntick; itick <= ntick; ++itick)
        GWline((float)(itick*dy), 0.0f, (float)(itick*dy), (float)(Lt));
}

void pattern(int nslit, double a, double L, double lambda, double ymax, int nbar)
{
    int islit, ipoint, npoint, t;
    double L2, amplitude, delta, dy, intensity, k, phase, r, r2, y, yslit;

    GWsetpen(-1, -1, -1, -1);
    k = 2*pi/lambda;                        // 波数
    npoint = 200;                           // y 軸の正の領域にある点の数
    dy = ymax/npoint;                       // プロットする点の間の距離
    delta = a/(nslit - 1);                  // スリット間の距離
    phase = 2*pi/nbar;
    L2 = L*L;
    for(ipoint = -npoint; ipoint <= npoint; ++ipoint)
    {
        y = ipoint*dy;                      // スクリーンの位置
        intensity = 0;                      // 光の強度      
        for(t = 1; t <= nbar; ++t)          // 振幅ではなく強度の平均
        {
            amplitude = 0;
            for(islit = 1; islit <= nslit; ++islit)
            {
                yslit = -0.5*a + (islit - 1)*delta; // スリットの位置
                yslit -= y;      // スリットとスクリーンの位置の縦の距離
                r2 = L2 + yslit*yslit;
                r = sqrt(r2);               // スリットのスクリーンからの距離
                amplitude += (1/r)*cos(k*r - phase*t);
            }
            intensity += amplitude*amplitude;
        }
        intensity /= nbar;
        GWline2((float)(y), (float)(intensity));
    }
}