/****************************************************************************/ /* UTIL: various useful routines and functions. */ /* */ /* Copyright (c) 1993 by Joshua E. Barnes, Honolulu, HI. */ /* It's free because it's yours. */ /****************************************************************************/ #include "stdinc.h" #include "real.h" #include "vectmath.h" #include #include #include #include #include #include #ifndef HZ # include # define HZ CLK_TCK #endif #define IM1 2147483563 #define IM2 2147483399 #define AM (1.0/IM1) #define IMM1 (IM1-1) #define IA1 40014 #define IA2 40692 #define IQ1 53668 #define IQ2 52774 #define IR1 12211 #define IR2 3791 #define NTAB 32 #define NDIV (1+IMM1/NTAB) #define EPSZU 1.2e-7 #define RNMX (1.0-EPSZU) /* * RAN2: long period random number generator. */ real ran2(long *idum) { int jj; long kk; static long idum2=123456789; static long iy=0; static long iv[NTAB]; real temp; if (*idum <= 0) { if (-(*idum) < 1) *idum=1; else *idum= -(*idum); idum2=(*idum); for (jj=NTAB+7;jj>=0;jj--) { kk=(*idum)/IQ1; *idum=IA1*(*idum-kk*IQ1)-kk*IR1; if (*idum < 0) *idum += IM1; if (jj < NTAB) iv[jj] = *idum; } iy=iv[0]; } kk=(*idum)/IQ1; *idum=IA1*(*idum-kk*IQ1)-kk*IR1; if (*idum < 0) *idum += IM1; kk=idum2/IQ2; idum2=IA2*(idum2-kk*IQ2)-kk*IR2; if (idum2 < 0) idum2 += IM2; jj=iy/NDIV; iy=iv[jj]-idum2; iv[jj]=*idum; if (iy < 1) iy += IMM1; if ((temp=AM*iy) > RNMX) return RNMX; else return temp; } /* * EXDEV: exponentially distributed, positive, random deviate of unit mean: exp(-aaa*r). */ real exdev(long *idum,real ll) { real ran2(long *idum); real dum,A,x1,rtild; A=1.0; x1=ran2(idum)*A; while (x1 < 1.0e-1) x1=ran2(idum)*A; rtild=-1.0/ll*log(x1); dum=rtild/(log(10.0)/ll); return dum; } /* * EXDEV2: exponentially distributed, positive, random deviate of unit mean: r*exp(-aaa*r). */ real exdev2(long *idum) { real ran2(long *idum); real dum,A,p,f,x1,x2,rtild,aaa; A=100.0; aaa=0.5; p=f=1.0; x2=100.0; while ( x2 >= p/f) { x1=ran2(idum)*A; if (x1 != 100.0) rtild=-100.0*log((100.0-x1)/100.0); f=10.0*exp(-rtild/10.0); p=rtild*exp(-aaa*rtild); x2=ran2(idum); } dum=rtild; return dum; } /* * HVELDEV: Maxwell Distribution of form v^2*exp(v^2/vr^2) */ real hveldev(long *idum,real vr) { real ran2(long *idum); real dum,A,p,f,x1,x2,rtild; A=1.0; p=f=1.0; x2=100.0; while ( x2 >= p/f) { x1=ran2(idum)*vr; if (x1 < (1.0e-4*vr)) x1=1.0e-4*vr; rtild=-vr*log(1.0-x1/vr); f=exp(-rtild/vr); p=4.0*PI*pow(1./(2.0*PI*vr*vr),1.5)*rtild*rtild*exp(-(rtild*rtild)/(2.0*vr*vr)); x2=ran2(idum); } dum=rtild; return dum; } /* * GASDEV: normally distributed deviate with zero mean and unit variance. */ real gasdev(long *idum) { real ran2(long *idum); static int iset; static real gset; real fac,rsq,v1,v2; if (iset == 0) { do { v1=2.0*ran2(idum)-1.0; v2=2.0*ran2(idum)-1.0; rsq=v1*v1+v2*v2; } while (rsq >= 1.0 || rsq == 0.0); fac=sqrt(-2.0*log(rsq)/rsq); gset=v1*fac; iset=1; return v2*fac; } else { iset=0; return gset; } } /* * MASSDEV: mass distribution with M^-1.5. */ real massdev(long *idum) { real ran2(long *idum); real dum; do dum=ran2(idum); while (dum == 0.0); return 4.0/(dum*dum)-4.0; } /* * HALODEV: 1/(r*(1+r)^2) distribution for dark halo */ real halodev(long *idum,real ofset) { real ran2(long *idum); real dum,A,p,f,x1,x2,rtild,rc; A=1.0/(ofset+1.0)-0.5; p=f=1.0; x2=10.0; while ( x2 >= p/f) { x1=ran2(idum)*A; rtild=1.0/(1.0/(ofset+1.0)-x1)-1.0; f=1.0/((rtild+1.0)*(rtild+1.0)); p=1.0/((rtild+0.1)*(5.0+rtild)*(5.0+rtild)); x2=ran2(idum); } dum=rtild; return dum; } /* * BULGEDEV: 1/(r*(1+r)^3) distribution for the bulge */ real bulgedev(long *idum,real ofset) { real ran2(long *idum); real dum,A,p,f,x1,x2,rtild,rc; A=1.0/(ofset+1.0)-0.5; p=f=1.0; x2=10.0; while ( x2 >= p/f) { x1=ran2(idum)*A; rtild=1.0/(1.0/(ofset+1.0)-x1)-1.0; f=1.0/((rtild+1.0)*(rtild+1.0)); p=1.0/((rtild+0.1)*(5.0+rtild)*(5.0+rtild)*(5.0+rtild)); x2=ran2(idum); } dum=rtild; return dum; } /* * quicksort subroutines (rhocmp in derivs1.c) */ int intcmp(int*a,int*b) { return(*a-*b); } void error(string, ...); extern float frand48(void); /* should be in math.h */ /* * XRANDOM: generate floating-point random number. */ real xrandom(real xl, real xh) { return (xl + (xh - xl) * drand48()); } /* * RSQR: compute x*x. */ real rsqr(real x) { return (x * x); } /* * DISTV: subtract vectors and return distance between. */ real distv(vector v, vector u) { real s, d; int n = NDIM; s = 0.0; while (--n >= 0) { d = (*v++) - (*u++); s += d * d; } return (rsqrt(s)); } /* * STREQ: test for equality of strings. */ bool streq(string a, string b) { return (strcmp(a, b) == 0); } /* * SCANOPT: scan string of the form "word1,word2,..." for match. Warning: * words must be separated by exactly one comma -- no spaces allowed! */ bool scanopt(string opt, string key) { char *op, *kp; op = (char *) opt; /* start scan of options */ while (*op != NULL) { /* loop over words in opt */ kp = (char *) key; /* start at front of key */ while ((*op != ',' ? *op : NULL) /* loop while this word */ == *kp) { /* ...matches text of key */ if (*kp++ == NULL) /* at end of key word? */ return (TRUE); /* keyword found */ op++; /* else advance ptrs */ } while (*op != NULL && *op++ != ',') /* loop till end of word, */ continue; /* passing "," at end */ } return (FALSE); /* keyword not found */ } /* * CPUTIME: compute CPU time in minutes. */ real cputime() { struct tms buffer; if (times(&buffer) == -1) error("times() call failed\n"); return (buffer.tms_utime / (60.0 * HZ)); } /* * ALLOCATE: memory allocation with error checking. */ void *allocate(nb) int nb; { void *mem; mem = (void *) calloc(nb, 1); /* calloc zeros memory */ if (mem == NULL) error("allocate: not enuf memory (%d bytes)\n", nb); return (mem); } /* * ERROR: print error message and exit. */ void error(string msg, ...) { va_list args; va_start(args, msg); vfprintf(stderr, msg, args); va_end(args); exit(-1); /* quit with error status */ } /* * EPRINTF: print error message, but don't exit. */ void eprintf(string msg, ...) { va_list args; va_start(args, msg); vfprintf(stderr, msg, args); va_end(args); }