// PROGRAM scatter
// 微分断面積と全断面積の計算

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

int main();
void initial(double *x0, double *vx0, double *R, double *delta_b, double *dt, 
             int *nbin, double *dtheta);
void trajectory(double *vx, double *vy, double vx0, double x0, double b, 
                double dt);
void detector(double dN[], double vx, double vy, double b, double dtheta, 
              double *N);
void EulerRichardson(double *x, double *y, double *vx, double *vy, double dt);
void force(double x, double y, double *fx, double *fy);
void output(double dN[], double N, int nbin, double dtheta, double R);
double arctan(double x, double y);

int main()
{
    int nbin, Itmp1_;
    double N, R, b, dN[361], delta_b, dt, dtheta, vx, vx0, vy, x0;
    for(Itmp1_ = 0; Itmp1_ < 361; ++Itmp1_) dN[Itmp1_] = 0;

    GWopen(0);

    initial(&x0, &vx0, &R, &delta_b, &dt, &nbin, &dtheta);
    N = 0;                                  // 入射粒子の数
    b = 0;
    // b = R としていても，b が R をわずかに越える可能性があることに注意
    while(b <= R + delta_b/2)
    {
        trajectory(&vx, &vy, vx0, x0, b, dt);
        detector(dN, vx, vy, b, dtheta, &N);
        b += delta_b;
    }
//  printf("Hit any key for the differential cross section.\n");
//  do
//  {
//      GWsleep(10);
//  } while(!TBkeyinput());
    output(dN, N, nbin, dtheta, R);

    GWquit();
    return 0;
}

void initial(double *x0, double *vx0, double *R, double *delta_b, double *dt, 
             int *nbin, double *dtheta)
{
    double E;
    float Rs, xwin, ywin;
    char Stmp1_[_LBUFF_];

    printf("kinetic energy of incident particles = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", &E);
    printf("incremental change in impact parameter = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", delta_b);
    (*R) = 1.0;                             // ビームの半径
    (*vx0) = sqrt(2*E);                     // 質量 = 1
    (*dt) = 0.01;                           // 分割時間
    (*nbin) = 9;                            // 検出器の数
    (*dtheta) = 180.0/(*nbin);              // 角度
    // 目標物から測ったビームの初期位置
    (*x0) = -3;
    compute_aspect_ratio((float)fabs(*x0), &xwin, &ywin);
    GWindow(-xwin, -ywin, xwin, ywin);
    GWsetpen(BLUE, -1, -1, -1);
    Rs = 1.0;                               // 散乱領域の半径
    GWellipse(-Rs, -Rs, Rs, Rs);
    GWsetpxl(0.0f, 0.0f, BLUE);
}

void trajectory(double *vx, double *vy, double vx0, double x0, double b, 
                double dt)
{
    double x, y;

    y = b;
    x = x0;
    (*vx) = vx0;
    (*vy) = 0;
    do
    {
        EulerRichardson(&x, &y, vx, vy, dt); // 分割時間の1ステップ
        GWsetpxl((float)(x), (float)(y), RED);
    } while(!(x*x + y*y > x0*x0 + b*b));
}

void detector(double dN[], double vx, double vy, double b, double dtheta, 
              double *N)
{
    int theta;

    theta = (int)(arctan(vx, vy)/dtheta);   // 散乱角
    dN[theta] += b;
    (*N) += b;
}

void EulerRichardson(double *x, double *y, double *vx, double *vy, double dt)
{
    double fx, fxmid, fy, fymid, vxmid, vymid, xmid, ymid;

    // 配列を使うとプログラムの行数を減らすことができる
    force(*x, *y, &fx, &fy);
    // プログラム中では質量は1である
    vxmid = (*vx) + 0.5*fx*dt;
    vymid = (*vy) + 0.5*fy*dt;
    xmid = (*x) + 0.5*(*vx)*dt;
    ymid = (*y) + 0.5*(*vy)*dt;
    force(xmid, ymid, &fxmid, &fymid);
    (*vx) += fxmid*dt;
    (*vy) += fymid*dt;
    (*x) += vxmid*dt;
    (*y) += vymid*dt;
}

void force(double x, double y, double *fx, double *fy)
{
    double f, r, r2;

    r2 = x*x + y*y;
    r = sqrt(r2);
    // 簡単な水素原子のモデルの場合の相互作用
    // 力の作用する範囲を 1 に設定
    if(r2 <= 1)
        f = 1/r2 - r;
    else
        f = 0;
    (*fx) = f*x/r;
    (*fy) = f*y/r;
}

void output(double dN[], double N, int nbin, double dtheta, double R)
{
    int i;
    double d_Omega, diffcross, dtheta_rad, target_density, theta, totalcross = 0;

//  GWclear(-1);
    dtheta_rad = dtheta*pi/180;             // ラジアン
    target_density = 1/(pi*R*R);
    printf("\n       theta         dN/N   sigma(theta)\n\n");
    for(i = 0; i <= nbin; ++i)
    {
        theta = (i + 0.5)*dtheta_rad;
        if(dN[i] > 0)
        {
            d_Omega = 2*pi*sin(theta)*dtheta_rad;
            diffcross = dN[i]/(d_Omega*N*target_density);
            totalcross = totalcross + diffcross*d_Omega;
            printf("%12f %12f %12f\n", (i+1)*dtheta, dN[i]/N, diffcross);
        }
    }
    printf("\ntotal cross section = %f\n", totalcross);
}

double arctan(double x, double y)
{
    double arctan_, theta;

    // -pi/2 ラジアンから pi/2 ラジアンを 0 度から 180 度に変換
    if(x > 0)
        theta = atan(fabs(y/x));
    else if(y > 0)
        theta = atan(y/x) + pi;
    else
        theta = pi - atan(y/x);
    arctan_ = theta*180/pi;
    return arctan_;
}