/**************************************************************************
*                   FRACTALS:   HENON'S  ATTRACTORS                       *
* ----------------------------------------------------------------------- *
* TUTORIAL:                                                               *
*                                                                         *
* The french politician and mathematician Poincaré highlighted the com-   *
* plexity of planet trajectories that simple laws of celestial mechanics  *
* do not show at first. A typical example is the movement of the Earth    *
* around the sun, everybody knows that the period is one year long. This  *
* is only true in first approximation, but not if you take into conside-  *
* ration all the small perturbations caused by the other planets !        *
*                                                                         *
* Periodic trajectories are useful for at least two reasons:              *
*                                                                         *
*  1) they can be predicted with accuracy by Kepler's laws,               *
*                                                                         *
*  2) it is possible to describe conveniently the small trajectory        *
*     variations by using the following:                                  *
*                                                                         *
* Let be T the basic periodic trajectory that we intercept by a fixed,    *
* perpendicular plane, and O the point, taken as origin, of the plane     *
* intercepted by the T trajectory. At next periods, the trajectory with   *
* its small variations will meet the same plane at nearby points A0, A1,  *
* A2 etc.                                                                 *
*                                                                         *
* This allows us to replace a complex problem of 3D trajectories by a     *
* simpler problem of 2D representation of points.                         *
*                                                                         *
* The french mathematician Hénon described in 1969 a dynamic system which *
* illustrates the proposed method:                                        *
*                                                                         *
* 	  x n+1 = x n cos a - ( y n - x n² ) sin a                        *
*	  y n+1 = x n sin a + ( y n - x n² ) cos a                        *
*                                                                         *
* a, being the angle of basic trajectory with respect to the observation  *
* plane (in radians).                                                     *
*                                                                         *
* The program randomly chooses an angle a (usually between 0.5 and 2.0    *
* radians). The number of iterations equals arbitrarily 200.              * 
*                                                                         *
* As one can see, caracreristic figures are obtained, the points have a   *
* noticeable tendancy to agglutinate themselves in particular "zones of   *
* attraction".                                                            *
* ----------------------------------------------------------------------- *
* REFERENCE:                                                              *
*     "Graphisme dans le plan et dans l'espace avec Turbo Pascal 4.0      *
*      By R. Dony - MASSON, Paris 1990" [BIBLI 12].                       *
*                                                                         *
*                Microsoft Visual C++ version in API style by J-P Moreau  *
*                      (Program to use with Henon.mak and Gr2D.cpp)       *
*                                   (www.jpmoreau.fr)                     *
**************************************************************************/
#include <windows.h>
#include <stdio.h>
#include <stdlib.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-5; rct.bottom=iy+15;
  DrawText(hdc, s, -1, &rct,
                DT_SINGLELINE | DT_CENTER | DT_VCENTER);
}

//global variables of Henon program
double angle,cosinus,sinus;
int    nbreiter;

 void Data(void)
 {
   angle=1.854;  //can be changed by user
   cosinus=cos(angle); sinus=sin(angle);
   nbreiter=500;
 }

 //start Henon trajectory at point (x,y)
 void Henon(double x, double y) {
   double aux; int i;
   for (i=0; i<nbreiter; i++) {
     if (x<1.1 && y<0.8) {
       MoveXY(x,y);
       LineXY(x+0.005,y);
     }
     aux=x;
     if (x < 1E12) {
       x=x*cosinus-(y-x*x)*sinus;
       y=aux*sinus+(y-aux*aux)*cosinus;
     }
   }
 }

// Draw Henon's attractor diagram (angle=1.85 rad)
void Draw_Henon() {
  int count;
  double xpos,ypos;
  char ch[15];
  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.1,1.1,-0.8,0.8);
  //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.2,0.2);
  //return to current (blue) pen
  SelectObject(hdc,oldpen);
  //draw a frame around drawing zone
  Bordure();
  //graduate axes
  Gradue(0.2,0.2);
  //initialize
  Data();
  count=0;
  //main Henon loop (start Henon at random point (xpos,ypos))
  do {
    xpos=-1.1+rand()*2.2/RAND_MAX;
    ypos=-0.8+rand()*1.6/RAND_MAX;
    Henon(xpos,ypos);
    count++;
  } while (count<=200);
  //print angle value in bottom right corner
  sprintf(ch,"Angle = %5.2f",angle);
  OutText(MaxX-150,MaxY-60,ch);
} //Draw_Henon()


//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);

            Draw_Henon();
               
            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[] = "Henon";
    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
        " HENON  ATTRACTOR",        // 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 henon.cpp