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  * anmerkung frank: bisschen was ging so: ./dfe 1 5 2 1 1 10
  13  *                  ist auch schotter, aber mit richtigen parametern
  14  *                  wirds. viel glueck!!
  15  *
  16  * author: christian leirer, frank zirkelbach
  17  *
  18  * need api: cvs -d:pserver:anonymous@hackdaworld.dyndns.org:/my-code co api
  19  *
  20  */
  21
  22 #include "list.h"
  23 #include "bmp.h"
  24
  25 typedef struct s_gp {
  26   int x,y;
  27   unsigned char status;
  28 #define FLUESSIG 0
  29 #define FEST 1
  30   double fx,fy;
  31   double op;
  32   double pp;
  33 } t_gp;
  34
  35 #define RANGE 4
  36
  37 #define X_ 100
  38 #define Y_ 20
  39
  40 #define STEPS 1000
  41 #define PLOT 10
  42
  43 int main(int argc,char **argv) {
  44
  45   int x,y=0;
  46   int dx,dy;
  47   int mx=0,my=0;
  48   double lf[2];
  49   double sf[2];
  50   double mf;
  51   double lpp;
  52   int i;
  53   int size;
  54   t_gp gp,cgp;
  55   t_gp *gp_ptr;
  56   t_list list;
  57   t_list_element *tmp;
  58   t_bmp bmp;
  59   char string[128];
  60   int fd;
  61
  62   double by;
  63   double pp;
  64   double xi;
  65   double rs;
  66   double co;
  67   double bf;
  68
  69   /* parse argv */
  70
  71   if(argc!=7) {
  72     printf("usage: %s <b-feld> <p-pot> <xi> <rand-scale> <coulomb> <b-force>\n",
  73            argv[0]);
  74     return -1;
  75   }
  76
  77   by=atof(argv[1]);
  78   pp=atof(argv[2]);
  79   xi=atof(argv[3]);
  80   rs=atof(argv[4]);
  81   co=atof(argv[5]);
  82   bf=atof(argv[6]);
  83
  84   /* init */
  85
  86   fd=open("/dev/null",O_WRONLY);
  87
  88   list_init(&list,1);
  89
  90   size=X_*Y_;
  91   bmp_init(&bmp,1);
  92   bmp.width=X_;
  93   bmp.height=Y_;
  94   bmp.mode=WRITE;
  95   bmp_alloc_map(&bmp);
  96
  97   srand(time(NULL));
  98
  99   for(x=0;x<X_;x++) {
 100     memset(&gp,0,sizeof(t_gp));
 101     //gp.status=FEST;
 102     gp.x=x; gp.y=0;
 103     list_add_element(&list,&gp,sizeof(t_gp));
 104   }
 105
 106   /* mainloop */
 107   for(i=0;i<STEPS;i++) {
 108
 109     /* determine status of vortex, liquid or solid */
 110     list_reset(&list);
 111     do {
 112       lpp=0;
 113       gp_ptr=(t_gp *)list.current->data;
 114       tmp=list.current;
 115       for(dx=-RANGE;dx<=RANGE;dx++) {
 116         for(dy=-RANGE;dy<=RANGE;dy++) {
 117           if((dx!=0)||(dy!=0)) {
 118             x=(gp_ptr->x+X_+dx)%X_;
 119             y=gp_ptr->y+dy;
 120             if(y>=0&&y<Y_) {
 121               gp.x=x; gp.y=y;
 122               if(list_search_data(&list,&gp,2*sizeof(int))==L_SUCCESS)
 123                 lpp+=exp(-1.0*sqrt(x*x+y*y)/xi);
 124             }
 125           }
 126         }
 127       }
 128       gp_ptr->pp=pp-lpp;
 129       if(lpp>(0.8*pp)) gp_ptr->status=FLUESSIG;
 130       else gp_ptr->status=FEST;
 131       list.current=tmp;
 132     }
 133     while(list_next(&list)!=L_NO_NEXT_ELEMENT);
 134
 135     /* force on vortex */
 136     list_reset(&list);
 137     mf=0;
 138     do {
 139       lf[0]=0; lf[1]=0;
 140       sf[0]=0; sf[1]=0;
 141
 142       tmp=list.current;
 143       gp_ptr=(t_gp *)list.current->data;
 144       for(dx=-RANGE;dx<=RANGE;dx++) {
 145         for(dy=-RANGE;dy<=RANGE;dy++) {
 146           if((dx!=0)||(dy!=0)) {
 147             x=(gp_ptr->x+X_+dx)%X_;
 148             y=gp_ptr->y+dy;
 149             if(y>=0&&y<Y_) {
 150               gp.x=x; gp.y=y;
 151               if(list_search_data(&list,&gp,2*sizeof(int))==L_SUCCESS) {
 152                 memcpy(&gp,list.current->data,sizeof(t_gp));
 153                 if(gp.status&FEST) {
 154                   if(dy==0) sf[0]+=(-1.0/(dx*dx)*co*(dx/abs(dx)));
 155                   else if (dx==0) sf[1]+=(-1.0/(dy*dy)*co*(dy/abs(dy)));
 156                   else {
 157                     sf[0]+=(-1.0/(dx*dx+dy*dy)*co*1.0*dx/abs(dy));
 158                     sf[1]+=(-1.0/(dx*dx+dy*dy)*co*1.0*dy/abs(dx));
 159                   }
 160                 }
 161                 else {
 162                   if(dy==0) lf[0]+=(-1.0*(dx/abs(dx))*(by+co));
 163                   else if(dx==0) lf[1]+=(-1.0*(dy/abs(dy))*(by+co));
 164                   else{
 165                     lf[0]+=(-1.0*dx/abs(dy)*(by+co));
 166                     lf[1]+=(-1.0*dy/abs(dx)*(by+co));
 167                   }
 168                 }
 169               }
 170             }
 171           }
 172         }
 173       }
 174       gp_ptr->fx=sf[0]+lf[0]+((1.0*rand()/RAND_MAX)-0.5)/rs;
 175       gp_ptr->fy=sf[1]+lf[1]+((1.0*rand()/RAND_MAX)-0.5)/rs;
 176       if(gp_ptr->y==0) gp_ptr->fy+=by;
 177       if(gp_ptr->status&FEST) {
 178         gp_ptr->fy+=(bf*exp(-y*1.0/1000));
 179         if(gp_ptr->fx>0) {
 180           gp_ptr->fx-=gp_ptr->pp;
 181           if(gp_ptr->fx<0) gp_ptr->fx=0;
 182         }
 183         else {
 184           gp_ptr->fx+=gp_ptr->pp;
 185           if(gp_ptr->fx>0) gp_ptr->fx=0;
 186         }
 187         if(gp_ptr->fy>0) {
 188           gp_ptr->fy-=gp_ptr->pp;
 189           if(gp_ptr->fy<0) gp_ptr->fy=0;
 190         }
 191         else {
 192           gp_ptr->fy+=gp_ptr->pp;
 193           if(gp_ptr->fy>0) gp_ptr->fy=0;
 194         }
 195       }
 196       if((mf<gp_ptr->fx)||(mf<gp_ptr->fy)) {
 197         mf=gp_ptr->fx>gp_ptr->fy?gp_ptr->fx:gp_ptr->fy;
 198         mx=gp_ptr->x;
 199         my=gp_ptr->y;
 200       }
 201       list.current=tmp;
 202     }
 203     while(list_next(&list)!=L_NO_NEXT_ELEMENT);
 204
 205     /* move vortex with highest force */
 206     gp.x=mx; gp.y=my;
 207     printf("step %d: want to move vortex %d %d,",i,mx,my);
 208     list_search_data(&list,&gp,2*sizeof(int));
 209     gp_ptr=(t_gp *)list.current->data;
 210     dx=0; dy=0;
 211     if(fabs(gp_ptr->fx)>fabs(gp_ptr->fy)) {
 212       if(gp_ptr->fx>0) dx=1;
 213       else dx=-1;
 214     }
 215     else {
 216       if(gp_ptr->fy>0) dy=1;
 217       else dy=-1;
 218     }
 219     printf(" with direction dx=%d dy=%d | force: %f %f\n",dx,dy,gp_ptr->fx,gp_ptr->fy);
 220     cgp.x=gp_ptr->x+dx;
 221     cgp.y=gp_ptr->y+dy;
 222     if(list_search_data(&list,&cgp,2*sizeof(int))==L_SUCCESS) {
 223       printf("but there is allready a flussfaden! this parameter set suX!\n");
 224     }
 225     else if(gp_ptr->y==0&&dy==-1) {
 226       printf("but this wuld be stupid (out of target), skipped!\n");
 227     }
 228     else {
 229       printf("and i am doing it now! :)\n");
 230       gp_ptr->x+=dx;
 231       gp_ptr->y+=dy;
 232       if(gp.y==0) {
 233         printf("adding moved flussfaden (y=0)\n");
 234         list_add_element(&list,&gp,sizeof(t_gp));
 235       }
 236     }
 237     if(gp_ptr->y==Y_-1) {
 238       printf("flussfaden reached top of target, ending simulation!\n");
 239       return 1;
 240     }
 241
 242     /*
 243     printf("alle flussfaden:\n");
 244     list_reset(&list);
 245     do {
 246       gp_ptr=(t_gp *)list.current->data;
 247       printf("flussfaden bei %d %d\n",gp_ptr->x,gp_ptr->y);
 248     } while(list_next(&list)!=L_NO_NEXT_ELEMENT);
 249     */
 250
 251     /* plot every PLOT steps */
 252     if(i%PLOT==0) {
 253       memset(bmp.map,0,3*size*sizeof(unsigned char));
 254       list_reset(&list);
 255       do {
 256         gp_ptr=(t_gp *)list.current->data;
 257         if(gp_ptr->status&FEST)
 258           memset(bmp.map+gp_ptr->y*X_+gp_ptr->x,0xff,1);
 259         else
 260           memset(bmp.map+gp_ptr->y*X_+gp_ptr->x,0xff,3);
 261       } while(list_next(&list)!=L_NO_NEXT_ELEMENT);
 262       sprintf(string,"dfe_%d_of_%d.bmp",i,STEPS);
 263       strcpy(bmp.file,string);
 264       bmp_write_file(&bmp);
 265     }
 266
 267   }
 268
 269   close(fd);
 270
 271   return 1;
 272 }