Back to index

wims  3.65+svn20090927
drawode.c
Go to the documentation of this file.
00001 /*
00002   #define _ISOC99_SOURCE
00003 
00004   This does not work on glibc 2.0.x --- so sorry.
00005 */
00006 #include <sys/time.h>
00007 #include <errno.h>
00008 #include <stdio.h>
00009 #include <stdlib.h>
00010 #include <math.h>
00011 #include <string.h>
00012 #include <unistd.h>
00013 #include "gd.h"
00014 #include "drawode.h"
00015 
00016 #define NTEST 50
00017 #define NPARAM 4
00018 #define MINDIV 2
00019 #define MAXCOUNT 50
00020 #define XDIV 32
00021 #define YDIV 32
00022 #define XMIN (-zoom)
00023 #define YMIN (-zoom)
00024 #define XMAX (zoom)
00025 #define YMAX (zoom)
00026 #define XSTEP ((XMAX-XMIN)/XDIV)
00027 #define YSTEP ((YMAX-YMIN)/YDIV)
00028 #define XDEN (1.0/XSTEP)
00029 #define YDEN (1.0/YSTEP)
00030 #define ARROW_LEN 8
00031 #define ARROW_ANGLE 0.3
00032 #define MIN_STEP 2 /* minimum steps to draw arrows */
00033 
00034 static char pr[XDIV][YDIV]; /* footprint */
00035 static char xr[XDIV][YDIV], nr[XDIV][YDIV]; /* used by get_point() */
00036 
00037 static double dir_angle; /* origin is x0, y0 */
00038 double values[NPARAM];
00039 static double zoom=1.0;
00040 static int width=200, height=200;
00041 static gdImagePtr img;
00042 static int black, white;
00043 static const char *out_fname=NULL;
00044 static int func_num=0;
00045 
00046 static inline int transx(double x)
00047 {
00048   return lrint((x-XMIN)*XDEN);
00049 }
00050 
00051 static inline int transy(double y)
00052 {
00053   return lrint((y-YMIN)*YDEN);
00054 }
00055 
00056 /* returns 0 for success, -1 for out-of-bound or already-gone-to (x,y) */
00057 static inline int footprint(int px, int py)
00058 {
00059   if (px>=0 && px<XDIV && py>=0 && py<YDIV && !pr[px][py])
00060     {
00061       pr[px][py]=1;
00062       return 0;
00063     }
00064   else return -1;
00065 }
00066 
00067 /* returns 0 for success, -1 for fail */
00068 static int get_point(double *x, double *y)
00069 {
00070   int div, i, j;
00071   int mdiv, mi=0, mj=0;
00072   int found;
00073   
00074   mdiv=0;
00075   div=1;
00076   memcpy(xr, pr, sizeof(pr));
00077   while (1)
00078     {
00079       if (div>XDIV) break;
00080       found=0;
00081       for (i=0; i<=XDIV-div; i++)
00082        for (j=0; j<=YDIV-div; j++)
00083          if (xr[i][j]==0)
00084            {
00085              found=1;
00086              mdiv=div;
00087              mi=i;
00088              mj=j;
00089              break;
00090            }
00091       if (!found) break;
00092       for (i=0; i<XDIV-div; i++)
00093        for (j=0; j<YDIV-div; j++)
00094          nr[i][j]=(xr[i][j]+xr[i][j+1]+xr[i+1][j]+xr[i+1][j+1]) ? 1 : 0;
00095       memcpy(xr, nr, sizeof(nr));
00096       div++;
00097     }
00098   if (mdiv>=MINDIV)
00099     {
00100       *x=XMIN+(mi+(mdiv-1)*0.5)*XSTEP;
00101       *y=YMIN+(mj+(mdiv-1)*0.5)*YSTEP;
00102       return 0;
00103     }
00104   else return -1;
00105 }
00106 
00107 static void draw_init(void)
00108 {
00109   int i, j;
00110   for (i=0; i<XDIV; i++)
00111     for (j=0; j<YDIV; j++)
00112       pr[i][j]=0;
00113 }
00114 
00115 static double frand(double min, double max)
00116 {
00117   return min+(max-min)*(double)rand()/RAND_MAX;
00118 }
00119 
00120 static double get_dt(void)
00121 {
00122   func_t dx, dy;
00123   double dt;
00124   double x, y, xp, yp;
00125   int i;
00126   
00127   dx=funcs[func_num][0];
00128   dy=funcs[func_num][1];
00129   dt=1.0; /* enough, eh? */
00130 
00131   srand(time(NULL));
00132   for (i=0; i<NTEST; i++)
00133     {
00134       x=frand(XMIN, XMAX);
00135       y=frand(YMIN, YMAX);
00136       xp=fabs(dx(x, y));
00137       yp=fabs(dy(x, y));
00138       if (xp*dt>XSTEP)
00139        dt=XSTEP/xp;
00140       if (yp*dt>YSTEP)
00141        dt=YSTEP/yp;
00142     }
00143   return dt*0.5; /* for precision */
00144 }
00145 
00146 static int do_loop(double x0, double y0, double dt)
00147 {
00148   func_t dx, dy;
00149   int i, count=0, steps=0;
00150   int px, py;
00151   double x, y, lx, ly;
00152   double xa[5], ya[5];
00153   double gxden, gyden;
00154 
00155   dx=funcs[func_num][0];
00156   dy=funcs[func_num][1];
00157   x=x0;
00158   y=y0;
00159   gxden=width/(XMAX-XMIN);
00160   gyden=height/(YMAX-YMIN);
00161   /*  fprintf(stderr, "Start: (%0.2f, %0.2f)\n", x0, y0); */
00162   while (1)
00163     {
00164       steps++;
00165       px=transx(x);
00166       py=transy(y);
00167       if (footprint(px, py)<0 && steps>1) break;
00168       /* fprintf(stderr, "Step %d: footprint(%d, %d)\n", steps, px, py); */
00169       lx=x;
00170       ly=y;
00171       count=0;
00172       while (count<MAXCOUNT)
00173        {
00174          count++;
00175          xa[0]=x;
00176          ya[0]=y;
00177          for (i=1; i<=4; i++)
00178            {
00179              xa[i]=xa[0]+dx(xa[i-1], ya[i-1])*dt;
00180              ya[i]=ya[0]+dy(xa[i-1], ya[i-1])*dt;
00181            }
00182          x=(xa[1]+2.0*xa[2]+xa[3]-xa[4])*(1/3.0);
00183          y=(ya[1]+2.0*ya[2]+ya[3]-ya[4])*(1/3.0);
00184          if (transx(x)!=px) break;
00185          if (transy(y)!=py) break;
00186        }
00187       if (dt>0 && steps==1)
00188        dir_angle=atan2(y-y0, x-x0);
00189       gdImageLine(img,
00190                   lrint((lx-XMIN)*gxden), height-1-lrint((ly-YMIN)*gyden),
00191                   lrint((x-XMIN)*gxden), height-1-lrint((y-YMIN)*gyden), black);
00192     }
00193   return steps;
00194 }
00195 
00196 static void draw_main(double dt)
00197 {
00198   double x0, y0;
00199   int px0, py0, px1, py1, px2, py2, px3, py3;
00200   int step1, step2;
00201 
00202   x0=0;
00203   y0=0;
00204   draw_init();
00205   do {
00206       step1=do_loop(x0, y0, dt);
00207       step2=do_loop(x0, y0, -dt);
00208       if (step1>=MIN_STEP+1 && step2>=MIN_STEP) /* draw an arrow */
00209        {
00210          px0=lrint((x0-XMIN)/(XMAX-XMIN)*width);
00211          py0=lrint((y0-YMIN)/(YMAX-YMIN)*height);
00212          px1=px0+lrint(ARROW_LEN*cos(dir_angle));
00213          py1=py0+lrint(ARROW_LEN*sin(dir_angle));         
00214          px2=px1-lrint(ARROW_LEN*cos(dir_angle-ARROW_ANGLE));
00215          py2=py1-lrint(ARROW_LEN*sin(dir_angle-ARROW_ANGLE));
00216          px3=px1-lrint(ARROW_LEN*cos(dir_angle+ARROW_ANGLE));
00217          py3=py1-lrint(ARROW_LEN*sin(dir_angle+ARROW_ANGLE));
00218          gdImageLine(img, px1, height-1-py1, px2, height-1-py2, black);
00219          gdImageLine(img, px1, height-1-py1, px3, height-1-py3, black);
00220        }
00221     }
00222   while (get_point(&x0, &y0)>=0);
00223 }
00224 
00225 void show_usage(void)
00226 {
00227   fprintf(stderr,
00228          "Usage: drawode [arguments...]\n\n"
00229          "Arguments:\n"
00230          " -h            Help\n"
00231          " -o <fname>    Set output file name\n"
00232          " -w <width>    Set width\n"
00233          " -l <height>   Set height\n"
00234          " -s <size>     Set width and height\n"
00235          " -f <num>      Set function number\n"
00236          " -p <params>   Set parameters (comma separated)\n"
00237          "\nValid function numbers:\n"
00238          " 0             x'=y y'=x(x+y)+ax+by\n"
00239          " 1             x'=ax+by y'=bx+ay\n"
00240          " 2             x'=ax+by y'=cx+dy\n"
00241          );
00242 }
00243 
00244 void read_params(const char *args)
00245 {
00246   int i;
00247   char *p;
00248 
00249   i=0;
00250   while (i<NPARAM && (*args))
00251     {
00252       values[i]=strtod(args, &p);
00253       if (*p!=',') break;
00254       args=p+1;
00255       i++;
00256     }
00257   if (! finite(values[i]))
00258     {
00259       fprintf(stderr, "Invalid parameter value.\n");
00260       exit(1);
00261     }
00262 }
00263 
00264 int main(int argc, char **argv)
00265 {
00266   FILE *file;
00267   int c;
00268   int i;
00269 
00270   for (i=0; i<NPARAM; i++) values[i]=0.0;
00271   /* read arguments */
00272   while ((c=getopt(argc, argv, "ho:w:l:s:f:p:v")) != -1)
00273     switch (c)
00274       {
00275       case 'h':
00276        show_usage();
00277        return 0;
00278        break;
00279       case 'v':
00280        fprintf(stderr, "DrawODE version 1.0\n");
00281        return 0;
00282        break;
00283       case 'o':
00284        out_fname=optarg;
00285        break;
00286       case 'w':
00287        width=atoi(optarg);
00288        break;
00289       case 'l':
00290        height=atoi(optarg);
00291        break;
00292       case 's':
00293        width=height=atoi(optarg);
00294        break;
00295       case 'f':
00296        func_num=atoi(optarg);
00297        if (func_num<0 || func_num>=nfunc)
00298          {
00299            fprintf(stderr, "Invalid function number %d (max %d).\n", func_num, nfunc-1);
00300            return 1;
00301          }
00302        break;
00303       case 'p':
00304        read_params(optarg);
00305        break;
00306       }
00307   
00308   if (out_fname==NULL)
00309     {
00310       fprintf(stderr, "You need to specify a file name.\n");
00311       return 1;
00312     }
00313   if (width<0 || height<0)
00314     {
00315       fprintf(stderr, "Invalid size.\n");
00316       return 1;
00317     }
00318   file=fopen(out_fname, "wb");
00319   if (file==NULL)
00320     {
00321       perror(out_fname);
00322       return 1;
00323     }
00324   img=gdImageCreate(width, height);
00325   white=gdImageColorAllocate(img, 255, 255, 255);
00326   black=gdImageColorAllocate(img, 0, 0, 0);
00327   draw_main(get_dt());
00328   gdImageGif(img, file);
00329   fclose(file);
00330   return 0;
00331 }
00332 
00333 
00334