Subversion Repositories wimsdev

Rev

Rev 7615 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
10 reyssat 1
/*
2
  #define _ISOC99_SOURCE
3
 
4
  This does not work on glibc 2.0.x --- so sorry.
5
*/
6
#include <sys/time.h>
3819 kbelabas 7
#include <time.h>
10 reyssat 8
#include <errno.h>
9
#include <stdio.h>
10
#include <stdlib.h>
11
#include <math.h>
12
#include <string.h>
13
#include <unistd.h>
7615 bpr 14
#include <gd.h>
10 reyssat 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;
7615 bpr 74
 
10 reyssat 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++)
7676 bpr 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
          }
10 reyssat 92
      if (!found) break;
93
      for (i=0; i<XDIV-div; i++)
7676 bpr 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;
10 reyssat 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;
7615 bpr 127
 
10 reyssat 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)
7676 bpr 140
      dt=XSTEP/xp;
10 reyssat 141
      if (yp*dt>YSTEP)
7676 bpr 142
      dt=YSTEP/yp;
10 reyssat 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)
7676 bpr 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
      }
10 reyssat 188
      if (dt>0 && steps==1)
7676 bpr 189
      dir_angle=atan2(y-y0, x-x0);
10 reyssat 190
      gdImageLine(img,
7676 bpr 191
                lrint((lx-XMIN)*gxden), height-1-lrint((ly-YMIN)*gyden),
192
                lrint((x-XMIN)*gxden), height-1-lrint((y-YMIN)*gyden), black);
10 reyssat 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 */
7676 bpr 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
      }
10 reyssat 222
    }
223
  while (get_point(&x0, &y0)>=0);
224
}
225
 
226
void show_usage(void)
227
{
228
  fprintf(stderr,
7676 bpr 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
        );
10 reyssat 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 (! finite(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':
7676 bpr 277
        show_usage();
278
        return 0;
279
        break;
10 reyssat 280
      case 'v':
7676 bpr 281
        fprintf(stderr, "DrawODE version 1.0\n");
282
        return 0;
283
        break;
10 reyssat 284
      case 'o':
7676 bpr 285
        out_fname=optarg;
286
        break;
10 reyssat 287
      case 'w':
7676 bpr 288
        width=atoi(optarg);
289
        break;
10 reyssat 290
      case 'l':
7676 bpr 291
        height=atoi(optarg);
292
        break;
10 reyssat 293
      case 's':
7676 bpr 294
        width=height=atoi(optarg);
295
        break;
10 reyssat 296
      case 'f':
7676 bpr 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;
10 reyssat 304
      case 'p':
7676 bpr 305
        read_params(optarg);
306
        break;
10 reyssat 307
      }
7615 bpr 308
 
10 reyssat 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