home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / src / X11 / endo / comp.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-05-12  |  21.5 KB  |  821 lines
  1. /*************************************************************************
  2.  *                                                                       *
  3.  *  Copyright (c) 1992, 1993 Ronald Joe Record                           *
  4.  *                                                                       *
  5.  *  All rights reserved. No part of this program or publication may be   *
  6.  *  reproduced, transmitted, transcribed, stored in a retrieval system,  *
  7.  *  or translated into any language or computer language, in any form or *
  8.  *  by any means, electronic, mechanical, magnetic, optical, chemical,   *
  9.  *  biological, or otherwise, without the prior written permission of:   *
  10.  *                                                                       *
  11.  *      Ronald Joe Record (408) 458-3718                                 *
  12.  *      212 Owen St., Santa Cruz, California 95062 USA                   *
  13.  *                                                                       *
  14.  *************************************************************************/
  15.  
  16. #include "endo.h"
  17.  
  18. int maxcrit = 4, numattrs = 0, use_avg_lyap=1;
  19. int begin, first, middle, second;
  20. double lastx, lasty;
  21. extern int idown;
  22.  
  23. void
  24. initial_kk(p)
  25. double *p;
  26. {
  27.     static double t;
  28.     extern double sqrt();
  29.  
  30.     if (p[4] == 0.0)
  31.         return;
  32.     t = ((p[1]*p[3])-(p[0]*p[2]))/(3.0*p[4]);
  33.     if (t > 0.0) {
  34.         start_x = sqrt(t);
  35.         if (p[1] > 0.0)
  36.             start_y = -start_x; /* initial condition on the off-diagonal */
  37.         else
  38.             start_y = start_x; /* initial condition on the diagonal */
  39.     }
  40. }
  41.  
  42. void
  43. initial_q2(p)
  44. double *p;
  45. {
  46.     if ((p[2] == 0.0) || (p[1] == 0.0))
  47.         return;
  48. /*    To find a critical point on the Y-axis, use this */
  49. /*
  50.  *    start_x = 0.0;
  51.  *    start_y = ((p[0]*p[5]/p[2]) - p[3])/p[1];
  52.  */
  53. /*    To find a critical point on the X-axis, use this */
  54. /*
  55.  *    start_y = 0.0;
  56.  *    start_x = ((p[0]*p[5]) - (p[2]*p[3]))/((2.0*p[2]*p[4])-(p[0]*p[1]));
  57.  */
  58. /*    To find a critical point with x=1, use this */
  59. /*
  60.  *    start_x = 1.0;
  61.  *    start_y = ((p[0]*(p[5]+p[1])/p[2]) - p[3] - (2.0*p[4]))/p[1];
  62.  */
  63. /*    To find a critical point with x=-1, use this */
  64. /*
  65.  *    start_x = -1.0;
  66.  *    start_y = ((p[0]*(p[5]-p[1])/p[2]) - p[3] + (2.0*p[4]))/p[1];
  67.  */
  68. /*    To find a critical point on the diagonal, use this */
  69.      start_y = (p[1]*p[2])+(2.0*p[2]*p[4])-(p[0]*p[1]);
  70.     if (start_y == 0.0)
  71.         start_x = 0.0;
  72.     else {
  73.          start_x = ((p[0]*p[5]) - (p[2]*p[3]))/start_y;
  74.         start_y = start_x;
  75.     }
  76. }
  77.  
  78. void
  79. initial_alex(p)
  80. double *p;
  81. {
  82.     static double t;
  83.     extern double sqrt();
  84.  
  85.     while (((t = 
  86.     (p[0]*p[1])+1.0-(2.0*(p[0]-p[1])*start_y)-(4.0*start_y*start_y))<0.0) &&
  87.     (start_y > 1.0e-12))
  88.         start_y *= 0.5;
  89.     if (t >= 0.0)
  90.         start_x = 0.5 * sqrt(t);
  91. }
  92.  
  93. void
  94. initial_secant(p)
  95. double *p;
  96. {
  97.     extern double poly();
  98.  
  99.     start_x = poly(start_y, p);
  100. }
  101.  
  102. void
  103. initial_julia(p)
  104. double *p;
  105. {
  106.     extern double sqrt();
  107.  
  108.     start_x = start_y = 0.0;
  109.     if ((p[3] >= 0.0) || (p[2] == 0.0))
  110.         return;
  111.     start_y = (p[2]/2.0) * sqrt(-p[3]);
  112. }
  113.  
  114. void
  115. initial_brussel(p)
  116. double *p;
  117. {
  118.     static double a, b;
  119.     extern double sqrt();
  120.  
  121.     if ((p[0] == 0.0) || ((p[0] + 1.0) == 0.0))
  122.         return;
  123.     a = (p[1] + 1.0) - (1.0 / p[0]);
  124.     b = a / (p[0] + 1.0);
  125.     if (b > 0.0)
  126.         start_x = start_y = ABS(sqrt(b));
  127. }
  128.  
  129. void
  130. initial_gucken(p)
  131. double *p;
  132. {
  133.     if ((p[2] == 0.0) || ((p[0] + p[1]) == 0.0))
  134.         return;
  135.     start_x = start_y = p[1] / (p[2] * (p[0] + p[1]));
  136. }
  137.  
  138. /* Pick a point on the critical arc as the initial condition */
  139. void
  140. initial_critical(p)
  141. double *p;
  142. {
  143.     static int i, j, done;
  144.     static pair c;
  145.  
  146.     done = 0;
  147.     for (i=0; i<250; i++) {
  148.         for (j=0; j<250; j++) {
  149.             if (ABS(((*deriv)((double)i, (double)j, p)).x) < cdelt) {
  150.                 start_x = (double)i; start_y = (double)j;
  151.                 c = (*map)((double)i, (double)j, p);
  152.                 if (ABS(((*deriv)(c.x, c.y, p)).x) < cdelt) {
  153.                     done = 1;
  154.                     break;
  155.                 }
  156.             }
  157.         }
  158.         if (done) break;
  159.     }
  160. }
  161.  
  162. /* Pick a point on the critical curve where x = y/2 */
  163. void
  164. initial_logistic(p)
  165. double *p;
  166. {
  167.     static double a, b, c, q;
  168.     extern double sqrt();
  169.  
  170.     if ((p[0] == 0.0) || (p[1] == 0.0))
  171.         return;
  172.     a = 128.0 * p[0] * p[1];
  173.     b = -96.0 * p[0] * p[1];
  174.     c = (15.0 * p[0] * p[1]) + p[0] + p[1] - 1.0;
  175.     if (b < 0)
  176.         q = -0.5 * (b - (sqrt((b*b) - (4.0*a*c))));
  177.     else
  178.         q = -0.5 * (b + (sqrt((b*b) - (4.0*a*c))));
  179.     if ((start_x = Max(q/a, c/q)) > 0.5) {
  180.         if ((start_x = Min(q/a, c/q)) <= 0.0) {
  181.             start_x = 0.3; 
  182.             start_y = 0.6;
  183.         }
  184.         else {
  185.             if (start_x <= 0.5)
  186.                 start_y = 2.0 * start_x;
  187.             else {
  188.                 start_x = 0.3; 
  189.                 start_y = 0.6;
  190.             }
  191.         }
  192.     }
  193.     else {
  194.         if (start_x > 0.0)
  195.             start_y = 2.0 * start_x;
  196.         else {
  197.             start_x = 0.3; 
  198.             start_y = 0.6;
  199.         }
  200.     }
  201. }
  202.  
  203. int
  204. preiterate(c, i, x1, y1)
  205. int c, i;
  206. double x1, y1;
  207. {
  208.     static int n;
  209.     extern void draw_precrit();
  210.  
  211.     if (i > maxcrit)
  212.         return(0);
  213.     if (precrit == 1) {
  214.         if (ABS(((*deriv)(x1, y1, params)).x) < cdelt) {
  215.             draw_precrit(x,y,(i%first)+begin);
  216.             return(0);
  217.         }
  218.       }
  219.     else if (precrit == 2) {/* compute pre-iterates of the diagonal */
  220.         if (ABS(x1 - y1) < cdelt) {
  221.             draw_precrit(x,y,(i%first)+begin);
  222.             return(0);
  223.         }
  224.       }
  225.       else if (precrit == 4) {/* compute pre-iterates of the origin */
  226.         if ((ABS(x1) < cdelt) && (ABS(y1) < cdelt)) {
  227.             draw_precrit(x,y,(i%second)+middle);
  228.             return(0);
  229.         }
  230.       }
  231.       else if (precrit == 5) {/* compute pre-iterates of the origin */
  232.         n = 1;
  233.         if (c < 2) {/* and diagonal, color origin differently */
  234.             if (ABS(x1 - y1) < cdelt) {
  235.                 draw_precrit(x,y,(i%second)+middle);
  236.                 n = 2;
  237.             }
  238.         }
  239.         if ((ABS(x1) < cdelt) && (ABS(y1) < cdelt)) {
  240.             draw_precrit(x,y,(i%first)+begin);
  241.             n = 0;
  242.         }
  243.         return(n);
  244.     }
  245.       else if (precrit == 6) {/* compute pre-iterates of the origin */
  246.         n = 1;
  247.         if (c < 2) {/* and critical, color origin differently */
  248.             if (ABS(((*deriv)(x1, y1, params)).x) < cdelt) {
  249.                 draw_precrit(x,y,(i%second)+middle);
  250.                 n = 2;
  251.             }
  252.         }
  253.         if ((ABS(x1) < delta) && (ABS(y1) < delta)) {
  254.             draw_precrit(x,y,(i%first)+begin);
  255.             n = 0;
  256.         }
  257.         return(n);
  258.       }
  259. }
  260.  
  261. double
  262. set_lyapunov(a, b)
  263. double a, b;
  264. {
  265.     double s;
  266.     extern double sqrt();
  267.  
  268.     if (use_avg_lyap)
  269.         return(a);
  270.     else {
  271.         s = (b*b) - (4.0*a);
  272.         if (s < 0.0) /* 4a > b^2 => a > 0 so return the positive average */
  273.             return(a);
  274.         s = sqrt(s);
  275.         return((b + s)/2.0); /* the smaller exponent is (b - s)/2 */
  276.     }
  277. }
  278.  
  279. void
  280. calc_lyapunov(prod,p_total,tr_prod,p_tr_total,p_alpha,p_beta,p_lyapunov)
  281. double tr_prod, prod;
  282. double *p_total;
  283. double *p_tr_total;
  284. double *p_alpha;
  285. double *p_beta;
  286. double *p_lyapunov;
  287. {
  288.     if (prod == 0.0)
  289.         *p_total += LN_MINDOUBLE;
  290.     else if (prod < 0.0)
  291.         *p_total += LN_MAXDOUBLE;
  292.     else
  293.         *p_total += log(prod);
  294.     if (tr_prod == 0.0)
  295.         *p_tr_total += LN_MINDOUBLE;
  296.     else if (tr_prod < 0.0)
  297.         *p_tr_total += LN_MAXDOUBLE;
  298.     else
  299.         *p_tr_total += log(tr_prod);
  300.     *p_alpha = (*p_total * M_LOG2E) / (double)dwell;
  301.     *p_beta = (*p_tr_total * M_LOG2E) / (double)dwell;
  302.     *p_lyapunov = set_lyapunov(*p_alpha, *p_beta);
  303. }
  304.  
  305. compendo()
  306. {
  307.     register i, j;
  308.     static int foundit, crit_hit, index, started=0;
  309.     static int color=MAXCOLOR; 
  310.     static double radius;
  311.     double tr_total, total, tr_prod, prod, alpha, beta, lyapunov; 
  312.     double attrs[MAXATTRS][4];
  313. #ifdef mips
  314.     register double d, x1, y1;
  315. #else
  316.     double d, x1, y1; 
  317. #endif
  318.     pair coordinates;
  319.     int found_crit;
  320.     extern int lowrange, lower, upper;
  321.     extern int xmargin, ymargin, animate, mandel, axes;
  322. #ifdef NorthSouth
  323.     extern double singularity, S_singularity, pB;
  324.     extern void calculate_params();
  325.     extern void solve_quadratic(), gamma_quadratic();
  326. #endif
  327.     extern void save_to_file(),draw_portrait();
  328.     extern void draw_critical(), Updt_Therm();
  329.     extern double drand48(), atan2();
  330.  
  331.     foundit=0;        /* whether we've detected an attractor */
  332.     found_crit=0;        /* whether we've detected a critical point */
  333.     crit_hit = 1;    /* whether a critical pre-iterate has been found */
  334.     if (!run)
  335.         return TRUE;
  336.     if (lyap) {
  337.         if (thermometer)
  338.             Updt_Therm(trajec, pixtra, trawidth, traheight, params[p1], x);
  339.         params[p1] = x; params[p2] = y;
  340.         if ((!Xflag) && (!Yflag)) {
  341. #ifdef NorthSouth
  342.             if (mapindex == 7)  /* North-South model */
  343.                 calculate_params(mapindex, params);
  344. #endif
  345.             /*double logistic*/
  346.             if ((mapindex == (9+MAP_OFF)) || (mapindex == (10+MAP_OFF)))
  347.                 initial_logistic(params);
  348.             if (mapindex == 2) /* Guckenheimer Leslie Matrix */
  349.                 initial_gucken(params);
  350.             if (mapindex == (14+MAP_OFF)) /* Brusselator */
  351.                 initial_julia(params);
  352.             if (mapindex == (15+MAP_OFF)) /* Brusselator */
  353.                 initial_brussel(params);
  354.             if (mapindex == (16+MAP_OFF)) /* Brusselator */
  355.                 initial_secant(params);
  356.             if (mapindex == (17+MAP_OFF)) /* Alexander */
  357.                 initial_alex(params);
  358.             if (mapindex == (22+MAP_OFF)) /* Mira-Agg (Q2) */
  359.                 initial_q2(params);
  360.             if (mapindex == (24+MAP_OFF)) /* Kawakami-Kawasaki (KK) */
  361.                 initial_kk(params);
  362.         }
  363.         if (lyap == 3) {
  364.             if (started) {
  365.                 start_x = lastx; start_y = lasty;
  366.             }
  367.             else
  368.                 started = 1;
  369.         }
  370.         if (lyap == 4)
  371.             initial_critical(params);
  372.         x1 = start_x; y1 = start_y;
  373.     }
  374.     else {
  375. #ifdef NorthSouth
  376.         if (mapindex == 7)
  377.             if ((x == singularity) || (y == S_singularity))/* N-S model #1*/
  378.                 return FALSE;
  379. #endif
  380.         x1 = x; y1 = y;
  381.     }
  382.     if (animate) {
  383.         if (thermometer)
  384.             i = trawidth - THERMWIDTH;
  385.         else
  386.             i = trawidth;
  387.         XCopyArea(dpy,pixtra,trajec,Data_GC[1],xmargin, ymargin,
  388.                 i,traheight,xmargin,ymargin);
  389.         XFillRectangle(dpy, pixtra, Data_GC[0], xmargin, ymargin, 
  390.                 i, traheight);
  391.         if (axes)
  392.             Draw_Axes(trajec, axes);
  393.     }
  394.     if (mapindex == (7+MAP_OFF)) /* Rotor map */
  395.         if (ABS(params[1]) < 1.0)
  396.             return FALSE;
  397.     /* initialize color bands to begin---first-----middle---second
  398.      * based on value of negative where begin is the initial color
  399.      * index, first a range of color indices, middle divides the
  400.      * color wheel, and second is another range of indices.
  401.      */
  402.     if (negative) {
  403.         begin = MINCOLINDEX;
  404.         first = numfreecols-1;
  405.         middle = STARTCOLOR;
  406.         second = lowrange-1;
  407.     }
  408.     else {
  409.         middle = MINCOLINDEX;
  410.         second = numfreecols-1;
  411.         begin = STARTCOLOR;
  412.         first = lowrange-1;
  413.     }
  414.     /*
  415.      * calculate the determinant of the jacobian and color the point
  416.      * if it's on the critical curve (that is, if the determinant
  417.      * is "close" to zero).
  418.      */
  419.     if (critical) {
  420.         d = ABS(((*deriv)(x1, y1, params)).x);
  421.         if (d < cdelt) {
  422.             draw_critical(crijec, x1,y1,begin);
  423.             crit_pts[0][numcrits] = x1;    
  424.             n_crit_pts[0][numcrits] = x1;    
  425.             crit_pts[1][numcrits] = y1;    
  426.             n_crit_pts[1][numcrits++] = y1;    
  427.             found_crit=1;
  428.         }
  429.     }
  430.     /*
  431.      * initialize periods (how many iterates till within delta of periodic)
  432.      * and indices (index into the color map based upon value of periods)
  433.      */
  434.     periods[frame][perind[frame]] = indices[frame][perind[frame]] = 0;
  435.     for (i=0;i<settle;i++) {
  436.         coordinates = (*map)(x1, y1, params);
  437.         x1 = coordinates.x;
  438.         y1 = coordinates.y;
  439.         if ((found_crit) && (critical) && (i < maxcrit)) {
  440.             draw_critical(crijec, x1,y1,(i+1)%first+begin);
  441.             if ((i == 0) && (d = ABS(((*deriv)(x1, y1, params)).x)) < cdelt)
  442.                 set_arc(x1, y1);
  443.         }
  444.         if (crit_hit)
  445.             crit_hit = preiterate(crit_hit, i, x1, y1);
  446.         if (mandel) {    /* we test to see if going to infinity */
  447.             radius = (x1*x1) + (y1*y1);
  448.             if (radius > maxradius) {
  449.                 if (mandel == 4)
  450.                       indices[frame][perind[frame]] = (i % lowrange) + STARTCOLOR;
  451.                 else
  452.                       indices[frame][perind[frame]] = i % 2;
  453.                 periods[frame][perind[frame]] = -i;
  454.                 if (find)
  455.                     basins[frame][perind[frame]] = 0;
  456.                   BufferPoint(dpy,which,pixmap,Data_GC,&Points,
  457.                       indices[frame][perind[frame]++],point.x,height-point.y-1);
  458.                 lastx = x1; lasty = y1;
  459.                 return FALSE;
  460.             }
  461.         }
  462.     }
  463.     if (portrait)
  464.         if (lyap)
  465.             color=(((((params[p1]-min_x)/x_range)+
  466.                 ((params[p2]-min_y)/y_range))/2.0)*second) + middle;
  467.         else
  468.             if (++color >= numcolors)
  469.                 color = begin;
  470.         prod = 1.0; total = 0.0;
  471.         tr_prod = 1.0; tr_total = 0.0;
  472.         for (i=0;i<dwell;i++) {
  473.             coords[frame][0][i] = x1;
  474.             coords[frame][1][i] = y1;
  475.             coordinates = (*map)(x1, y1, params);
  476.             x1 = coordinates.x;
  477.             y1 = coordinates.y;
  478.             if (mandel) {    /* we test to see if going to infinity */
  479.                 radius = (x1*x1) + (y1*y1);
  480.                 if (radius > maxradius) {
  481.                     if (mandel == 4)
  482.                           indices[frame][perind[frame]] = 
  483.                                 ((i+settle) % lowrange) + STARTCOLOR;
  484.                     else
  485.                           indices[frame][perind[frame]]=(i+settle)%2;
  486.                     periods[frame][perind[frame]] = -(i+settle);
  487.                     if (find)
  488.                         basins[frame][perind[frame]] = 0;
  489.                       BufferPoint(dpy,which,pixmap,Data_GC,&Points,
  490.                       indices[frame][perind[frame]++],point.x,height-point.y-1);
  491.                     if (portrait) {
  492.                         draw_portrait(x1,y1,color);
  493.                         FlushBuffer(dpy, trajec, pixtra, Data_GC, &Orbits,
  494.                                         color, color+1);
  495.                     }
  496.                     lastx = x1; lasty = y1;
  497.                     return FALSE;
  498.                 }
  499.             }
  500.         /*
  501.          * calculate the determinant of the jacobian and color the point
  502.          * if one of its iterates falls on the critical curve. In this way,
  503.          * we display the pre-iterates of the critical curve (colored 
  504.          * according to how many iterates it took to get to the critical
  505.          * curve). If precrit is set to 2 (-P 4 argument), then we calculate
  506.         * instead the pre-iterates of the diagonal. If precrit is set to
  507.         * 4 (-P 6 argument), then we calculate the pre-images of the origin.
  508.         * If precrit is set to 5 (-P 7 argument), then we calculate both
  509.         * pre-images of the diagonal and origin. If precrit == 6, then
  510.         * the pre-images of the critical curve and the origin are calculated
  511.          */
  512.         if (crit_hit)
  513.             preiterate(crit_hit, i+settle, x1, y1);
  514.         if ((found_crit) && (critical) && ((settle + i) < maxcrit)) {
  515.           draw_critical(crijec, x1,y1,(i+settle+1)%first+begin);
  516.           if ((settle==0)&&(i==0)&&(d=ABS(((*deriv)(x1,y1,params)).x))<cdelt)
  517.               set_arc(x1, y1);
  518.         }
  519.         if (((lyap >= 2) && (mandel < 2)) || calclyap) {
  520.             /* first, sum the determinants */
  521.             if ((d = ABS(((*deriv)(x1, y1, params)).x))>=MAXFLOAT) {
  522.                 total += log(prod);
  523.                 prod = 1.0;
  524.                 total += log(d);
  525.             }
  526.             else {
  527.                 prod *= d;
  528.                 /* we need to prevent overflow and underflow */
  529.                 if ((prod > 1.0e12) || (prod < 1.0e-12)) {
  530.                     if (prod == 0.0)
  531.                         total += LN_MINDOUBLE;
  532.                     else if (prod < 0.0)
  533.                         total += LN_MAXDOUBLE;
  534.                     else
  535.                         total += log(prod);
  536.                     prod = 1.0;
  537.                 }
  538.             }
  539.             /* next, sum the traces */
  540.             if ((d = ABS(((*deriv)(x1, y1, params)).y))>=MAXFLOAT) {
  541.                 tr_total += log(tr_prod);
  542.                 tr_prod = 1.0;
  543.                 tr_total += log(d);
  544.             }
  545.             else {
  546.                 tr_prod *= d;
  547.                 /* we need to prevent overflow and underflow */
  548.                 if ((tr_prod > 1.0e12) || (tr_prod < 1.0e-12)) {
  549.                     if (tr_prod == 0.0)
  550.                         tr_total += LN_MINDOUBLE;
  551.                     else if (tr_prod < 0.0)
  552.                         tr_total += LN_MAXDOUBLE;
  553.                     else
  554.                         tr_total += log(tr_prod);
  555.                     tr_prod = 1.0;
  556.                 }
  557.             }
  558.         }
  559.         if (find) {    /* if we are looking for basins of attraction */
  560.           coordinates = (*map)(x1, y1, params);
  561.           for (j=0;j<numattrs;j++) {
  562.             if ((ABS(attrs[j][0] - x1) < delta) && 
  563.                 (ABS(attrs[j][1] - y1) < delta) &&
  564.                 (ABS(attrs[j][2] - coordinates.x) < cdelt) && 
  565.                 (ABS(attrs[j][3] - coordinates.y) < cdelt)) {
  566.                 if (foundit) {    /* we're close to two attractors */
  567.                     foundit = 0;
  568.                     break;
  569.                 }
  570.                 else /* we're within delta of a previous point */
  571.                     foundit=j+1;
  572.             }
  573.           }
  574.           if (foundit) {
  575.             basins[frame][perind[frame]] = foundit;
  576.             periods[frame][perind[frame]] = i + settle;
  577.           }
  578.         }
  579.         if (portrait) {
  580.           if (attractors) {     /* if we are drawing the attractors */
  581.             if (foundit)
  582.                   for (j=i;j<dwell;j++) {
  583.                     coordinates = (*map)(x1, y1, params);
  584.                     x1 = coordinates.x;
  585.                     y1 = coordinates.y;
  586.                     draw_portrait(x1,y1,color);
  587.                   }
  588.           }
  589.           else {    /* we are drawing orbits */
  590.             if (upper) 
  591.                 if (x1 < y1)
  592.                     draw_portrait(x1,y1,color);
  593.             if (lower) 
  594.                 if (x1 > y1)
  595.                     draw_portrait(x1,y1,color);
  596.             if (x1 == y1)
  597.                 draw_portrait(x1,y1,color);
  598.           }
  599.         }
  600.         if (foundit && (!attractors) && (!lyap)) {
  601.           if (calclyap) {
  602.             calc_lyapunov(prod,&total,tr_prod,&tr_total,&alpha,&beta,&lyapunov);
  603.             fprintf(stdout,"%lf\n",lyapunov);
  604.             calclyap = 0;
  605.             if (find == 2)
  606.                 find = 0;
  607.           }
  608.           break;
  609.         }
  610.     }
  611.     if (precrit == 3) { /* diagonal divides attractors */
  612.         if (x1 > y1)
  613.             indices[frame][perind[frame]] = 
  614.                     (periods[frame][perind[frame]] % second) + middle;
  615.         else {
  616.             indices[frame][perind[frame]] = 
  617.                     (periods[frame][perind[frame]] % first) + begin;
  618.             periods[frame][perind[frame]] *= -1;
  619.         }
  620.     }
  621.     else if (precrit == 7) { /* y-axis divides attractors */
  622.         if (x1 > 0)
  623.             indices[frame][perind[frame]] = 
  624.                 (periods[frame][perind[frame]] % second) + middle;
  625.         else {
  626.             indices[frame][perind[frame]] = 
  627.                 (periods[frame][perind[frame]] % first) + begin;
  628.             periods[frame][perind[frame]] *= -1;
  629.         }
  630.     }
  631.     else if (precrit == 8) { /* x-axis divides attractors */
  632.         if (y1 > 0)
  633.             indices[frame][perind[frame]] = 
  634.             (periods[frame][perind[frame]] % second) + middle;
  635.         else {
  636.             indices[frame][perind[frame]] = 
  637.             (periods[frame][perind[frame]] % first) + begin;
  638.                 periods[frame][perind[frame]] *= -1;
  639.         }
  640.     }
  641.     else if (precrit == 9)                /* color rate of attraction */
  642.         indices[frame][perind[frame]] = /* according to dist from origin */
  643.                 ((int)(radius*numfreecols/maxradius) % second) + middle;
  644.     else if (precrit == 10) {             /* color rate of attraction */
  645.         indices[frame][perind[frame]] = /* according to angle with x-axis */
  646.                 ((int)((atan2(y1,x1)+M_PI)*numfreecols/M_2PI) % second) 
  647.                     + middle - 1;
  648.     }
  649.     else if (precrit == 11) { /* x and y axes divide attractors */
  650.         if (y1 > 0) {
  651.               if (x1 > 0)
  652.                 indices[frame][perind[frame]] = 
  653.                           (periods[frame][perind[frame]] % second) + middle;
  654.               else 
  655.                 indices[frame][perind[frame]] = (second - 
  656.                         (periods[frame][perind[frame]]%second)-1)+middle;
  657.         }
  658.         else {
  659.               if (x1 > 0)
  660.                 indices[frame][perind[frame]] = 
  661.                         (periods[frame][perind[frame]] % first) + begin;
  662.               else
  663.                 indices[frame][perind[frame]] = 
  664.                         (first-(periods[frame][perind[frame]]%first)-1) + begin;
  665.         }
  666.     }
  667.     else if (!(find && (!foundit))) {
  668.         if (basins[frame][perind[frame]] && numattrs) {
  669.             j = numfreecols/numattrs;
  670.             indices[frame][perind[frame]]=(periods[frame][perind[frame]]%j) +
  671.                     ((basins[frame][perind[frame]]-1)*j) + STARTCOLOR;
  672.         }
  673.         else
  674.             indices[frame][perind[frame]] = 
  675.                     (periods[frame][perind[frame]] % second) + middle;
  676.     }
  677.     if ((!lyap) && find && (!foundit)) { /* didn't find a known attractor */
  678.         if (!((x1 == y1) && (precrit == 3))) { /* disallow diagonal points */
  679.           if (numattrs < MAXATTRS) {
  680.             attrs[numattrs][0] = x1;
  681.             attrs[numattrs][1] = y1;
  682.             coordinates = (*map)(x1, y1, params);
  683.             attrs[numattrs][2] = coordinates.x;
  684.             attrs[numattrs][3] = coordinates.y;
  685.             numattrs++;
  686.             fprintf(stderr,"Adding pair #%d ((%lf, %lf), (%lf, %lf))\n",
  687.                     numattrs,x1,y1,coordinates.x,coordinates.y);
  688.           }
  689.           else
  690.             fprintf(stderr,"Warning: found more than %d attractors\n",MAXATTRS);
  691.         }
  692.     }
  693.     lastx = x1; lasty = y1;
  694.     if ((lyap >= 2) && (mandel < 2)) {
  695.         calc_lyapunov(prod,&total,tr_prod,&tr_total,&alpha,&beta,&lyapunov);
  696.         if (lyapunov > 0) {
  697.             index = periods[frame][perind[frame]] = (int)(lyapunov*first);
  698.             indices[frame][perind[frame]]=(index % first)+begin;
  699.         }
  700.         else {
  701.             index = periods[frame][perind[frame]] = (int)(-lyapunov*second);
  702.             indices[frame][perind[frame]] = (index % second) + middle;
  703.         }
  704.         if (lyapunov > maxexp)
  705.             maxexp = lyapunov;
  706.         if (lyapunov < minexp)
  707.             minexp = lyapunov;
  708.           BufferPoint(dpy,which,pixmap,Data_GC,&Points,
  709.                 indices[frame][perind[frame]++],point.x,height-point.y-1);
  710.     }
  711.     else if ((lyap >= 2) && (mandel == 2)) { /* draw Mandelbrot black */
  712.         periods[frame][perind[frame]]=indices[frame][perind[frame]] = 0;
  713.           BufferPoint(dpy,which,pixmap,Data_GC,&Points,
  714.                 indices[frame][perind[frame]++],point.x,height-point.y-1);
  715.     }
  716.     else if ((lyap >= 2) && (mandel == 3)) { /* display periods */
  717.         if (foundit)
  718.             periods[frame][perind[frame]] = settle + i;
  719.         else
  720.                periods[frame][perind[frame]] = 0;
  721.         indices[frame][perind[frame]] = 
  722.             (periods[frame][perind[frame]]%second)+middle;
  723.           BufferPoint(dpy,which,pixmap,Data_GC,&Points,
  724.                 indices[frame][perind[frame]++],point.x,height-point.y-1);
  725.     }
  726.     else if (!lyap) {
  727.           BufferPoint(dpy,which,pixmap,Data_GC,&Points,
  728.             indices[frame][perind[frame]++],point.x,height-point.y-1);
  729.     }
  730.     else
  731.         perind[frame]++;
  732.     if (portrait)
  733.         FlushBuffer(dpy, trajec, pixtra, Data_GC, &Orbits, color, color+1);
  734.     lastx = x1; lasty = y1;
  735.     return FALSE;
  736. }
  737.  
  738. set_arc(x1, y1)
  739. double x1, y1;
  740. {
  741.     extern void draw_critical();
  742.  
  743.     draw_critical(crijec, x1,y1,1);
  744.     draw_critical(crijec, x,y,1);
  745.     A0.x = x; A0.y = y;
  746.     A1.x = x1; A1.y = y1;
  747.     if (!found_arc) {
  748.         B0.x = x; B0.y = y;
  749.         B1.x = x1; B1.y = y1;
  750.         found_arc = 1;
  751.         if (find_arcs(B0, B1) == TRUE) {
  752.             found_arc = 1;
  753.             if (!idown)
  754.                 Show_Info();
  755.         }
  756.         else
  757.             found_arc = 0;
  758.     }
  759.     else {
  760.         if ((B0.x == A0.x) && (B0.y == A0.y)) {
  761.             if (((x - B1.x)*(x - B1.x) + (y - B1.y)*(y - B1.y)) >
  762.                ((B0.x-B1.x)*(B0.x-B1.x) + (B0.y-B1.y)*(B0.y-B1.y))) {
  763.                 B0.x = x; B0.y = y;
  764.                 B1.x = x1; B1.y = y1;
  765.                 if (find_arcs(B0, B1) == TRUE) {
  766.                     found_arc = 1;
  767.                     if (!idown)
  768.                         Show_Info();
  769.                 }
  770.                 return;
  771.             }
  772.         }
  773.         if (find_arcs(A0, A1) == TRUE) {
  774.             found_arc++;
  775.             if (!idown)
  776.                 Show_Info();
  777.         }
  778.     }
  779. }
  780.  
  781. find_arcs(p0, p1)
  782. pair p0, p1;
  783. {
  784.     static int i;
  785.     static double mix=0.0, max=0.0, miy=0.0, may=0.0;
  786.     extern int found_arc;
  787.     extern void draw_critical();
  788.  
  789.     if (!found_arc) {
  790.         fprintf(stderr,"No intersections detected\n");
  791.         return FALSE;
  792.     }
  793.     if (p0.x < p1.x) {
  794.         mix = p0.x; max = p1.x;
  795.     }
  796.     else if (p0.x > p1.x) {
  797.         mix = p1.x; max = p0.x;
  798.     }
  799.     if (p0.y < p1.y) {
  800.         miy = p0.y; may = p1.y;
  801.     }
  802.     else if (p0.y > p1.y) {
  803.         miy = p1.y; may = p0.y;
  804.     }
  805.     if ((p0.x == p1.x) && (p0.y == p1.y)) {
  806.         fprintf(stderr,"Found a fixed point at (%lf, %lf)\n",p0.x, p0.y);
  807.         return FALSE;
  808.     }
  809.     numarcs = 0;
  810.     for (i=0;i<numcrits;i++) {
  811.         if ((crit_pts[0][i] >= mix) && (crit_pts[0][i] <= max) &&
  812.             (crit_pts[1][i] >= miy) && (crit_pts[1][i] <= may)) {
  813.             draw_critical(crijec, crit_pts[0][i],crit_pts[1][i],1);
  814.             n_crit_arc[0][numarcs] = crit_arc[0][numarcs] = crit_pts[0][i];
  815.             n_crit_arc[1][numarcs] = crit_arc[1][numarcs++] = crit_pts[1][i];
  816.         }
  817.     }
  818.     return TRUE;
  819. }
  820.  
  821.