From: hackbard Date: Tue, 28 Mar 2006 15:18:41 +0000 (+0000) Subject: nearly finished init stuff (probs with rand function!) X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=3ffe2a08e25fc091b6241885055450009267e2d8;p=physik%2Fposic.git nearly finished init stuff (probs with rand function!) --- diff --git a/Makefile b/Makefile index b6e31c0..ff3fbe3 100644 --- a/Makefile +++ b/Makefile @@ -1,12 +1,12 @@ CC=gcc CFLAGS=-Wall -OBJS=init/init.o visual/visual.o math/math.o moldyn.o +OBJS=init/init.o visual/visual.o math/math.o random/random.o moldyn.o all: moldyn.o posic posic: moldyn.o $(OBJS) - $(CC) $(CFLAGS) -o $@ $(OBJS) $(LIBS) posic.c + $(CC) $(CFLAGS) -lm -o $@ $(OBJS) $(LIBS) posic.c clean: rm -f *.o posic diff --git a/math/math.c b/math/math.c index 6d707e9..9cfaff6 100644 --- a/math/math.c +++ b/math/math.c @@ -7,6 +7,8 @@ #include +#include + #include "math.h" int v3_add(t_3dvec *sum,t_3dvec *a,t_3dvec *b) { @@ -64,3 +66,13 @@ int v3_cmp(t_3dvec *a,t_3dvec *b) { return(memcmp(a,b,sizeof(t_3dvec))); } +double v3_absolute_square(t_3dvec *a) { + + return(a->x*a->x+a->y*a->y+a->z*a->z); +} + +double v3_norm(t_3dvec *a) { + + return(sqrt(v3_absolute_square(a))); +} + diff --git a/math/math.h b/math/math.h index df8f6af..53e2114 100644 --- a/math/math.h +++ b/math/math.h @@ -25,5 +25,7 @@ int v3_zero(t_3dvec *vec); int v3_set(t_3dvec *vec,double *ptr); int v3_copy(t_3dvec *trg,t_3dvec *src); int v3_cmp(t_3dvec *a,t_3dvec *b); +double v3_absolute_square(t_3dvec *a); +double v3_norm(t_3dvec *a); #endif diff --git a/moldyn.c b/moldyn.c index e96cdd4..b507d65 100644 --- a/moldyn.c +++ b/moldyn.c @@ -9,9 +9,11 @@ #include #include +#include #include "math/math.h" #include "init/init.h" +#include "random/random.h" int create_lattice(unsigned char type,int element,double mass,double lc, @@ -63,38 +65,83 @@ int create_lattice(unsigned char type,int element,double mass,double lc, return ret; } -int thermal_init(t_atom *atom,int count,double t) { +int destroy_lattice(t_atom *atom) { + + if(atom) free(atom); + + return 0; +} + +int thermal_init(t_atom *atom,t_random *random,int count,double t) { /* * - gaussian distribution of velocities + * - zero total momentum * - velocity scaling (E = 3/2 N k T), E: kinetic energy */ int i; - double e,c,v; - - e=.0; + double v,sigma; + t_3dvec p_total,delta; + /* gaussian distribution of velocities */ + v3_zero(&p_total); for(i=0;i + * + */ + +#define _GNU_SOURCE +#include +#include +#include +#include +#include +#include +#include + +#include "random.h" + +int rand_init(t_random *random,char *randomfile,int logfd) { + + random->status=0; + + if(logfd) random->logfd=logfd; + + if(randomfile==NULL) { + random->fd=open("/dev/urandom",O_RDONLY); + if(random->fd<0) { + perror("open urandom"); + return -1; + } + random->status|=RAND_STAT_UDEV; + } else { + random->fd=open(randomfile,O_RDONLY); + if(random->fd<0) { + perror("open randomfile"); + return -1; + } + } + random->buffer=malloc(RAND_BUFSIZE*sizeof(unsigned int)); + if(random->buffer==NULL) { + perror("malloc random buffer"); + return -1; + } + random->b_ptr=random->buffer+RAND_BUFSIZE; + + dprintf(random->logfd,"[random] init successfull\n"); + + return 0; +} + + +int rand_close(t_random *random) { + + dprintf(random->logfd,"[random] shutdown\n"); + + if(random->buffer) free(random->buffer); + if(random->fd) close(random->fd); + if(random->logfd) close(random->logfd); + + return 0; +} + +unsigned int rand_get(t_random *random) { + + if(random->b_ptr==random->buffer+RAND_BUFSIZE) { + if(random->status&RAND_STAT_VERBOSE) + dprintf(random->logfd, + "[random] getting new random numbers\n"); + random->b_ptr=random->buffer; + if(!(random->status&RAND_STAT_UDEV)) { + lseek(random->fd,0,SEEK_SET); + dprintf(random->logfd, + "[random] warning, rereading random file\n"); + } + read(random->fd,random->b_ptr, + RAND_BUFSIZE*sizeof(unsigned int)); + if(random->status&RAND_STAT_VERBOSE) + dprintf(random->logfd, + "[random] got new random numbers\n"); + } + + return(*(random->b_ptr++)); +} + +double rand_get_double(t_random *random) { + + return(1.0*rand_get(random)/(double)(URAND_MAX+1)); +} + +double rand_get_double_linear(t_random *random,double a,double b) { + + return((-b+sqrt(b*b+2*a*(b+a/2)*rand_get_double(random)))/a); +} + +double rand_get_gauss(t_random *random) { + + double a,b,w; + + if(random->status&RAND_STAT_GAUSS) { + random->status&=~RAND_STAT_GAUSS; + return random->gauss; + } + + w=0; + while((w>=1.0)||(w==0.0)) { + a=-2.0*rand_get_double(random)+1.0; + b=-2.0*rand_get_double(random)+1.0; + w=a*a+b*b; + } + + w=sqrt(-2.0*log(w)/w); + + random->gauss=b*w; + random->status|=RAND_STAT_GAUSS; + + return(a*w); +} + +unsigned int rand_get_int_linear(t_random *random,unsigned int max, + double a,double b) { + + return((int)(rand_get_double_linear(random,a,b)*max)); +} + +unsigned int rand_get_reject(t_random *random,unsigned int max_x, + unsigned int max_y,unsigned int *graph) { + + unsigned int x,y; + unsigned char ok; + + ok=0; + while(!ok) { + x=(max_x*rand_get_double(random)); + y=(max_y*rand_get_double(random)); + if(y<=graph[x]) ok=1; + } + + return x; +} + diff --git a/random/random.h b/random/random.h new file mode 100644 index 0000000..8a21276 --- /dev/null +++ b/random/random.h @@ -0,0 +1,48 @@ +/* + * random.c - random functions header file + * + * author: Frank Zirkelbach + * + */ + +#ifndef RANDOM_H +#define RANDOM_H + +/* datatypes */ + +typedef struct s_random { + int fd; + int logfd; + unsigned char status; + unsigned int *buffer; + unsigned int *b_ptr; + double gauss; +} t_random; + +/* defines */ + +#define RAND_STAT_INITIALIZED 0x01 +#define RAND_STAT_VERBOSE 0x02 +#define RAND_STAT_UDEV 0x04 +#define RAND_STAT_GAUSS 0x08 + +#define RAND_BUFSIZE (1024) /* 16 mega byte */ + +#define URAND_MAX 0xffffffff +#define URAND_MAX_PLUS_ONE 0x100000000 + +/* function prototypes */ + +int rand_init(t_random *random,char *randomfile,int logfd); +int rand_close(t_random *random); +unsigned int rand_get(t_random *random); +double rand_get_double(t_random *random); +double rand_get_double_linear(t_random *random,double a,double b); +double rand_get_gauss(t_random *random); +unsigned int rand_get_int_linear(t_random *radnom,unsigned int max, + double a,double b); +unsigned int rand_get_reject(t_random *random,unsigned int max_x, + unsigned int max_y,unsigned int *graph); + + +#endif /* RANDOM_H */