Subversion Repositories wimsdev

Rev

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