// PROGRAM bifurcate
// plot values of x for different values of r

#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;
    }
    // maximum value of r done separately to avoid r > 1
    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);
    // important that r not be greater than 1
    printf("maximum value of r = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", rmax);
    // suggest 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));   // number of r values plotted
    (*x0) = 0.5;                            // initial value
    xmax = 1;                               // maximum value of x
    mx = 0.05*xmax;                         // margin
    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);                         // erase previous output 
    GWsetogn(0, 1);
    for(i = 1; i <= ntransient; ++i)     // x values not plotted
        (*x) = f((*x), r);
    for(i = 1; i <= nplot/2; ++i)
    {
        (*x) = f((*x), r);
        // show different x-values for given value of r
        GWsetpxl((float)(r), (float)((*x)), RED);
    }
    // change color to see if values of x have converged
    for(i = (nplot/2 + 1); i <= nplot; ++i)
    {
        (*x) = f((*x), r);
        GWsetpxl((float)(r), (float)((*x)), BLUE);
    }
}