// PROGRAM bifurcate
// いろいろな r について x の値をプロットする

#include "TrueBASIC.h"
#define f(x,r) (4*(r)*(x)*(1 - (x)))

int main();
void initial(double *x0, double *r0, double *rmax, int *nvalues, double *dr, 
             int *ntransient, int *nplot);
void output(double *x, double r, int ntransient, int nplot, float r0);

int main()
{
    int ir, nplot, ntransient, nvalues;
    double dr, r, rmax, x;
    float r0;

    GWopen(0);

    initial(&x, &r, &rmax, &nvalues, &dr, &ntransient, &nplot);
    r0 = (float)r;
    for(ir = 0; ir <= nvalues; ++ir)
    {
        output(&x, r, ntransient, nplot, r0);
        r += dr;
    }
    // r > 1 を避けるために r の最大値については独立に処理
    output(&x, rmax, ntransient, nplot, r0);

    GWquit();
    return 0;
}

void initial(double *x0, double *r0, double *rmax, int *nvalues, double *dr, 
             int *ntransient, int *nplot)
{
    double mx, xmax;
    char Stmp1_[_LBUFF_];

    printf("initial value of control parameter r = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", r0);
    // r が 1 より大きくならないことが大切
    printf("maximum value of r = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", rmax);
    // dr <= 0.01 としてみよ
    printf("incremental change of r = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", dr);
    printf("number of iterations not plotted = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", ntransient);
    printf("number of iterations plotted = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", nplot);
    (*nvalues) = (int)(((*rmax) - (*r0))/(*dr));   // プロットされる r の値の数
    (*x0) = 0.5;                            // 初期値
    xmax = 1;                               // x の最大値
    mx = 0.05*xmax;                         // 余白
    GWsetogn(0, 1);
    GWindow((float)((*r0)-(*dr)), (float)(-mx), (float)((*rmax)+(*dr)), 
            (float)(xmax + mx));
    GWrect((float)((*r0)), 0.0f, (float)((*rmax)), 1.0f);
}

void output(double *x, double r, int ntransient, int nplot, float r0)
{
    int i;
    char Stmp1_[_LBUFF_];

    GWsetogn(0, -2);
    sprintf(Stmp1_, "r = %f", r);
    GWputtxt(r0, 1.0f, Stmp1_);
    GWflush(-2);                         // 前の出力を消去 
    GWsetogn(0, 1);
    for(i = 1; i <= ntransient; ++i)     // プロットされない x の値
        (*x) = f((*x), r);
    for(i = 1; i <= nplot/2; ++i)
    {
        (*x) = f((*x), r);
        // 与えられた r について種々の x の値をプロットする
        GWsetpxl((float)(r), (float)((*x)), RED);
    }
    // x の値が収束したかどうかが分かるように色を変える
    for(i = (nplot/2 + 1); i <= nplot; ++i)
    {
        (*x) = f((*x), r);
        GWsetpxl((float)(r), (float)((*x)), BLUE);
    }
}