Subversion Repositories wimsdev

Rev

Rev 7676 | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. /*
  2.   #define _ISOC99_SOURCE
  3.  
  4.   This does not work on glibc 2.0.x --- so sorry.
  5. */
  6. #include <sys/time.h>
  7. #include <time.h>
  8. #include <errno.h>
  9. #include <stdio.h>
  10. #include <stdlib.h>
  11. #include <math.h>
  12. #include <string.h>
  13. #include <unistd.h>
  14. #include <gd.h>
  15. #include "drawode.h"
  16.  
  17. #define NTEST 50
  18. #define NPARAM 4
  19. #define MINDIV 2
  20. #define MAXCOUNT 50
  21. #define XDIV 32
  22. #define YDIV 32
  23. #define XMIN (-zoom)
  24. #define YMIN (-zoom)
  25. #define XMAX (zoom)
  26. #define YMAX (zoom)
  27. #define XSTEP ((XMAX-XMIN)/XDIV)
  28. #define YSTEP ((YMAX-YMIN)/YDIV)
  29. #define XDEN (1.0/XSTEP)
  30. #define YDEN (1.0/YSTEP)
  31. #define ARROW_LEN 8
  32. #define ARROW_ANGLE 0.3
  33. #define MIN_STEP 2 /* minimum steps to draw arrows */
  34.  
  35. static char pr[XDIV][YDIV]; /* footprint */
  36. static char xr[XDIV][YDIV], nr[XDIV][YDIV]; /* used by get_point() */
  37.  
  38. static double dir_angle; /* origin is x0, y0 */
  39. double values[NPARAM];
  40. static double zoom=1.0;
  41. static int width=200, height=200;
  42. static gdImagePtr img;
  43. static int black, white;
  44. static const char *out_fname=NULL;
  45. static int func_num=0;
  46.  
  47. static inline int transx(double x)
  48. {
  49.   return lrint((x-XMIN)*XDEN);
  50. }
  51.  
  52. static inline int transy(double y)
  53. {
  54.   return lrint((y-YMIN)*YDEN);
  55. }
  56.  
  57. /* returns 0 for success, -1 for out-of-bound or already-gone-to (x,y) */
  58. static inline int footprint(int px, int py)
  59. {
  60.   if (px>=0 && px<XDIV && py>=0 && py<YDIV && !pr[px][py])
  61.     {
  62.       pr[px][py]=1;
  63.       return 0;
  64.     }
  65.   else return -1;
  66. }
  67.  
  68. /* returns 0 for success, -1 for fail */
  69. static int get_point(double *x, double *y)
  70. {
  71.   int div, i, j;
  72.   int mdiv, mi=0, mj=0;
  73.   int found;
  74.  
  75.   mdiv=0;
  76.   div=1;
  77.   memcpy(xr, pr, sizeof(pr));
  78.   while (1)
  79.     {
  80.       if (div>XDIV) break;
  81.       found=0;
  82.       for (i=0; i<=XDIV-div; i++)
  83.       for (j=0; j<=YDIV-div; j++)
  84.         if (xr[i][j]==0)
  85.           {
  86.             found=1;
  87.             mdiv=div;
  88.             mi=i;
  89.             mj=j;
  90.             break;
  91.           }
  92.       if (!found) break;
  93.       for (i=0; i<XDIV-div; i++)
  94.       for (j=0; j<YDIV-div; j++)
  95.         nr[i][j]=(xr[i][j]+xr[i][j+1]+xr[i+1][j]+xr[i+1][j+1]) ? 1 : 0;
  96.       memcpy(xr, nr, sizeof(nr));
  97.       div++;
  98.     }
  99.   if (mdiv>=MINDIV)
  100.     {
  101.       *x=XMIN+(mi+(mdiv-1)*0.5)*XSTEP;
  102.       *y=YMIN+(mj+(mdiv-1)*0.5)*YSTEP;
  103.       return 0;
  104.     }
  105.   else return -1;
  106. }
  107.  
  108. static void draw_init(void)
  109. {
  110.   int i, j;
  111.   for (i=0; i<XDIV; i++)
  112.     for (j=0; j<YDIV; j++)
  113.       pr[i][j]=0;
  114. }
  115.  
  116. static double frand(double min, double max)
  117. {
  118.   return min+(max-min)*(double)rand()/RAND_MAX;
  119. }
  120.  
  121. static double get_dt(void)
  122. {
  123.   func_t dx, dy;
  124.   double dt;
  125.   double x, y, xp, yp;
  126.   int i;
  127.  
  128.   dx=funcs[func_num][0];
  129.   dy=funcs[func_num][1];
  130.   dt=1.0; /* enough, eh? */
  131.  
  132.   srand(time(NULL));
  133.   for (i=0; i<NTEST; i++)
  134.     {
  135.       x=frand(XMIN, XMAX);
  136.       y=frand(YMIN, YMAX);
  137.       xp=fabs(dx(x, y));
  138.       yp=fabs(dy(x, y));
  139.       if (xp*dt>XSTEP)
  140.       dt=XSTEP/xp;
  141.       if (yp*dt>YSTEP)
  142.       dt=YSTEP/yp;
  143.     }
  144.   return dt*0.5; /* for precision */
  145. }
  146.  
  147. static int do_loop(double x0, double y0, double dt)
  148. {
  149.   func_t dx, dy;
  150.   int i, count=0, steps=0;
  151.   int px, py;
  152.   double x, y, lx, ly;
  153.   double xa[5], ya[5];
  154.   double gxden, gyden;
  155.  
  156.   dx=funcs[func_num][0];
  157.   dy=funcs[func_num][1];
  158.   x=x0;
  159.   y=y0;
  160.   gxden=width/(XMAX-XMIN);
  161.   gyden=height/(YMAX-YMIN);
  162.   /*  fprintf(stderr, "Start: (%0.2f, %0.2f)\n", x0, y0); */
  163.   while (1)
  164.     {
  165.       steps++;
  166.       px=transx(x);
  167.       py=transy(y);
  168.       if (footprint(px, py)<0 && steps>1) break;
  169.       /* fprintf(stderr, "Step %d: footprint(%d, %d)\n", steps, px, py); */
  170.       lx=x;
  171.       ly=y;
  172.       count=0;
  173.       while (count<MAXCOUNT)
  174.       {
  175.         count++;
  176.         xa[0]=x;
  177.         ya[0]=y;
  178.         for (i=1; i<=4; i++)
  179.           {
  180.             xa[i]=xa[0]+dx(xa[i-1], ya[i-1])*dt;
  181.             ya[i]=ya[0]+dy(xa[i-1], ya[i-1])*dt;
  182.           }
  183.         x=(xa[1]+2.0*xa[2]+xa[3]-xa[4])*(1/3.0);
  184.         y=(ya[1]+2.0*ya[2]+ya[3]-ya[4])*(1/3.0);
  185.         if (transx(x)!=px) break;
  186.         if (transy(y)!=py) break;
  187.       }
  188.       if (dt>0 && steps==1)
  189.       dir_angle=atan2(y-y0, x-x0);
  190.       gdImageLine(img,
  191.                 lrint((lx-XMIN)*gxden), height-1-lrint((ly-YMIN)*gyden),
  192.                 lrint((x-XMIN)*gxden), height-1-lrint((y-YMIN)*gyden), black);
  193.     }
  194.   return steps;
  195. }
  196.  
  197. static void draw_main(double dt)
  198. {
  199.   double x0, y0;
  200.   int px0, py0, px1, py1, px2, py2, px3, py3;
  201.   int step1, step2;
  202.  
  203.   x0=0;
  204.   y0=0;
  205.   draw_init();
  206.   do {
  207.       step1=do_loop(x0, y0, dt);
  208.       step2=do_loop(x0, y0, -dt);
  209.       if (step1>=MIN_STEP+1 && step2>=MIN_STEP) /* draw an arrow */
  210.       {
  211.         px0=lrint((x0-XMIN)/(XMAX-XMIN)*width);
  212.         py0=lrint((y0-YMIN)/(YMAX-YMIN)*height);
  213.         px1=px0+lrint(ARROW_LEN*cos(dir_angle));
  214.         py1=py0+lrint(ARROW_LEN*sin(dir_angle));
  215.         px2=px1-lrint(ARROW_LEN*cos(dir_angle-ARROW_ANGLE));
  216.         py2=py1-lrint(ARROW_LEN*sin(dir_angle-ARROW_ANGLE));
  217.         px3=px1-lrint(ARROW_LEN*cos(dir_angle+ARROW_ANGLE));
  218.         py3=py1-lrint(ARROW_LEN*sin(dir_angle+ARROW_ANGLE));
  219.         gdImageLine(img, px1, height-1-py1, px2, height-1-py2, black);
  220.         gdImageLine(img, px1, height-1-py1, px3, height-1-py3, black);
  221.       }
  222.     }
  223.   while (get_point(&x0, &y0)>=0);
  224. }
  225.  
  226. void show_usage(void)
  227. {
  228.   fprintf(stderr,
  229.         "Usage: drawode [arguments...]\n\n"
  230.         "Arguments:\n"
  231.         " -h            Help\n"
  232.         " -o <fname>    Set output file name\n"
  233.         " -w <width>    Set width\n"
  234.         " -l <height>   Set height\n"
  235.         " -s <size>     Set width and height\n"
  236.         " -f <num>      Set function number\n"
  237.         " -p <params>   Set parameters (comma separated)\n"
  238.         "\nValid function numbers:\n"
  239.         " 0             x'=y y'=x(x+y)+ax+by\n"
  240.         " 1             x'=ax+by y'=bx+ay\n"
  241.         " 2             x'=ax+by y'=cx+dy\n"
  242.         );
  243. }
  244.  
  245. void read_params(const char *args)
  246. {
  247.   int i;
  248.   char *p;
  249.  
  250.   i=0;
  251.   while (i<NPARAM && (*args))
  252.     {
  253.       values[i]=strtod(args, &p);
  254.       if (*p!=',') break;
  255.       args=p+1;
  256.       i++;
  257.     }
  258.   if (! isfinite(values[i]))
  259.     {
  260.       fprintf(stderr, "Invalid parameter value.\n");
  261.       exit(1);
  262.     }
  263. }
  264.  
  265. int main(int argc, char **argv)
  266. {
  267.   FILE *file;
  268.   int c;
  269.   int i;
  270.  
  271.   for (i=0; i<NPARAM; i++) values[i]=0.0;
  272.   /* read arguments */
  273.   while ((c=getopt(argc, argv, "ho:w:l:s:f:p:v")) != -1)
  274.     switch (c)
  275.       {
  276.       case 'h':
  277.         show_usage();
  278.         return 0;
  279.         break;
  280.       case 'v':
  281.         fprintf(stderr, "DrawODE version 1.0\n");
  282.         return 0;
  283.         break;
  284.       case 'o':
  285.         out_fname=optarg;
  286.         break;
  287.       case 'w':
  288.         width=atoi(optarg);
  289.         break;
  290.       case 'l':
  291.         height=atoi(optarg);
  292.         break;
  293.       case 's':
  294.         width=height=atoi(optarg);
  295.         break;
  296.       case 'f':
  297.         func_num=atoi(optarg);
  298.         if (func_num<0 || func_num>=nfunc)
  299.         {
  300.           fprintf(stderr, "Invalid function number %d (max %d).\n", func_num, nfunc-1);
  301.           return 1;
  302.         }
  303.         break;
  304.       case 'p':
  305.         read_params(optarg);
  306.         break;
  307.       }
  308.  
  309.   if (out_fname==NULL)
  310.     {
  311.       fprintf(stderr, "You need to specify a file name.\n");
  312.       return 1;
  313.     }
  314.   if (width<0 || height<0)
  315.     {
  316.       fprintf(stderr, "Invalid size.\n");
  317.       return 1;
  318.     }
  319.   file=fopen(out_fname, "wb");
  320.   if (file==NULL)
  321.     {
  322.       perror(out_fname);
  323.       return 1;
  324.     }
  325.   img=gdImageCreate(width, height);
  326.   white=gdImageColorAllocate(img, 255, 255, 255);
  327.   black=gdImageColorAllocate(img, 0, 0, 0);
  328.   draw_main(get_dt());
  329.   gdImageGif(img, file);
  330.   fclose(file);
  331.   return 0;
  332. }
  333.  
  334.  
  335.  
  336.