/**************************************************************************
*                   FRACTALS:  FEIGENBAUM  DIAGRAM                        *
* ----------------------------------------------------------------------- *
* TUTORIAL:                                                               *
*                                                                         *
* The dynamic process that allows to very simply simulate "chaos" has the *
* following mathematical form:                                            *
*	  x    = f (x , C ).                                              *
*          n+1       n                                                    *
* The only condition is that the relation between input x n and output    *
* x n+1 be non linear. If this iterative process starts with an arbitrary * 
* value x0, it will give a series of values:  x1, x2,...,x n... The long  *
* term behaviour of this series is interesting.                           *
*                                                                         *  
* Let us consider the classical example of the growth of a population over*
* several years. Let be an initial population x0 which after n years has  *
* become x n. The growth ratio is:                                        *
*                 r = ( x n+1 - x n ) / x n.                              *
* If this ratio remains constant, the dynamic law is:                     *
*                                                                         *
*                 x n+1 = ( 1 + r )  x n                                  *
*                                                                         *
* After n years, the population will be:    x n = ( 1 + r )  x0.          *
*                                                                         *
* To clarify the situation, let us consider the case where x0 = 0,001.    *
* If the growth ratio equals 5%, the population will roughly double       *
* every 15 years. As a matter of fact:                                    *
*                                                                         *
*   x0 = 0,001	       x1 = 0,00105	 x2 = 0,0011025                   *
*   x3 = 0,00115763    x4 = 0,00121551	 x5 = 0,00127628                  *
*   x6 = 0,00134010    x7 = 0,00140710	 x8 = 0,00147746                  *
*   x9 = 0,00155133    x10 = 0,00162889	 x11 = 0,00171034                 *
*   x12 = 0,00179586   x13 = 0,00188565	 x14 = 0,00197993                 *
*   x15 = 0,00207893   ...                                                *
*                                                                         *
* This kind of growth is exponential. But this dynamic law is linear,     *
* which is not judicious. Actually, the real growth of a poulation is     *
* necessarily limited by the ressources of its habitat, which are not in  *
* infinite quantity, such as food, for example. The belgian Verhulst was  *
* one of the first to study this phenomenon in 1845.                      *
* He advised to consider a variable growth ratio r, taking the form       *
* r = r - C x n.  The dynamic law of growth then takes the form:          * 
*                                                                         *
*		  x n+1 = ( 1 + r ) x n - C x² n                          * 
*                                                                         *
* By having C = r / X, the population increases up to the value X, then   *
* stabilizes at this value. At least, this remains true so long as the    *
* growth ratio is < 200 %. A human population has never such a high growth*  
* ratio, but in the case of insects, for example, this can be quite       *
* possible. For a growth ratio even higher, one can observe surprising    *
* results (see verhulst.pas program).                                     *
*                                                                         *
* The calculation begins at x0 = 0,1 X.                                   *
*                                                                         *
* Case r = 1.8                                                            *
*                                                                         *
* The response curve climbs up rapidly and after some oscillations reaches*
* a limit that is an "attractor".                                         *
*                                                                         *
* Case r = 2.3                                                            *
*                                                                         *
* The curve oscilates rapidly between two levels that frame the value X.  *
* The suite has two attractors.                                           *  
*                                                                         *
* Case r = 2.5                                                            *
*                                                                         *
* The suite has four attractors.                                          *
*                                                                         *
* Case r = 3.0                                                            *
*                                                                         *
* The numbers x jump from one value to another one, without having twice  *
* the same value. Such a behaviour can be qualified as "chaotic".         *
*                                                                         *
* The Feigenbaum diagram:                                                 *
*                                                                         *
* To better observe the behaviour of such suites x when r varies, we now  *
* only consider the "attractors" for each r value. The first 100 transient*
* values are skipped then at each r value 300 points are displayed.       *
* For r < 2,57, the behaviour is non chaotic: the attractors are in       *
* limited number. When r > 2,57, the attractors become queer and the dia- *
* gram has more and more ramifications until being fully inextricable:    *
* now we have a chaotic bahaviour!                                        *
*                                                                         *
*  The obtained picture is called the bifurcation diagram or Feigenbaum   *
* tree. By an accurate analysis of the bifurcation points, the mathema-   *
* tician Feigenbaum discovered a new universal constant. The length of the*
* r intervals for which a stable period is obtained,  is shortened, when  *
* the period is doubled, by a factor that tends toward the universal      *
* constant k = 4,669201660910...                                          *
*                                                                         *
* This constant, called the Feigenbaum Constant, can be found in other    *
* chaotic phenomenons, such as fluidic turbulences, chemical reactions,   *
* and even in human heart!                                                *
* ----------------------------------------------------------------------- *
* From "Graphisme dans le plan et dans l'espace avec Turbo Pascal 4.0     *
* By R. Dony - MASSON, Paris 1990, pages 189-192 " [BIBLI 12].            *
*                                                                         *
*               Microsoft C++ Version in API WINDOWS Style By J-P Moreau  *
*                    (Program to use with Chaos.mak and Gr2D.cpp)         *
*                              (www.jpmoreau.fr)                          *
**************************************************************************/
#include <windows.h>
#include <stdio.h>
#include <math.h>

HDC  hdc;
RECT rect;

//Functions of module Gr2D.cpp
void Fenetre(double,double,double,double); 
void Cloture(int,int,int,int); 
void Bordure();
void Axes();
void Gradue(double,double);
void Grille(double,double);
void MoveXY(double,double);
void LineXY(double,double);

//"home made" graphic commands for hdc environment used by above functions
void DrawPixel(int ix,int iy) {
  //sorry, no other available command found
  Rectangle(hdc,rect.left+ix,rect.top+iy,
	            rect.left+ix+2,rect.top+iy+1);
}

void Swap(int *i1,int *i2) {
  int it;
  it=*i1; *i1=*i2; *i2=it;
}

void DrawLine(int ix1,int iy1,int ix2,int iy2) {
  int i,il,ix,iy;
  float dx,dy;
  if (ix2<ix1) {
	Swap(&ix1,&ix2); Swap(&iy1,&iy2);
  }
  il=(int) sqrt((ix2-ix1)*(ix2-ix1)+(iy2-iy1)*(iy2-iy1));
  dx=(float) (ix2-ix1)/il; 
  dy=(float) (iy2-iy1)/il;
  ix=ix1; iy=iy1;
  for (i=1;i<il+1; i++) {
	DrawPixel(ix,iy);
	ix=(int) (ix1+i*dx); 
	iy=(int) (iy1+i*dy);
  }
}

void OutText(int ix,int iy,char *s) {
  RECT rct;
  rct.left=ix; rct.right=ix+10*strlen(s);
  rct.top=iy; rct.bottom=iy+10;
  DrawText(hdc, s, -1, &rct,
                DT_SINGLELINE | DT_CENTER | DT_VCENTER);
}

//draw Feigenbaum diagram
void Feigenbaum() {
  double  dr,r,x;
  int    i;
  extern int MaxX,MaxY;

  HPEN redpen,oldpen;
  //define red pen for grid
  redpen = CreatePen(PS_DOT,1,RGB(255,0,0));
  //client drawing zone in pixels
  MaxX=(rect.right-rect.left);
  MaxY=(rect.bottom-rect.top);
  //define window in physical coordinates
  Fenetre(1.9,3.0,0.0,1.4);
  //drawing zone for curve 2D in pixels
  Cloture(50,MaxX-25,40,MaxY-10);
  //draw ox and oy axes
  Axes();
  //select red pen
  oldpen = SelectObject(hdc,redpen);
  //draw grid
  Grille(0.1,0.1);
  //return to current (blue) pen
  SelectObject(hdc,oldpen);
  //draw a frame around drawing zone
  Bordure();
  //graduate axes
  Gradue(0.2,0.2);
  //main Feigenbaum loop
  r=1.9; dr=0.0025; 
  do {
    x=0.3;
    for (i=0; i<200; i++)  x=(1+r)*x-r*x*x;
    for (i=0; i<300; i++) {
      x=(1+r)*x-r*x*x;
      if (r>=1.95)  {
	MoveXY(r,x);
	LineXY(r+dr,x);
      }
    }
    r+=dr;
  } while (r <= 3);
  //print graph title
  OutText(60,20,"FEIGENBAUM DIAGRAM");
}


//Handle window Paint and Destroy messages
long FAR PASCAL /*_export*/ WndProc(HWND hwnd,
    UINT message, UINT wParam, LONG lParam) 
{
    HPEN        hpen, hpenOld;
    PAINTSTRUCT ps;

    switch (message) {
        case WM_PAINT:
	    hdc = BeginPaint(hwnd, &ps);
	    GetClientRect(hwnd, &rect);
            hpen = CreatePen(PS_SOLID, 1, RGB(0, 0, 255));
            hpenOld = SelectObject(hdc,hpen);
               
            Feigenbaum();
               
            SelectObject(hdc,hpenOld);
            DeleteObject(hpen);
            EndPaint(hwnd, &ps);
            return 0;
        case WM_DESTROY:
            PostQuitMessage(0);
            return 0;
    }
    return DefWindowProc(hwnd, message, wParam, lParam);
}

// main program
int PASCAL WinMain(HANDLE hInstance, HANDLE hPrevInstance,
    LPSTR lpszCmdParam, int nCmdShow)
{
    static char szAppName[] = "Feigenbaum";
    HWND        hwnd;
    MSG         msg;
    WNDCLASS    wndclass;

    if (!hPrevInstance) {
        wndclass.style         = CS_HREDRAW | CS_VREDRAW;
        wndclass.lpfnWndProc   = WndProc;
        wndclass.cbClsExtra    = 0;
        wndclass.cbWndExtra    = 0;
        wndclass.hInstance     = hInstance;
        wndclass.hIcon         = LoadIcon(NULL, IDI_APPLICATION);
        wndclass.hCursor       = LoadCursor(NULL, IDC_ARROW);
        wndclass.hbrBackground = GetStockObject(WHITE_BRUSH);
        wndclass.lpszMenuName  = NULL;
        wndclass.lpszClassName = szAppName;

        RegisterClass(&wndclass);
    }

    hwnd = CreateWindow(szAppName,  // window class name
        "FEIGENBAUM",               // window caption
        WS_OVERLAPPEDWINDOW,        // window style                       
        CW_USEDEFAULT,              // initial x position
        CW_USEDEFAULT,              // initial y position
        CW_USEDEFAULT,              // initial x size
        CW_USEDEFAULT,              // initial y size
        NULL,                       // parent window handle
        NULL,                       // window menu handle
        hInstance,                  // program instance handle
        NULL);                      // creation parameters

   ShowWindow(hwnd, nCmdShow);
   UpdateWindow(hwnd);

   while (GetMessage(&msg, NULL, 0, 0)) {
       TranslateMessage(&msg);
       DispatchMessage(&msg);
   }
   return msg.wParam;
}

//end of file chaos.cpp