// PROGRAM nuclear_decay
// 不安定核の崩壊のシミュレーション

#include "TrueBASIC.h"
int main();
void initial(int *N0, double *p, int *ntrial, int *tmax, int ncum[]);
void decay(int N0, double p, int ntrial, int tmax, int ncum[]);
void output(int ntrial, int tmax, int ncum[]);

int main()
{
    int N0, ncum[1001], ntrial, tmax;
    double p;

    initial(&N0, &p, &ntrial, &tmax, ncum);
    decay(N0, p, ntrial, tmax, ncum);
    output(ntrial, tmax, ncum);

    return 0;
}

void initial(int *N0, double *p, int *ntrial, int *tmax, int ncum[])
{
    int t;
    char Stmp1_[_LBUFF_];

    rnd(-1);
    printf("initial number of unstable nuclei = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", N0);
    printf("decay probability for unstable nucleus = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", p);
    printf("number of time intervals per trial = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", tmax);
    printf("number of trials = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", ntrial);
    for(t = 0; t <= (*tmax); ++t)
        ncum[t] = 0;                        // 不安定核の数を積算する
}

void decay(int N0, double p, int ntrial, int tmax, int ncum[])
{
    int N[5001], i, t, itrial, n_unstable;

    for(itrial = 1; itrial <= ntrial; ++itrial)
    {
        for(i = 1; i <= N0; ++i)
            N[i] = 1;                       // 不安定な場合に nucleus = 1 とする
        ncum[0] += N0;
        n_unstable = N0;                    // 不安定核の数
        for(t = 1; t <= tmax; ++t)
        {
            for(i = 1; i <= N0; ++i)
            {
                if(N[i] == 1)
                {
                    if(rnd(0) <= p)
                    {
                        N[i] = 0;           // 核の崩壊
                        --n_unstable;
                    }
                }
            }
            ncum[t] += n_unstable;          // データを積算する
        }
    }
}

void output(int ntrial, int tmax, int ncum[])
// 別のプログラムで読めるようにデータをファイルに出力する
{
    int t;
    FILE *FP1_;

    FP1_ = fopen("decay.dat", "w");
    fprintf(FP1_, "        time mean number of unstable nuclei\n");
    for(t = 0; t <= tmax; ++t)
        fprintf(FP1_, "%12d %12f\n", t, (double)ncum[t]/ntrial);
    fclose(FP1_);
}