dfe.c

   1 #define _GNU_SOURCE
   2 #include <stdio.h>
   3 #include <math.h>
   4 #include <time.h>
   5 #include <stdlib.h>
   6
   7 /*
   8  * compile: gcc -Wall -o dfe dfe.c bmp.c list.c -O3 -lm
   9  *
  10  * usage: ./dfe 100 50 2 1 1 10
  11  *
  12  * author: christian leirer, frank zirkelbach
  13  *
  14  * need api: cvs -d:pserver:anonymous@hackdaworld.dyndns.org:/my-code co api
  15  *
  16  */
  17
  18 #include "list.h"
  19 #include "bmp.h"
  20
  21 typedef struct s_gp {
  22   int x,y;
  23   unsigned char status;
  24 #define FLUESSIG 0
  25 #define FEST 1
  26   double fx,fy;
  27   double op;
  28   double pp;
  29 } t_gp;
  30
  31 #define RANGE 4
  32
  33 #define X_ 50
  34 #define Y_ 10
  35
  36 #define STEPS 100
  37 #define PLOT 1
  38
  39 int main(int argc,char **argv) {
  40
  41   int x,y;
  42   int dx,dy;
  43   int mx,my;
  44   double lf[2];
  45   double sf[2];
  46   double mf;
  47   double lpp;
  48   int i;
  49   int size;
  50   t_gp gp;
  51   t_gp *gp_ptr;
  52   t_list list;
  53   t_list_element *tmp;
  54   t_bmp bmp;
  55   char string[128];
  56   int fd;
  57
  58   double by;
  59   double pp;
  60   double xi;
  61   double rs;
  62   double co;
  63   double bf;
  64
  65   /* parse argv */
  66
  67   if(argc!=7) {
  68     printf("usage: %s <b-feld> <p-pot> <xi> <rand-scale> <coulomb> <b-force>\n",
  69            argv[0]);
  70     return -1;
  71   }
  72
  73   by=atof(argv[1]);
  74   pp=atof(argv[2]);
  75   xi=atof(argv[3]);
  76   rs=atof(argv[4]);
  77   co=atof(argv[5]);
  78   bf=atof(argv[6]);
  79
  80   /* init */
  81
  82   fd=open("/dev/null",O_WRONLY);
  83
  84   list_init(&list,fd);
  85
  86   size=X_*Y_;
  87   bmp_init(&bmp,1);
  88   bmp.width=X_;
  89   bmp.height=Y_;
  90   bmp.mode=WRITE;
  91   bmp_alloc_map(&bmp);
  92
  93   srandom(time(NULL));
  94
  95   for(x=0;x<X_;x++) {
  96     memset(&gp,0,sizeof(t_gp));
  97     //gp.status=FEST;
  98     gp.x=x; gp.y=0;
  99     list_add_element(&list,&gp,sizeof(t_gp));
 100   }
 101
 102   /* mainloop */
 103   for(i=0;i<STEPS;i++) {
 104     /* determine status of vortex, liquid or solid */
 105     list_reset(&list);
 106     do {
 107       lpp=0;
 108       gp_ptr=(t_gp *)list.current->data;
 109       tmp=list.current;
 110       for(dx=-RANGE;dx<=RANGE;dx++) {
 111         for(dy=-RANGE;dy<=RANGE;dy++) {
 112           if((dx!=0)||(dy!=0)) {
 113             x=(gp_ptr->x+X_+dx)%X_;
 114             y=gp_ptr->y+dy;
 115             if(y>=0&&y<Y_) {
 116               gp.x=x; gp.y=y;
 117               if(list_search_data(&list,&gp,2*sizeof(int))==L_SUCCESS)
 118                 gp_ptr->pp+=exp(-1.0*sqrt(x*x+y*y)/xi);
 119
 120             }
 121           }
 122         }
 123       }
 124       if(lpp>(0.8*pp)) gp_ptr->status=FLUESSIG;
 125       else gp_ptr->status=FEST;
 126       list.current=tmp;
 127     }
 128     while(list_next(&list)!=L_NO_NEXT_ELEMENT);
 129
 130     /* force on vortex */
 131     list_reset(&list);
 132     do {
 133       lf[0]=0; lf[1]=0;
 134       sf[0]=0; sf[1]=0;
 135
 136       tmp=list.current;
 137       gp_ptr=(t_gp *)list.current->data;
 138       for(dx=-RANGE;dx<=RANGE;dx++) {
 139         for(dy=-RANGE;dy<=RANGE;dy++) {
 140           if((dx!=0)||(dy!=0)) {
 141             x=(gp_ptr->x+X_+dx)%X_;
 142             y=gp_ptr->y+dy;
 143             if(y>=0&&y<Y_) {
 144               gp.x=x; gp.y=y;
 145               if(list_search_data(&list,&gp,2*sizeof(int))==L_SUCCESS) {
 146                 memcpy(&gp,list.current->data,sizeof(t_gp));
 147                 if(gp.status&FEST) {
 148                   if(dy==0) sf[0]+=(-1.0/(dx*dx+dy*dy)*co*(dx/abs(dx)));
 149                   else if (dx==0) sf[1]+=(-1.0/(dx*dx+dy*dy)*co*(dy/abs(dy)));
 150                   /* hier ist die scheisse drin!! */
 151                   else {
 152                     sf[0]+=(-1.0/(dx*dx+dy*dy)*co*dx/abs(dy));
 153                     sf[1]+=(-1.0/(dx*dx+dy*dy)*co*dy/abs(dx));
 154                   }
 155                 }
 156                 else {
 157                   /* und hier wahrscheinlich auch */
 158                   if(dy==0) lf[0]+=(-1.0*(dx/abs(dx))*(by+co));
 159                   else if(dx==0) lf[1]+=(-1.0*(dy/abs(dy))*(by+co));
 160                   else{
 161                     lf[0]+=(-1.0*dx/abs(dy)*(by+co));
 162                     lf[1]+=(-1.0*dy/abs(dx)*(by+co));
 163                   }
 164                 }
 165               }
 166             }
 167           }
 168         }
 169       }
 170       gp_ptr->fx=sf[0]+lf[0]+((1.0*rand()/RAND_MAX)-0.5)/rs;
 171       gp_ptr->fy=sf[1]+lf[1]+((1.0*rand()/RAND_MAX)-0.5)/rs;
 172       if(gp_ptr->status&FEST) {
 173         gp_ptr->fy+=(bf*exp(-y*1.0/100));
 174         if(gp_ptr->fx>0) {
 175           gp_ptr->fx-=gp_ptr->pp;
 176           if(gp_ptr->fx<0) gp_ptr->fx=0;
 177         }
 178         else {
 179           gp_ptr->fx+=gp_ptr->pp;
 180           if(gp_ptr->fx>0) gp_ptr->fx=0;
 181         }
 182         if(gp_ptr->fy>0) {
 183           gp_ptr->fy-=gp_ptr->pp;
 184           if(gp_ptr->fy<0) gp_ptr->fy=0;
 185         }
 186         else {
 187           gp_ptr->fy+=gp_ptr->pp;
 188           if(gp_ptr->fy>0) gp_ptr->fy=0;
 189         }
 190       }
 191       if(gp_ptr->y==0) gp_ptr->fy+=by; /* wenn rand by dazu, arschloch !!1 */
 192       if((mf<gp_ptr->fx)||(mf<gp_ptr->fy)) {
 193         mf=gp_ptr->fx>gp_ptr->fy?gp_ptr->fx:gp_ptr->fy;
 194         mx=gp_ptr->x;
 195         my=gp_ptr->y;
 196       }
 197       list.current=tmp;
 198     }
 199     while(list_next(&list)!=L_NO_NEXT_ELEMENT);
 200
 201     /* move vortex with highest force */
 202     gp.x=mx; gp.y=my;
 203     printf("step %d: move vortex %d %d,",i,mx,my);
 204     list_search_data(&list,&gp,2*sizeof(int));
 205     gp_ptr=(t_gp *)list.current->data;
 206     dx=0; dy=0;
 207     if(fabs(gp_ptr->fx)>fabs(gp_ptr->fy)) {
 208       if(gp_ptr->fx>0) dx=1;
 209       else dx=-1;
 210     }
 211     else {
 212       if(gp_ptr->fy>0) dy=1;
 213       else dy=-1;
 214     }
 215     printf(" with direction dx=%d dy=%d | force: %f %f\n",dx,dy,gp_ptr->fx,gp_ptr->fy);
 216     gp_ptr->x+=dx;
 217     gp_ptr->y+=dy;
 218     if(gp.y==0) list_add_element(&list,&gp,sizeof(t_gp));
 219
 220     /* plot every PLOT steps */
 221     if(i%PLOT==0) {
 222       memset(bmp.map,0,3*size*sizeof(unsigned char));
 223       list_reset(&list);
 224       do {
 225         gp_ptr=(t_gp *)list.current->data;
 226         if(gp_ptr->status&FEST)
 227           memset(bmp.map+gp_ptr->y*X_+gp_ptr->x,0xff,1);
 228         else
 229           memset(bmp.map+gp_ptr->y*X_+gp_ptr->x,0xff,3);
 230       } while(list_next(&list)!=L_NO_NEXT_ELEMENT);
 231       sprintf(string,"dfe_%d_of_%d.bmp",i,STEPS);
 232       strcpy(bmp.file,string);
 233       bmp_write_file(&bmp);
 234     }
 235
 236   }
 237
 238   close(fd);
 239
 240   return 1;
 241 }