home *** CD-ROM | disk | FTP | other *** search
/ Rat's Nest 1 / ratsnest1.iso / prgmming / c / atractor.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1996-02-19  |  8.8 KB  |  281 lines
  1. //
  2.  
  3. //+---------------------------------------------------------------------+
  4.  
  5. //+ Program ATRACTOR.CPP  16/05/94                                      +
  6.  
  7. //+ By Rafael GOMEZ {gomez@gosc.enst-bretagne.fr} (France)              +
  8.  
  9. //+                                                                     +
  10.  
  11. //+ PLOTS A ROTATING Lorenz ATRACTOR IN 3D.                             +
  12.  
  13. //+                                                                     +
  14.  
  15. //+ Based on the program LORENZ.CPP                                     +
  16.  
  17. //+ By Ramiro Perez {RPEREZ@UTPVM1.BITNET}, (Panama)                    +
  18.  
  19. //+ and Fausto A. A. Barbuto {BJ06@C53000.PETROBRAS.ANRJ.BR}, (Brazil). +
  20.  
  21. //+ C++ 3.1 programme's creator: Fausto A. A. Barbuto, April 14, 1994.  +
  22.  
  23. //+ After a TURBO BASIC program by Ramiro Perez.                        +
  24.  
  25. //+                                                                     +
  26.  
  27. //+ You can rotate it with the numeric keyboard (NumLock = ON)          +
  28.  
  29. //+ Keys:        (NumLock = ON)                                         +
  30.  
  31. //+      4 and 6 rotate the atractor around the Y axis (azimuth).       +
  32.  
  33. //+      2 and 8 change the elevation of the looking point              +
  34.  
  35. //+              (Rotation around the screen's X axis)                  +
  36.  
  37. //+         5    returns the azimuth and elevation to 0.0               +
  38.  
  39. //+      - and + decrease or increase the step angle of the rotations.  +
  40.  
  41. //+      / and * decrease or increase the number of iterations.         +
  42.  
  43. //+      A or  a turns on/of the reference axis                         +
  44.  
  45. //+       SPACE  turns on/off the repetition of the last command        +
  46.  
  47. //+              (with this you can have continuous rotation)           +
  48.  
  49. //+      < and > decrease or increase the integration interval dt       +
  50.  
  51. //+      Q or q  Quit                                                   +
  52.  
  53. //+                                                                     +
  54.  
  55. //+ VGA 16 colours, 640x480 points version.                             +
  56.  
  57. //+---------------------------------------------------------------------+
  58.  
  59. //
  60.  
  61.  
  62.  
  63. #include <graphics.h>
  64.  
  65. #include <conio.h>
  66.  
  67. #include <stdio.h>
  68.  
  69. #include <dos.h>
  70.  
  71. #include <math.h>
  72.  
  73.  
  74.  
  75. #define centx 320
  76.  
  77. #define centy 240
  78.  
  79. #define axesx 35
  80.  
  81. #define axesy 440
  82.  
  83. #define longaxe 30
  84.  
  85. #define pi 3.141592654
  86.  
  87. #define twopi 6.283185307
  88.  
  89. #define halfpi 1.570796327
  90.  
  91.  
  92.  
  93. int get_cmd (int *);
  94.  
  95. double bound (double);
  96.  
  97. void update_param (int, double *, int *, int *, int *);
  98.  
  99. void axe_x ();
  100.  
  101. void axe_y ();
  102.  
  103. void axe_z ();
  104.  
  105. void draw_axes ();
  106.  
  107. void print_info (int, int, double);
  108.  
  109.  
  110.  
  111. static double azim, elev, fxx, fxz, fyx, fyy, fyz;
  112.  
  113.  
  114.  
  115. void main()
  116.  
  117. {
  118.  
  119.         int c, b, a, l, n, kb, cl, step, iter, axis = 1, rep = 0;
  120.  
  121.         double dt, xa[2], ya[2], za[2], x, y, z, x1, y1, z1, xd, yd;
  122.  
  123.         int graphdriver = VGA, graphmode = VGAHI;
  124.  
  125.  
  126.  
  127.         initgraph (&graphdriver,&graphmode,"");
  128.  
  129.         settextstyle (SMALL_FONT, HORIZ_DIR, 4);
  130.  
  131.  
  132.  
  133.         azim = 0.0;
  134.  
  135.         elev = 0.0;
  136.  
  137.         step = 9;
  138.  
  139.         iter = 1000;
  140.  
  141.         dt = 0.03;
  142.  
  143.         a  = 5;
  144.  
  145.         b  = 15;
  146.  
  147.         c  = 1;
  148.  
  149.  
  150.  
  151.         do {
  152.  
  153.                 cleardevice ();
  154.  
  155. //
  156.  
  157. //  Start-up values for xa, ya and za.
  158.  
  159. //
  160.  
  161.                 xa[0] = 3.051522;
  162.  
  163.                 ya[0] = 1.592542;
  164.  
  165.                 za[0] = 15.62388;
  166.  
  167.                 xa[1] = 3.051522;
  168.  
  169.                 ya[1] = 1.582542;
  170.  
  171.                 za[1] = 15.62388;
  172.  
  173.  
  174.  
  175.                 fxx = cos(azim);
  176.  
  177.                 fxz = -sin(azim);
  178.  
  179.                 fyx = sin(elev)*fxz;
  180.  
  181.                 fyy = cos(elev);
  182.  
  183.                 fyz = -sin(elev)*fxx;
  184.  
  185.                 for (n=1; n<=iter; n++) {
  186.  
  187.                         for (l=0;l<=1;l++) {
  188.  
  189.                                 x = xa[l];
  190.  
  191.                                 y = ya[l];
  192.  
  193.                                 z = za[l];
  194.  
  195.                                 x1 = x - a*x*dt + a*y*dt;
  196.  
  197.                                 y1 = y + b*x*dt - y*dt - z*x*dt;
  198.  
  199.                                 z1 = z - c*z*dt + x*y*dt;
  200.  
  201.                                 xa[l] = x1;
  202.  
  203.                                 ya[l] = y1;
  204.  
  205.                                 za[l] = z1;
  206.  
  207.                                 z1 -= 16.0;
  208.  
  209.                                 xd = 19.3*(x1*fxx + z1*fxz) + centx;
  210.  
  211.                                 yd = -11.0*(x1*fyx + y1*fyy + z1*fyz) + centy;
  212.  
  213.                                 putpixel((int)xd, (int)yd, (l==0) ? LIGHTBLUE :
  214.  
  215. LIGHTGREEN);
  216.  
  217.                         }
  218.  
  219.                 }
  220.  
  221.                 if (axis) draw_axes ();
  222.  
  223.                 if (!rep) print_info (step, iter, dt);
  224.  
  225.                 kb = get_cmd (&rep);
  226.  
  227.                 update_param (kb, &dt, &step, &iter, &axis);
  228.  
  229.         } while (kb!='q' && kb!='Q');
  230.  
  231.         closegraph ();
  232.  
  233. }
  234.  
  235.  
  236.  
  237.  
  238.  
  239. int get_cmd (int *rep)
  240.  
  241. {
  242.  
  243.         static int lkb;
  244.  
  245.         int kb;
  246.  
  247.  
  248.  
  249.         if (*rep) {
  250.  
  251.                 if (kbhit ()) kb = getch ();
  252.  
  253.                 if (kb==' ') *rep = 0;
  254.  
  255.                 else kb = lkb;
  256.  
  257.         }
  258.  
  259.         else {
  260.  
  261.                 for ( ; !kbhit (); );
  262.  
  263.                 kb = getch ();
  264.  
  265.                 if (kb==' ') {
  266.  
  267.                         *rep = 1;
  268.  
  269.                         kb = lkb;
  270.  
  271.                 }
  272.  
  273.                 else lkb = kb;
  274.  
  275.         }
  276.  
  277.         return (kb);
  278.  
  279. }
  280.  
  281.  
  282.  
  283.  
  284.  
  285. void axe_x ()
  286.  
  287. {
  288.  
  289.         int br, xd, yd;
  290.  
  291.         br = azim>=0.0 && azim<=pi || elev==halfpi || elev==3*halfpi;
  292.  
  293.         if (elev>halfpi && elev<3*halfpi) br = !br;
  294.  
  295.         setcolor (br ? CYAN : BLUE);
  296.  
  297.         xd = longaxe*fxx + axesx;
  298.  
  299.         yd = -longaxe*fyx + axesy;
  300.  
  301.         line (axesx, axesy, xd, yd);
  302.  
  303.         outtextxy (xd, yd, "X");
  304.  
  305.         return;
  306.  
  307. }
  308.  
  309.  
  310.  
  311. void axe_y ()
  312.  
  313. {
  314.  
  315.         int xd, yd;
  316.  
  317.         setcolor ((elev<=pi) ? YELLOW : GREEN);
  318.  
  319.         xd = axesx;
  320.  
  321.         yd = -longaxe*fyy + axesy;
  322.  
  323.         line (axesx, axesy, xd, yd);
  324.  
  325.         outtextxy (xd + 2, yd, "Y");
  326.  
  327.         return;
  328.  
  329. }
  330.  
  331.  
  332.  
  333. void axe_z ()
  334.  
  335. {
  336.  
  337.         int br, xd, yd;
  338.  
  339.         br = azim>=(3*halfpi) || azim<=halfpi || elev==halfpi || elev==3*halfpi;
  340.  
  341.         if (elev>halfpi && elev<3*halfpi) br = !br;
  342.  
  343.         setcolor (br ? LIGHTRED : RED);
  344.  
  345.         xd = longaxe*fxz + axesx;
  346.  
  347.         yd = -longaxe*fyz + axesy;
  348.  
  349.         line (axesx, axesy, xd, yd);
  350.  
  351.         outtextxy (xd, yd, "Z");
  352.  
  353.         return;
  354.  
  355. }
  356.  
  357.  
  358.  
  359. void draw_axes ()
  360.  
  361. {
  362.  
  363.         int dr = 1, ord;
  364.  
  365.  
  366.  
  367. //These if's are for line-hiding in the reference axis.
  368.  
  369. //If you want more speed you can delete them.
  370.  
  371.  
  372.  
  373.         if (elev>=pi && elev<=twopi) {
  374.  
  375.                 axe_y ();
  376.  
  377.                 dr = 0;
  378.  
  379.         }
  380.  
  381.         ord = azim>=3*halfpi && azim<=twopi;
  382.  
  383.         if (elev>halfpi) ord = !ord;
  384.  
  385.         if (ord) {
  386.  
  387.                 axe_x ();
  388.  
  389.                 axe_z ();
  390.  
  391.         }
  392.  
  393.         else {
  394.  
  395.                 axe_z ();
  396.  
  397.                 axe_x ();
  398.  
  399.         }
  400.  
  401.         if (dr) axe_y ();
  402.  
  403.         return;
  404.  
  405. }
  406.  
  407.  
  408.  
  409. void print_info (int step, int iter, double dt)
  410.  
  411. {
  412.  
  413.         char mess[100];
  414.  
  415.         setcolor (WHITE);
  416.  
  417.         outtextxy (216, 5, "L O R E N Z    A T R A C T O R");
  418.  
  419.         sprintf (mess, "Elev= %4.0f deg    Azim= %4.0f deg    Step= %i deg
  420.  
  421.                         Iterations= %i    dt= %4.3f", elev*180/pi, azim*180/pi,
  422.  
  423. step, iter, dt);
  424.  
  425.         outtextxy (80, 460, mess);
  426.  
  427.         return;
  428.  
  429. }
  430.  
  431.  
  432.  
  433.  
  434.  
  435. double bound (double angle)
  436.  
  437. {
  438.  
  439.         double angb;
  440.  
  441.  
  442.  
  443.         if (angle>=twopi) angb = angle - twopi;
  444.  
  445.         else if (angle<0.0) angb = angle + twopi;
  446.  
  447.         else angb = angle;
  448.  
  449.         return (angb);
  450.  
  451. }
  452.  
  453.  
  454.  
  455.  
  456.  
  457. void update_param (int kb, double *dt, int *step, int *iter, int *axis)
  458.  
  459. {
  460.  
  461.         double delta;
  462.  
  463.  
  464.  
  465.         delta = (*step)*pi/180.0;
  466.  
  467.  
  468.  
  469.         switch (kb) {
  470.  
  471.         case('4'):
  472.  
  473.                 azim += delta;
  474.  
  475.                 break;
  476.  
  477.         case('6'):
  478.  
  479.                 azim -= delta;
  480.  
  481.                 break;
  482.  
  483.         case('2'):
  484.  
  485.                 elev += delta;
  486.  
  487.                 break;
  488.  
  489.         case('8'):
  490.  
  491.                 elev -= delta;
  492.  
  493.                 break;
  494.  
  495.         case('5'):
  496.  
  497.                 azim = 0.0;
  498.  
  499.                 elev = 0.0;
  500.  
  501.                 break;
  502.  
  503.         case('*'):
  504.  
  505.                 (*iter) += 500;
  506.  
  507.                 break;
  508.  
  509.         case('/'):
  510.  
  511.                 if (*iter>500) (*iter) -= 500;
  512.  
  513.                 break;
  514.  
  515.         case('+'):
  516.  
  517.                 if (*step<90.0) (*step)++;
  518.  
  519.                 break;
  520.  
  521.         case('-'):
  522.  
  523.                 if (*step>0.0) (*step)--;
  524.  
  525.                 break;
  526.  
  527.         case('>'):
  528.  
  529.                 if (*dt<0.05) (*dt) += 0.005;
  530.  
  531.                 break;
  532.  
  533.         case('<'):
  534.  
  535.                 if (*dt>0.0) (*dt) -= 0.005;
  536.  
  537.                 break;
  538.  
  539.         case('a'):
  540.  
  541.         case('A'):
  542.  
  543.                 *axis = !(*axis);
  544.  
  545.                 break;
  546.  
  547.         default:
  548.  
  549.                 break;
  550.  
  551.         }
  552.  
  553.         azim = bound (azim);
  554.  
  555.         elev = bound (elev);
  556.  
  557.         return;
  558.  
  559. }
  560.  
  561.