// PROGRAM galaxy
// Schulman と Seiden によって提案されたモデル

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

int main();
void parameter();
void initial();
void evolve();
void create(int r, int a, int *newactive, int newactive_r[], int newactive_a[]);
void plot_spiral();

int active_a[10001];
int active_r[10001];
int cell[51][301];
int dt;
int nactive;
int nring;
double p;
double s;
int t;
double two_pi;
double v;
int ICMB;

int main()
{
    GWopen(0);

    parameter();
    initial();
    do
    {
        evolve();
        plot_spiral();
        t += dt;
    } while(!TBkeyinput());

    GWquit();
    return 0;
}

void parameter()
{
    nring = 50;                             // リングの数
    nactive = 200;                          // 初期の活性なセルの数
    v = 1;                                  // 角速度
    p = 0.18;                               // 星の形成確率
    dt = 10;                                // 分割時間
    s = 2*pi/6;                             // セルの幅
    two_pi = 2*pi;
}

void initial()
{
    int a, i, r, x, y, Itmp1_, Itmp2_;
    float xwin, ywin;

    rnd(-1);
    for(Itmp1_ = 0; Itmp1_ <= 50; ++Itmp1_)
        for(Itmp2_ = 0; Itmp2_ <= 300; ++Itmp2_)
            cell[Itmp1_][Itmp2_] = 0;
    // ランダムに活性化されたセル
    i = 0;
    do
    {
        do
        {
            x = (int)(nring*rnd(0)) + 1;
            y = (int)(nring*rnd(0)) + 1;
            r = (int)(sqrt(x*x + y*y)) + 1;
        } while(!(r <= nring));
        a = (int)(6*r*rnd(0)) + 1;          // 角度に対応した配列の指数
        if(cell[r][a] == 0)
        {
            ++i;
            active_a[i] = a;                // 活性な領域の位置
            active_r[i] = r;
            // 活性領域, 15 分割時間の間活性
            cell[r][a] = 15;
        }
    } while(!(i == nactive));
    t = 0;                                  // 初期時間
    compute_aspect_ratio((float)nring, &xwin, &ywin);
    GWindow(-xwin, -ywin, xwin, ywin);
//    GWcolor(WHITE, CLR_BACKGRND);
//    GWcolor(BLACK, CLR_TEXT);
    GWcolor(BLACK, CLR_BACKGRND);
    GWcolor(WHITE, CLR_TEXT);
    GWsetpen(-8, -1, -1, 4);
    GWsetbrs(-8, 1, -1);
    GWbegincmb(0, 1.0, 1.0);
    GWellipse(0.0,0.0,1.0,1.0);
    ICMB = GWendcmb();
}

void evolve()
{
    int a, am, ap, i, newactive, newactive_a[10001], newactive_r[10001], r, 
        Itmp1_;
    double angle, wt;

    // 次の単位時間の間活性な星のクラスターの数
    newactive = 0;
    for(i = 1; i <= nactive; ++i)
    {
        r = active_r[i];
        a = active_a[i];
        // 同じリングの中で活性な隣のセルの数
        create(r, a+1, &newactive, newactive_r, newactive_a);
        create(r, a-1, &newactive, newactive_r, newactive_a);
        angle = fmod((double)(a*s + v*t)/r, two_pi);
        if(r < nring)                       // 次の大きなリングの中の活性なセル
        {
            wt = fmod((double)v*t/(r+1), two_pi);
            ap = (int)((angle - wt)*(r+1)/s);
            ap = ((ap) % (6*(r+1)));
            create(r+1, ap, &newactive, newactive_r, newactive_a);
            create(r+1, ap+1, &newactive, newactive_r, newactive_a);
        }
        if(r > 1)                           // 次の小さなリングの中の活性なセル
        {
            wt = fmod(v*t/(r-1), two_pi);
            am = (int)((angle - wt)*(r-1)/s);
            am = ((am) % (6*(r-1)));
            create(r-1, am, &newactive, newactive_r, newactive_a);
            create(r-1, am+1, &newactive, newactive_r, newactive_a);
        }
    }
    nactive = newactive;
    for(Itmp1_ = 0; Itmp1_ <= 10000; ++Itmp1_) {
        active_r[Itmp1_] = newactive_r[Itmp1_];
        active_a[Itmp1_] = newactive_a[Itmp1_];
    }
}

void create(int r, int a, int *newactive, int newactive_r[], int newactive_a[])
// 星のクラスターの生成
{
    if(a < 1)
        a += 6*r;
    if(rnd(0) < p && cell[r][a] != 15)
    {
        ++(*newactive);
        newactive_a[(*newactive)] = a;
        newactive_r[(*newactive)] = r;
        cell[r][a] = 15;                    // 活性なセル
    }
}

void plot_spiral()
{
    int a, r;
    double plotsize, theta, x, y;
    char Stmp1_[_LBUFF_];

    GWsetogn(0, -2);
    sprintf(Stmp1_, "number of active star clusters = %d", nactive);
    GWputtxt((float)(-nring*GWaspect(1)), (float)(nring), Stmp1_);
    for(r = 1; r <= nring; ++r)
    {
        for(a = 1; a <= 6*r; ++a)
        {
            if(cell[r][a] > 0)
            {
                theta = (double)(a*s + v*t)/r;
                x = r*cos(theta);
                y = r*sin(theta);
                plotsize = cell[r][a]/15.0;
                GWputcmb(ICMB, x, y, plotsize, plotsize, 0);
                // 星のクラスターの寿命を減少させる
                --(cell[r][a]);
            }
        }
    }
    GWflush(-2);
}