guilhem_lavaux/vide_public
0
0
Fork 0
mirror of https://bitbucket.org/cosmicvoids/vide_public.git synced 2025-07-04 15:21:11 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Moved zobov into its private directory. New CMakeLists.txt to accomodate this

Browse source
This commit is contained in:
Guilhem Lavaux 2012-10-29 11:37:27 -04:00
parent 162fd7d9c4
commit a7a1e004e3
10 changed files with 0 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

81
zobov/findrtop.c Normal file
View file

@ -0,0 +1,81 @@
#include <math.h>
#include <stdlib.h>
/*------------------------------------------------------------------------------
Find nb elements of Real array a having the largest value.
Returns index iord of these elements, ordered so iord[0] corresponds
to element a[iord[0]] having the largest value.
If nb > na, then the last nb-na elements of iord are undefined.
Elements of a that are equal are left in their original order.
Courtesy Andrew Hamilton.
*/
void findrtop(double *a, int na, int *iord, int nb)
{
#undef ORDER
#define ORDER(ia, ja) a[ia] > a[ja] || (a[ia] == a[ja] && ia < ja)
int i, ia, ib, it, n;
n = (na <= nb)? na : nb;
if (n <= 0) return;
/* heap first n elements, so smallest element is at top of heap */
for (ib = (n >> 1); ib < n; ib++) {
iord[ib] = ib;
}
for (ia = (n >> 1) - 1; ia >= 0; ia--) {
i = ia;
for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) {
if (ib+1 < n) {
if (ORDER(iord[ib], iord[ib+1])) ib++;
}
if (ORDER(ia, iord[ib])) {
iord[i] = iord[ib];
i = ib;
} else {
break;
}
}
iord[i] = ia;
}
/* now compare rest of elements of array to heap */
for (ia = n; ia < na; ia++) {
/* if new element is greater than smallest, sift it into heap */
i = 0;
if (ORDER(ia, iord[i])) {
for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) {
if (ib+1 < n) {
if (ORDER(iord[ib], iord[ib+1])) ib++;
}
if (ORDER(ia, iord[ib])) {
iord[i] = iord[ib];
i = ib;
} else {
break;
}
}
iord[i] = ia;
}
}
/* unheap iord so largest element is at top */
for (ia = n - 1; ia > 0; ia--) {
it = iord[ia];
i = 0;
iord[ia] = iord[i];
for (ib = (i << 1) + 1; ib < ia; ib = (i << 1) + 1) {
if (ib+1 < ia) {
if (ORDER(iord[ib], iord[ib+1])) ib++;
}
if (ORDER(it, iord[ib])) {
iord[i] = iord[ib];
i = ib;
} else {
break;
}
}
iord[i] = it;
}
}

550
zobov/jozov.c Normal file
View file

@ -0,0 +1,550 @@
#include <assert.h>
/* jovoz.c by Mark Neyrinck */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define BIGFLT 1e30 /* Biggest possible floating-point number */
#define NLINKS 100 /* Number of possible links with the same rho_sl */
#define FF fflush(stdout)
typedef struct Particle {
float dens;
int nadj;
int ncnt;
int *adj;
} PARTICLE;
typedef struct Zone {
int core; /* Identity of peak particle */
int np; /* Number of particles in zone */
int npjoin; /* Number of particles in the joined void */
int nadj; /* Number of adjacent zones */
int nhl; /* Number of zones in final joined void */
float leak; /* Volume of last leak zone*/
int *adj; /* Each adjacent zone, with ... */
float *slv; /* Smallest Linking Volume */
float denscontrast; /* density contrast */
double vol; /* Total volume of all particles in the zone */
double voljoin; /* Total volume of all particles in the joined void */
} ZONE;
typedef struct ZoneT {
int nadj; /* Number of zones on border */
int *adj; /* Each adjacent zone, with ... */
float *slv; /* Smallest Linking Volume */
} ZONET;
void findrtop(double *a, int na, int *iord, int nb);
int main(int argc,char **argv) {
FILE *adj, *vol, *zon, *zon2, *txt;
PARTICLE *p;
ZONE *z;
ZONET *zt;
char *adjfile, *volfile, *zonfile, *zonfile2, *txtfile;
int i, j,k,l, h, h2,hl,n,np, np2, nzones, nhl, nhlcount, nhl2;
int *jumper, *jumped, *numinh;
int *zonenum, *zonelist, *zonelist2;
int link[NLINKS], link2, nl;
float lowvol, voltol, prob;
int q,q2;
int za, nin;
int testpart;
char already, interior, *inyet, *inyet2, added, beaten;
int *nm, **m;
float maxvol, minvol;
double *sorter, e1,maxdenscontrast;
int *iord;
e1 = exp(1.)-1.;
if (argc != 7) {
printf("Wrong number of arguments.\n");
printf("arg1: adjacency file\n");
printf("arg2: volume file\n");
printf("arg3: output file containing particles in each zone\n");
printf("arg4: output file containing zones in each void\n");
printf("arg5: output text file\n");
printf("arg6: Density threshold (0 for no threshold)\n\n");
exit(0);
}
adjfile = argv[1];
volfile = argv[2];
zonfile = argv[3];
zonfile2 = argv[4];
txtfile = argv[5];
if (sscanf(argv[6],"%f",&voltol) == 0) {
printf("Bad density threshold.\n");
exit(0);
}
if (voltol <= 0.) {
printf("Proceeding without a density threshold.\n");
voltol = 1e30;
}
adj = fopen(adjfile, "r");
if (adj == NULL) {
printf("Unable to open %s\n",adjfile);
exit(0);
}
fread(&np,1, sizeof(int),adj);
printf("adj: %d particles\n", np);
FF;
p = (PARTICLE *)malloc(np*sizeof(PARTICLE));
/* Adjacencies*/
for (i=0;i<np;i++) {
fread(&p[i].nadj,1,sizeof(int),adj);
/* The number of adjacencies per particle */
if (p[i].nadj > 0)
p[i].adj = (int *)malloc(p[i].nadj*sizeof(int));
else p[i].adj = 0;
p[i].ncnt = 0; /* Temporarily, it's an adj counter */
}
for (i=0;i<np;i++) {
fread(&nin,1,sizeof(int),adj);
if (nin > 0)
for (k=0;k<nin;k++) {
fread(&j,1,sizeof(int),adj);
if (j < np) {
/* Set both halves of the pair */
assert(i < j);
if (p[i].ncnt == p[i].nadj)
{
p[i].adj = (int *)realloc(p[i].adj, (p[i].nadj+1)*sizeof(int));
p[i].nadj++;
}
if (p[j].ncnt == p[j].nadj)
{
p[j].adj = (int *)realloc(p[j].adj, (p[j].nadj+1)*sizeof(int));
p[j].nadj++;
}
p[i].adj[p[i].ncnt] = j;
p[j].adj[p[j].ncnt] = i;
p[i].ncnt++; p[j].ncnt++;
} else {
printf("%d: adj = %d\n",i,j);
}
}
}
fclose(adj);
/* Check that we got all the pairs */
// adj = fopen(adjfile, "r");
// fread(&np,1, sizeof(int),adj);
for (i=0;i<np;i++) {
// fread(&nin,1,sizeof(int),adj); /* actually nadj */
if (p[i].ncnt != p[i].nadj) {
p[i].nadj = p[i].ncnt;
printf("We didn't get all of %d's adj's; %d != %d.\n",i,nin,p[i].nadj);
/*exit(0);*/
}
}
// fclose(adj);
/* Volumes */
vol = fopen(volfile, "r");
if (vol == NULL) {
printf("Unable to open volume file %s.\n\n",volfile);
exit(0);
}
fread(&np2,1, sizeof(int),adj);
if (np != np2) {
printf("Number of particles doesn't match! %d != %d\n",np,np2);
exit(0);
}
printf("%d particles\n", np);
FF;
for (i=0;i<np;i++) {
fread(&p[i].dens,1,sizeof(float),vol);
if ((p[i].dens < 1e-30) || (p[i].dens > 1e30)) {
printf("Whacked-out volume found, of particle %d: %f\n",i,p[i].dens);
p[i].dens = 1.;
}
p[i].dens = 1./p[i].dens; /* Get density from volume */
}
fclose(vol);
jumped = (int *)malloc(np*sizeof(int));
jumper = (int *)malloc(np*sizeof(int));
numinh = (int *)malloc(np*sizeof(int));
/* find jumper */
for (i = 0; i < np; i++) {
minvol = p[i].dens; jumper[i] = -1;
for (j=0; j< p[i].nadj; j++) {
if (p[p[i].adj[j]].dens < minvol) {
jumper[i] = p[i].adj[j];
minvol = p[jumper[i]].dens;
}
}
numinh[i] = 0;
}
printf("About to jump ...\n"); FF;
/* Jump */
for (i = 0; i < np; i++) {
jumped[i] = i;
while (jumper[jumped[i]] > -1)
jumped[i] = jumper[jumped[i]];
numinh[jumped[i]]++;
}
printf("Post-jump ...\n"); FF;
nzones = 0;
for (i = 0; i < np; i++)
if (numinh[i] > 0) nzones++;
printf("%d initial zones found\n", nzones);
z = (ZONE *)malloc(nzones*sizeof(ZONE));
if (z == NULL) {
printf("Unable to allocate z\n");
exit(0);
}
zt = (ZONET *)malloc(nzones*sizeof(ZONET));
if (zt == NULL) {
printf("Unable to allocate zt\n");
exit(0);
}
zonenum = (int *)malloc(np*sizeof(int));
if (zonenum == NULL) {
printf("Unable to allocate zonenum\n");
exit(0);
}
h = 0;
for (i = 0; i < np; i++)
if (numinh[i] > 0) {
z[h].core = i;
zonenum[i] = h;
h++;
} else {
zonenum[i] = -1;
}
/* Count border particles */
for (i = 0; i < np; i++)
for (j = 0; j < p[i].nadj; j++) {
testpart = p[i].adj[j];
if (jumped[i] != jumped[testpart])
zt[zonenum[jumped[i]]].nadj++;
}
for (h=0;h<nzones;h++) {
zt[h].adj = (int *)malloc(zt[h].nadj*sizeof(int));
if (zt[h].adj == NULL) {
printf("Unable to allocate %d adj's of zone %d\n",zt[h].nadj,h);
exit(0);
}
zt[h].slv = (float *)malloc(zt[h].nadj*sizeof(float));
if (zt[h].slv == NULL) {
printf("Unable to allocate %d slv's of zone %d\n",zt[h].nadj,h);
exit(0);
}
zt[h].nadj = 0;
}
/* Find "weakest links" */
for (i = 0; i < np; i++) {
h = zonenum[jumped[i]];
for (j = 0; j < p[i].nadj; j++) {
testpart = p[i].adj[j];
if (h != zonenum[jumped[testpart]]) {
if (p[testpart].dens > p[i].dens) {
/* there could be a weakest link through testpart */
already = 0;
for (za = 0; za < zt[h].nadj; za++)
if (zt[h].adj[za] == zonenum[jumped[testpart]]) {
already = 1;
if (p[testpart].dens < zt[h].slv[za]) {
zt[h].slv[za] = p[testpart].dens;
}
}
if (already == 0) {
zt[h].adj[zt[h].nadj] = zonenum[jumped[testpart]];
zt[h].slv[zt[h].nadj] = p[testpart].dens;
zt[h].nadj++;
}
} else { /* There could be a weakest link through i */
already = 0;
for (za = 0; za < zt[h].nadj; za++)
if (zt[h].adj[za] == zonenum[jumped[testpart]]) {
already = 1;
if (p[i].dens < zt[h].slv[za]) {
zt[h].slv[za] = p[i].dens;
}
}
if (already == 0) {
zt[h].adj[zt[h].nadj] = zonenum[jumped[testpart]];
zt[h].slv[zt[h].nadj] = p[i].dens;
zt[h].nadj++;
}
}
}
}
}
printf("Found zone adjacencies\n"); FF;
/* Free particle adjacencies */
for (i=0;i<np; i++) if (p[i].adj != 0) free(p[i].adj);
/* Use z instead of zt */
for (h=0;h<nzones;h++) {
/*printf("%d ",zt[h].nadj);*/
z[h].nadj = zt[h].nadj;
z[h].adj = (int *)malloc(zt[h].nadj*sizeof(int));
z[h].slv = (float *)malloc(zt[h].nadj*sizeof(float));
for (za = 0; za<zt[h].nadj; za++) {
z[h].adj[za] = zt[h].adj[za];
z[h].slv[za] = zt[h].slv[za];
}
free(zt[h].adj);
free(zt[h].slv);
z[h].np = numinh[z[h].core];
}
free(zt);
free(numinh);
m = (int **)malloc(nzones*sizeof(int *));
/* Not in the zone struct since it'll be freed up (contiguously, we hope)
soon */
nm = (int *)malloc(nzones*sizeof(int));
for (h=0; h<nzones; h++) {
m[h] = (int *)malloc(z[h].np*sizeof(int));
nm[h] = 0;
z[h].vol = 0.;
}
for (i=0; i<np; i++) {
h = zonenum[jumped[i]];
if (i == z[h].core) {
m[h][nm[h]] = m[h][0];
m[h][0] = i; /* Put the core particle at the top of the list */
} else {
m[h][nm[h]] = i;
}
z[h].vol += 1.0/(double)p[i].dens;
nm[h] ++;
}
free(nm);
zon = fopen(zonfile,"w");
if (zon == NULL) {
printf("Problem opening zonefile %s.\n\n",zonfile);
exit(0);
}
fwrite(&np,1,4,zon);
fwrite(&nzones,1,4,zon);
for (h=0; h<nzones; h++) {
fwrite(&(z[h].np),1,4,zon);
fwrite(m[h],z[h].np,4,zon);
free(m[h]);
}
free(m);
close(zon);
inyet = (char *)malloc(nzones*sizeof(char));
inyet2 = (char *)malloc(nzones*sizeof(char));
zonelist = (int *)malloc(nzones*sizeof(int));
zonelist2 = (int *)malloc(nzones*sizeof(int));
sorter = (double *)malloc((nzones+1)*sizeof(double));
for (h = 0; h< nzones; h++) {
inyet[h] = 0;
inyet2[h] = 0;
}
nhl = 0;
maxvol = 0.;
minvol = BIGFLT;
maxdenscontrast = 0.;
for(i=0;i<np; i++){
if (p[i].dens > maxvol) maxvol = p[i].dens;
if (p[i].dens < minvol) minvol = p[i].dens;
}
printf("Densities range from %e to %e.\n",minvol,maxvol);FF;
zon2 = fopen(zonfile2,"w");
if (zon2 == NULL) {
printf("Problem opening zonefile %s.\n\n",zonfile2);
exit(0);
}
fwrite(&nzones,1,4,zon2);
for (h = 0; h<nzones; h++) {
nhlcount = 0;
for (hl = 0; hl < nhl; hl++)
inyet[zonelist[hl]] = 0;
zonelist[0] = h;
inyet[h] = 1;
nhl = 1;
z[h].npjoin = z[h].np;
do {
/* Find the lowest-volume (highest-density) adjacency */
lowvol = BIGFLT; nl = 0; beaten = 0;
for (hl = 0; hl < nhl; hl++) {
h2 = zonelist[hl];
if (inyet[h2] == 1) { /* If it's not already identified as
an interior zone, with inyet=2 */
interior = 1;
for (za = 0; za < z[h2].nadj; za ++) {
if (inyet[z[h2].adj[za]] == 0) {
interior = 0;
if (z[h2].slv[za] == lowvol) {
link[nl] = z[h2].adj[za];
nl ++;
if (nl == NLINKS) {
printf("Too many links with the same rho_sl! Increase NLINKS from %d\n",nl);
exit(0);
}
}
if (z[h2].slv[za] < lowvol) {
lowvol = z[h2].slv[za];
link[0] = z[h2].adj[za];
nl = 1;
}
}
}
if (interior == 1) inyet[h2] = 2; /* No bordering exter. zones */
}
}
if (nl == 0) {
beaten = 1;
z[h].leak = maxvol;
continue;
}
if (lowvol > voltol) {
beaten = 1;
z[h].leak = lowvol;
continue;
}
for (l=0; l < nl; l++)
if (p[z[link[l]].core].dens < p[z[h].core].dens)
beaten = 1;
if (beaten == 1) {
z[h].leak = lowvol;
continue;
}
/* Add everything linked to the link(s) */
nhl2 = 0;
for (l=0; l < nl; l++) {
if (inyet2[link[l]] == 0) {
zonelist2[nhl2] = link[l];
inyet2[link[l]] = 1;
nhl2 ++;
added = 1;
while ((added == 1) && (beaten == 0)) {
added = 0;
for (hl = 0; (hl < nhl2) && (beaten == 0); hl++) {
h2 = zonelist2[hl];
if (inyet2[h2] == 1) {
interior = 1; /* Guilty until proven innocent */
for (za = 0; za < z[h2].nadj; za ++) {
link2 = z[h2].adj[za];
if ((inyet[link2]+inyet2[link2]) == 0) {
interior = 0;
if (z[h2].slv[za] <= lowvol) {
if (p[z[link2].core].dens < p[z[h].core].dens) {
beaten = 1;
break;
}
zonelist2[nhl2] = link2;
inyet2[link2] = 1;
nhl2++;
added = 1;
}
}
}
if (interior == 1) inyet2[h2] = 2;
}
}
}
}
}
for (hl = 0; hl < nhl2; hl++)
inyet2[zonelist2[hl]] = 0;
/* See if there's a beater */
if (beaten == 1) {
z[h].leak = lowvol;
} else {
for (h2 = 0; h2 < nhl2; h2++) {
zonelist[nhl] = zonelist2[h2];
inyet[zonelist2[h2]] = 1;
nhl++;
z[h].npjoin += z[zonelist2[h2]].np;
}
}
if (nhl/10000 > nhlcount) {
nhlcount = nhl/10000;
printf(" %d",nhl); FF;
}
} while((lowvol < BIGFLT) && (beaten == 0));
z[h].denscontrast = z[h].leak/p[z[h].core].dens;
if (z[h].denscontrast < 1.) z[h].denscontrast = 1.;
/* find biggest denscontrast */
if (z[h].denscontrast > maxdenscontrast) {
maxdenscontrast = (double)z[h].denscontrast;
}
/* Don't sort; want the core zone to be first */
if (nhlcount > 0) { /* Outputs the number of zones in large voids */
printf(" h%d:%d\n",h,nhl);
FF;
}
/* Calculate volume */
z[h].voljoin = 0.;
for (q = 0; q<nhl; q++) {
z[h].voljoin += z[zonelist[q]].vol;
}
z[h].nhl = nhl;
fwrite(&nhl,1,4,zon2);
fwrite(zonelist,nhl,4,zon2);
}
fclose(zon2);
printf("Maxdenscontrast = %f.\n",maxdenscontrast);
/* Assign sorter by probability (could use volume instead) */
for (h=0; h< nzones; h++)
sorter[h] = (double)z[h].denscontrast;
/* Text output file */
printf("about to sort ...\n");FF;
iord = (int *)malloc(nzones*sizeof(int));
findrtop(sorter, nzones, iord, nzones);
txt = fopen(txtfile,"w");
fprintf(txt,"%d particles, %d voids.\n", np, nzones);
fprintf(txt,"Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n");
for (h=0; h<nzones; h++) {
i = iord[h];
prob = exp(-5.12*(z[i].denscontrast-1.) - 0.8*pow(z[i].denscontrast-1.,2.8));
fprintf(txt,"%d %d %d %e %e %d %d %e %d %f %6.2e\n",
h+1, i, z[i].core, p[z[i].core].dens, z[i].vol, z[i].np, z[i].nhl, z[i].voljoin, z[i].npjoin, z[i].denscontrast, prob);
} /* h+1 to start from 1, not zero */
fclose(txt);
return(0);
}

609
zobov/jozov_persistent.c Normal file
View file

@ -0,0 +1,609 @@
#include <assert.h>
/* jovoz.c by Mark Neyrinck */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define BIGFLT 1e30 /* Biggest possible floating-point number */
#define NLINKS 100 /* Number of possible links with the same rho_sl */
#define FF fflush(stdout)
typedef struct Particle {
float dens;
int nadj;
int ncnt;
int *adj;
} PARTICLE;
typedef struct Zone {
int core; /* Identity of peak particle */
int np; /* Number of particles in zone */
int npjoin; /* Number of particles in the joined void */
int nadj; /* Number of adjacent zones */
int nhl; /* Number of zones in final joined void */
float leak; /* Volume of last leak zone*/
int *adj; /* Each adjacent zone, with ... */
float *slv; /* Smallest Linking Volume */
float denscontrast; /* density contrast */
double vol; /* Total volume of all particles in the zone */
double voljoin; /* Total volume of all particles in the joined void */
float smallest_saddle; /* Smallest saddle */
} ZONE;
typedef struct ZoneT {
int nadj; /* Number of zones on border */
int *adj; /* Each adjacent zone, with ... */
float *slv; /* Smallest Linking Volume */
float smallest_saddle; /* Smallest saddle point value */
} ZONET;
void findrtop(double *a, int na, int *iord, int nb);
int main(int argc,char **argv) {
FILE *adj, *vol, *zon, *zon2, *txt;
PARTICLE *p;
ZONE *z;
ZONET *zt;
char *adjfile, *volfile, *zonfile, *zonfile2, *txtfile;
int i, j,k,l, h, h2,hl,n,np, np2, nzones, nhl, nhlcount, nhl2;
int *jumper, *jumped, *numinh;
int *zonenum, *zonelist, *zonelist2;
int link[NLINKS], link2, nl;
float lowvol, voltol, prob;
int q,q2;
int za, nin;
int testpart;
char already, interior, *inyet, *inyet2, added, beaten;
int *nm, **m;
float maxvol, minvol;
double *sorter, e1,maxdenscontrast;
int *iord;
e1 = exp(1.)-1.;
if (argc != 7) {
printf("Wrong number of arguments.\n");
printf("arg1: adjacency file\n");
printf("arg2: volume file\n");
printf("arg3: output file containing particles in each zone\n");
printf("arg4: output file containing zones in each void\n");
printf("arg5: output text file\n");
printf("arg6: Density threshold (0 for no threshold)\n\n");
exit(0);
}
adjfile = argv[1];
volfile = argv[2];
zonfile = argv[3];
zonfile2 = argv[4];
txtfile = argv[5];
if (sscanf(argv[6],"%f",&voltol) == 0) {
printf("Bad density threshold.\n");
exit(0);
}
if (voltol <= 0.) {
printf("Proceeding without a density threshold.\n");
voltol = 1e30;
}
adj = fopen(adjfile, "r");
if (adj == NULL) {
printf("Unable to open %s\n",adjfile);
exit(0);
}
fread(&np,1, sizeof(int),adj);
printf("adj: %d particles\n", np);
FF;
p = (PARTICLE *)malloc(np*sizeof(PARTICLE));
/* Adjacencies*/
for (i=0;i<np;i++) {
fread(&p[i].nadj,1,sizeof(int),adj);
/* The number of adjacencies per particle */
if (p[i].nadj > 0)
p[i].adj = (int *)malloc(p[i].nadj*sizeof(int));
else p[i].adj = 0;
p[i].ncnt = 0; /* Temporarily, it's an adj counter */
}
for (i=0;i<np;i++) {
fread(&nin,1,sizeof(int),adj);
if (nin > 0)
for (k=0;k<nin;k++) {
fread(&j,1,sizeof(int),adj);
if (j < np) {
/* Set both halves of the pair */
assert(i < j);
if (p[i].ncnt == p[i].nadj)
{
p[i].adj = (int *)realloc(p[i].adj, (p[i].nadj+1)*sizeof(int));
p[i].nadj++;
}
if (p[j].ncnt == p[j].nadj)
{
p[j].adj = (int *)realloc(p[j].adj, (p[j].nadj+1)*sizeof(int));
p[j].nadj++;
}
p[i].adj[p[i].ncnt] = j;
p[j].adj[p[j].ncnt] = i;
p[i].ncnt++; p[j].ncnt++;
} else {
printf("%d: adj = %d\n",i,j);
}
}
}
fclose(adj);
/* Check that we got all the pairs */
// adj = fopen(adjfile, "r");
// fread(&np,1, sizeof(int),adj);
for (i=0;i<np;i++) {
// fread(&nin,1,sizeof(int),adj); /* actually nadj */
if (p[i].ncnt != p[i].nadj) {
p[i].nadj = p[i].ncnt;
printf("We didn't get all of %d's adj's; %d != %d.\n",i,nin,p[i].nadj);
/*exit(0);*/
}
}
// fclose(adj);
/* Volumes */
vol = fopen(volfile, "r");
if (vol == NULL) {
printf("Unable to open volume file %s.\n\n",volfile);
exit(0);
}
fread(&np2,1, sizeof(int),adj);
if (np != np2) {
printf("Number of particles doesn't match! %d != %d\n",np,np2);
exit(0);
}
printf("%d particles\n", np);
FF;
for (i=0;i<np;i++) {
fread(&p[i].dens,1,sizeof(float),vol);
if ((p[i].dens < 1e-30) || (p[i].dens > 1e30)) {
printf("Whacked-out volume found, of particle %d: %f\n",i,p[i].dens);
p[i].dens = 1.;
}
p[i].dens = 1./p[i].dens; /* Get density from volume */
}
fclose(vol);
jumped = (int *)malloc(np*sizeof(int));
jumper = (int *)malloc(np*sizeof(int));
numinh = (int *)malloc(np*sizeof(int));
/* find jumper */
for (i = 0; i < np; i++) {
minvol = p[i].dens; jumper[i] = -1;
for (j=0; j< p[i].nadj; j++) {
if (p[p[i].adj[j]].dens < minvol) {
jumper[i] = p[i].adj[j];
minvol = p[jumper[i]].dens;
}
}
numinh[i] = 0;
}
printf("About to jump ...\n"); FF;
/* Jump */
for (i = 0; i < np; i++) {
jumped[i] = i;
while (jumper[jumped[i]] > -1)
jumped[i] = jumper[jumped[i]];
numinh[jumped[i]]++;
}
printf("Post-jump ...\n"); FF;
nzones = 0;
for (i = 0; i < np; i++)
if (numinh[i] > 0) nzones++;
printf("%d initial zones found\n", nzones);
z = (ZONE *)malloc(nzones*sizeof(ZONE));
if (z == NULL) {
printf("Unable to allocate z\n");
exit(0);
}
zt = (ZONET *)malloc(nzones*sizeof(ZONET));
if (zt == NULL) {
printf("Unable to allocate zt\n");
exit(0);
}
zonenum = (int *)malloc(np*sizeof(int));
if (zonenum == NULL) {
printf("Unable to allocate zonenum\n");
exit(0);
}
h = 0;
for (i = 0; i < np; i++)
if (numinh[i] > 0) {
z[h].core = i;
zonenum[i] = h;
h++;
} else {
zonenum[i] = -1;
}
/* Count border particles */
for (i = 0; i < np; i++)
for (j = 0; j < p[i].nadj; j++) {
testpart = p[i].adj[j];
if (jumped[i] != jumped[testpart])
zt[zonenum[jumped[i]]].nadj++;
}
for (h=0;h<nzones;h++) {
zt[h].adj = (int *)malloc(zt[h].nadj*sizeof(int));
if (zt[h].adj == NULL) {
printf("Unable to allocate %d adj's of zone %d\n",zt[h].nadj,h);
exit(0);
}
zt[h].slv = (float *)malloc(zt[h].nadj*sizeof(float));
if (zt[h].slv == NULL) {
printf("Unable to allocate %d slv's of zone %d\n",zt[h].nadj,h);
exit(0);
}
zt[h].nadj = 0;
zt[h].smallest_saddle = BIGFLT;
}
/* Find "weakest links" */
for (i = 0; i < np; i++) {
h = zonenum[jumped[i]];
for (j = 0; j < p[i].nadj; j++) {
testpart = p[i].adj[j];
if (h != zonenum[jumped[testpart]]) {
if (p[testpart].dens > p[i].dens) {
/* there could be a weakest link through testpart */
already = 0;
for (za = 0; za < zt[h].nadj; za++)
if (zt[h].adj[za] == zonenum[jumped[testpart]]) {
already = 1;
if (p[testpart].dens < zt[h].slv[za]) {
zt[h].slv[za] = p[testpart].dens;
}
}
if (already == 0) {
zt[h].adj[zt[h].nadj] = zonenum[jumped[testpart]];
zt[h].slv[zt[h].nadj] = p[testpart].dens;
zt[h].nadj++;
}
} else { /* There could be a weakest link through i */
already = 0;
for (za = 0; za < zt[h].nadj; za++)
if (zt[h].adj[za] == zonenum[jumped[testpart]]) {
already = 1;
if (p[i].dens < zt[h].slv[za]) {
zt[h].slv[za] = p[i].dens;
}
}
if (already == 0) {
int q = zt[h].nadj;
zt[h].adj[q] = zonenum[jumped[testpart]];
zt[h].slv[q] = p[i].dens;
if (zt[h].smallest_saddle > zt[h].slv[q])
zt[h].smallest_saddle = zt[h].slv[q];
zt[h].nadj++;
}
}
}
}
}
printf("Found zone adjacencies\n"); FF;
/* Free particle adjacencies */
for (i=0;i<np; i++) if (p[i].adj != 0) free(p[i].adj);
/* Use z instead of zt */
for (h=0;h<nzones;h++) {
/*printf("%d ",zt[h].nadj);*/
z[h].nadj = zt[h].nadj;
z[h].adj = (int *)malloc(zt[h].nadj*sizeof(int));
z[h].slv = (float *)malloc(zt[h].nadj*sizeof(float));
for (za = 0; za<zt[h].nadj; za++) {
z[h].adj[za] = zt[h].adj[za];
z[h].slv[za] = zt[h].slv[za];
}
free(zt[h].adj);
free(zt[h].slv);
z[h].np = numinh[z[h].core];
z[h].smallest_saddle = zt[h].smallest_saddle;
}
free(zt);
free(numinh);
m = (int **)malloc(nzones*sizeof(int *));
/* Not in the zone struct since it'll be freed up (contiguously, we hope)
soon */
nm = (int *)malloc(nzones*sizeof(int));
for (h=0; h<nzones; h++) {
m[h] = (int *)malloc(z[h].np*sizeof(int));
nm[h] = 0;
z[h].vol = 0.;
}
for (i=0; i<np; i++) {
h = zonenum[jumped[i]];
if (i == z[h].core) {
m[h][nm[h]] = m[h][0];
m[h][0] = i; /* Put the core particle at the top of the list */
} else {
m[h][nm[h]] = i;
}
z[h].vol += 1.0/(double)p[i].dens;
nm[h] ++;
}
free(nm);
zon = fopen(zonfile,"w");
if (zon == NULL) {
printf("Problem opening zonefile %s.\n\n",zonfile);
exit(0);
}
fwrite(&np,1,4,zon);
fwrite(&nzones,1,4,zon);
for (h=0; h<nzones; h++) {
fwrite(&(z[h].np),1,4,zon);
fwrite(m[h],z[h].np,4,zon);
free(m[h]);
}
free(m);
close(zon);
inyet = (char *)malloc(nzones*sizeof(char));
inyet2 = (char *)malloc(nzones*sizeof(char));
zonelist = (int *)malloc(nzones*sizeof(int));
zonelist2 = (int *)malloc(nzones*sizeof(int));
sorter = (double *)malloc((nzones+1)*sizeof(double));
for (h = 0; h< nzones; h++) {
inyet[h] = 0;
inyet2[h] = 0;
}
nhl = 0;
maxvol = 0.;
minvol = BIGFLT;
maxdenscontrast = 0.;
for(i=0;i<np; i++){
if (p[i].dens > maxvol) maxvol = p[i].dens;
if (p[i].dens < minvol) minvol = p[i].dens;
}
printf("Densities range from %e to %e.\n",minvol,maxvol);FF;
zon2 = fopen(zonfile2,"w");
if (zon2 == NULL) {
printf("Problem opening zonefile %s.\n\n",zonfile2);
exit(0);
}
fwrite(&nzones,1,4,zon2);
#define CORE_DENSITY(a) p[z[a].core].dens
#define SADDLE_DENSITY(a) z[a].smallest_saddle
#define COMPARE_ZONES_BEATEN(a,b) (CORE_DENSITY(a) < CORE_DENSITY(b))
for (h = 0; h<nzones; h++) {
nhlcount = 0;
for (hl = 0; hl < nhl; hl++)
inyet[zonelist[hl]] = 0;
zonelist[0] = h;
inyet[h] = 1;
nhl = 1;
z[h].npjoin = z[h].np;
do {
/* Find the lowest-volume (highest-density) adjacency */
lowvol = BIGFLT; nl = 0; beaten = 0;
/* Look at all zones that have been investigated in all previous iterations. */
for (hl = 0; hl < nhl; hl++) {
h2 = zonelist[hl];
/* Is this zone been processed and deemed interior in this round ? */
if (inyet[h2] == 1) { /* If it's not already identified as
an interior zone, with inyet=2 */
/* Yes. Let's have a look at the adjacents... */
interior = 1;
for (za = 0; za < z[h2].nadj; za ++) {
/* This zones has not been deemed as interior */
if (inyet[z[h2].adj[za]] == 0) {
interior = 0;
/* Just at the threshold ? */
if (z[h2].slv[za] == lowvol) {
/* Yes. Queue it for adding the zones in bunch */
link[nl] = z[h2].adj[za];
nl ++;
if (nl == NLINKS) {
printf("Too many links with the same rho_sl! Increase NLINKS from %d\n",nl);
exit(0);
}
}
/* Just below the threshold ? */
if (z[h2].slv[za] < lowvol) {
/* Set the new threshold and reset the list because we may have beater */
lowvol = z[h2].slv[za];
link[0] = z[h2].adj[za];
nl = 1;
}
}
}
if (interior == 1)
inyet[h2] = 2; /* No bordering exter. zones */
}
}
if (nl == 0) {
beaten = 1;
z[h].leak = maxvol;
continue;
}
if (lowvol > voltol) {
beaten = 1;
z[h].leak = lowvol;
continue;
}
/* Check if among the found adjacent candidates there is a beater */
for (l=0; l < nl; l++)
if (COMPARE_ZONES_BEATEN(link[l], h))
beaten = 1;
if (beaten == 1) {
z[h].leak = lowvol;
continue;
}
/* Add everything linked to the link(s) */
nhl2 = 0;
for (l=0; l < nl; l++) {
if (inyet2[link[l]] == 0) {
zonelist2[nhl2] = link[l];
inyet2[link[l]] = 1;
nhl2 ++;
added = 1;
while ((added == 1) && (beaten == 0)) {
added = 0;
/* Loop over the newly found zones */
for (hl = 0; (hl < nhl2) && (beaten == 0); hl++) {
h2 = zonelist2[hl];
/* This hole has already been added to zonelist2 ? */
if (inyet2[h2] == 1) {
/* No, process it */
interior = 1; /* Guilty until proven innocent */
/* Take all adjacent zones of h2 */
for (za = 0; za < z[h2].nadj; za ++) {
link2 = z[h2].adj[za];
/* If this zone (link2) has not been already processed, process it */
if ((inyet[link2]+inyet2[link2]) == 0) {
interior = 0;
/* Does it have a low saddle point ? */
if (z[h2].slv[za] <= lowvol) {
/* Is this void beaten by the zone link2 ? */
if (COMPARE_ZONES_BEATEN(link2, h)) {
/* Yes. Stop here. */
beaten = 1;
break;
}
/* Insert it. */
zonelist2[nhl2] = link2;
inyet2[link2] = 1;
nhl2++;
/* Note that something was added so please continue */
added = 1;
}
}
}
if (interior == 1)
inyet2[h2] = 2;
}
}
}
}
}
for (hl = 0; hl < nhl2; hl++)
inyet2[zonelist2[hl]] = 0;
/* See if there's a beater */
if (beaten == 1) {
z[h].leak = lowvol;
} else {
for (h2 = 0; h2 < nhl2; h2++) {
zonelist[nhl] = zonelist2[h2];
inyet[zonelist2[h2]] = 1;
nhl++;
z[h].npjoin += z[zonelist2[h2]].np;
}
}
if (nhl/10000 > nhlcount) {
nhlcount = nhl/10000;
printf(" %d",nhl); FF;
}
} while((lowvol < BIGFLT) && (beaten == 0));
z[h].denscontrast = z[h].leak/p[z[h].core].dens;
if (z[h].denscontrast < 1.)
z[h].denscontrast = 1.;
/* find biggest denscontrast */
if (z[h].denscontrast > maxdenscontrast) {
maxdenscontrast = (double)z[h].denscontrast;
}
/* Don't sort; want the core zone to be first */
if (nhlcount > 0) { /* Outputs the number of zones in large voids */
printf(" h%d:%d\n",h,nhl);
FF;
}
/* Calculate volume */
z[h].voljoin = 0.;
for (q = 0; q<nhl; q++) {
z[h].voljoin += z[zonelist[q]].vol;
}
z[h].nhl = nhl;
fwrite(&nhl,1,4,zon2);
fwrite(zonelist,nhl,4,zon2);
}
fclose(zon2);
printf("Maxdenscontrast = %f.\n",maxdenscontrast);
/* Assign sorter by probability (could use volume instead) */
for (h=0; h< nzones; h++)
sorter[h] = (double)z[h].denscontrast;
/* Text output file */
printf("about to sort ...\n");FF;
iord = (int *)malloc(nzones*sizeof(int));
findrtop(sorter, nzones, iord, nzones);
txt = fopen(txtfile,"w");
fprintf(txt,"%d particles, %d voids.\n", np, nzones);
fprintf(txt,"Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n");
for (h=0; h<nzones; h++) {
i = iord[h];
prob = exp(-5.12*(z[i].denscontrast-1.) - 0.8*pow(z[i].denscontrast-1.,2.8));
fprintf(txt,"%d %d %d %e %e %d %d %e %d %f %6.2e\n",
h+1, i, z[i].core, p[z[i].core].dens, z[i].vol, z[i].np, z[i].nhl, z[i].voljoin, z[i].npjoin, z[i].denscontrast, prob);
} /* h+1 to start from 1, not zero */
fclose(txt);
return(0);
}

121
zobov/readfiles.c Normal file
View file

@ -0,0 +1,121 @@
#include <stdio.h>
#include <stdlib.h>
#define DL for (d=0;d<3;d++) /* Dimension loop */
#define BF 1e30
/* Positions */
/* Returns number of particles read */
int posread(char *posfile, float ***p, float fact) {
FILE *pos;
int np,dum,d,i;
float xmin,xmax,ymin,ymax,zmin,zmax;
float *ptemp;
pos = fopen(posfile, "r");
if (pos == NULL) {
printf("Unable to open position file %s\n\n",posfile);
exit(0);
}
/* Fortran77 4-byte headers and footers */
/* Delete "dum" statements if you don't need them */
/* Read number of particles */
fread(&dum,1,4,pos); fread(&np,1, sizeof(int),pos); fread(&dum,1,4,pos);
/* Allocate the arrays */
(*p) = (float **)malloc(np*sizeof(float *));
ptemp = (float *)malloc(np*sizeof(float));
printf("np = %d\n",np);
/* Fill the arrays */
fread(&dum,1,4,pos);
fread(ptemp,np,4,pos);
for (i=0; i<np; i++) {
(*p)[i] = (float *)malloc(3*sizeof(float));
if ((*p)[i] == NULL) {
printf("Unable to allocate particle array in readfiles!\n");
fflush(stdout);
exit(0);
}
(*p)[i][0] = ptemp[i];
}
fread(&dum,1,4,pos);
fread(&dum,1,4,pos);
fread(ptemp,np,4,pos);
for (i=0; i<np; i++) (*p)[i][1] = ptemp[i];
fread(&dum,1,4,pos);
fread(&dum,1,4,pos);
fread(ptemp,np,4,pos);
for (i=0; i<np; i++) (*p)[i][2] = ptemp[i];
fread(&dum,1,4,pos);
fclose(pos);
free(ptemp);
/* Get into physical units (Mpc/h) */
printf("%f\n",fact);fflush(stdout);
for (i=0; i<np; i++) DL (*p)[i][d] *= fact;
/* Test range -- can comment out */
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
for (i=0; i<np;i++) {
if ((*p)[i][0]<xmin) xmin = (*p)[i][0]; if ((*p)[i][0]>xmax) xmax = (*p)[i][0];
if ((*p)[i][1]<ymin) ymin = (*p)[i][1]; if ((*p)[i][1]>ymax) ymax = (*p)[i][1];
if ((*p)[i][2]<zmin) zmin = (*p)[i][2]; if ((*p)[i][2]>zmax) zmax = (*p)[i][2];
}
printf("np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np,xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout);
return(np);
}
/* Velocities */
/* Returns number of particles read */
/* Used in voboz, but not zobov, which doesn't use velocities */
int velread(char *velfile, float ***v, float fact) {
FILE *vel;
int np,dum,d,i;
float xmin,xmax,ymin,ymax,zmin,zmax;
vel = fopen(velfile, "r");
if (vel == NULL) {
printf("Unable to open velocity file %s\n\n",velfile);
exit(0);
}
/* Fortran77 4-byte headers and footers */
/* Delete "dum" statements if you don't need them */
/* Read number of particles */
fread(&dum,1,4,vel); fread(&np,1, sizeof(int),vel); fread(&dum,1,4,vel);
/* Allocate the arrays */
(*v) = (float **)malloc(3*sizeof(float*));
for (i=0;i<3;i++) (*v)[i] = (float *)malloc(np*sizeof(float));
/* Fill the arrays */
fread(&dum,1,4,vel); fread((*v)[0],np,4,vel); fread(&dum,1,4,vel);
fread(&dum,1,4,vel); fread((*v)[1],np,4,vel); fread(&dum,1,4,vel);
fread(&dum,1,4,vel); fread((*v)[2],np,4,vel); fread(&dum,1,4,vel);
fclose(vel);
/* Convert from code units into physical units (km/sec) */
for (i=0; i<np; i++) DL (*v)[d][i] *= fact;
/* Test range -- can comment out */
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
for (i=0; i<np;i++) {
if ((*v)[0][i] < xmin) xmin = (*v)[0][i]; if ((*v)[0][i] > xmax) xmax = (*v)[0][i];
if ((*v)[1][i] < ymin) ymin = (*v)[1][i]; if ((*v)[1][i] > ymax) ymax = (*v)[1][i];
if ((*v)[2][i] < zmin) zmin = (*v)[2][i]; if ((*v)[2][i] > zmax) zmax = (*v)[2][i];
}
printf("vx: %f,%f; vy: %f,%f; vz: %f,%f\n",xmin,xmax, ymin,ymax, zmin,zmax);
return(np);
}

8
zobov/voz.h Normal file
View file

@ -0,0 +1,8 @@
#define MAXVERVER 4000
#define NGUARD 42 /*Actually, the number of SPACES between guard points
in each dim */
typedef struct Partadj {
int nadj;
int *adj;
} PARTADJ;

337
zobov/voz1b1.c Normal file
View file

@ -0,0 +1,337 @@
#include "qhull_a.h"
#include "voz.h"
#define DL for (d=0;d<3;d++)
#define BF 1e30
#define MAX(a,b) ( ((a) < (b)) ? (a) : (b) )
int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs);
int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol);
int posread(char *posfile, float ***p, float fact);
int main(int argc, char *argv[]) {
int exitcode;
int i, j, np;
float **r;
coordT rtemp[3], *parts;
coordT deladjs[3*MAXVERVER], points[3*MAXVERVER];
pointT intpoints[3*MAXVERVER];
FILE *pos, *out;
char *posfile, outfile[80], *suffix;
PARTADJ *adjs;
float *vols;
float predict, xmin,xmax,ymin,ymax,zmin,zmax;
int *orig;
int isitinbuf;
char isitinmain, d;
int numdiv, nvp, nvpall, nvpbuf;
float width, width2, totwidth, totwidth2, bf, s, g;
float border, boxsize;
float c[3];
int b[3];
double totalvol;
if (argc != 9) {
printf("Wrong number of arguments.\n");
printf("arg1: position file\n");
printf("arg2: border size\n");
printf("arg3: boxsize\n");
printf("arg4: suffix\n");
printf("arg5: number of divisions\n");
printf("arg6-8: b[0-2]\n\n");
exit(0);
}
posfile = argv[1];
if (sscanf(argv[2],"%f",&border) != 1) {
printf("That's no border size; try again.\n");
exit(0);
}
if (sscanf(argv[3],"%f",&boxsize) != 1) {
printf("That's no boxsize; try again.\n");
exit(0);
}
suffix = argv[4];
if (sscanf(argv[5],"%d",&numdiv) != 1) {
printf("%s is no number of divisions; try again.\n",argv[5]);
exit(0);
}
if (numdiv == 1) {
printf("Only using one division; should only use for an isolated segment.\n");
}
if (numdiv < 1) {
printf("Cannot have a number of divisions less than 1. Resetting to 1.\n");
numdiv = 1;
}
if (sscanf(argv[6],"%d",&b[0]) != 1) {
printf("That's no b index; try again.\n");
exit(0);
}
if (sscanf(argv[7],"%d",&b[1]) != 1) {
printf("That's no b index; try again.\n");
exit(0);
}
if (sscanf(argv[8],"%d",&b[2]) != 1) {
printf("That's no b index; try again.\n");
exit(0);
}
/* Boxsize should be the range in r, yielding a range 0-1 */
np = posread(posfile,&r,1./boxsize);
printf("%d particles\n",np);fflush(stdout);
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
for (i=0; i<np;i++) {
if (r[i][0]<xmin) xmin = r[i][0]; if (r[i][0]>xmax) xmax = r[i][0];
if (r[i][1]<ymin) ymin = r[i][1]; if (r[i][1]>ymax) ymax = r[i][1];
if (r[i][2]<zmin) zmin = r[i][2]; if (r[i][2]>zmax) zmax = r[i][2];
}
printf("np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np,xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout);
width = 1./(float)numdiv;
width2 = 0.5*width;
if (border > 0.) bf = border;
else bf = 0.1;
/* In units of 0-1, the thickness of each subregion's buffer*/
totwidth = width+2.*bf;
totwidth2 = width2 + bf;
s = width/(float)NGUARD;
if ((bf*bf - 2.*s*s) < 0.) {
printf("bf = %f, s = %f.\n",bf,s);
printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.",
sqrt(2.)*width/bf);
exit(0);
}
g = (bf / 2.)*(1. + sqrt(1 - 2.*s*s/(bf*bf)));
printf("s = %f, bf = %f, g = %f.\n",s,bf,g);
fflush(stdout);
adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ));
if (adjs == NULL) {
printf("Unable to allocate adjs\n");
exit(0);
}
DL c[d] = ((float)b[d]+0.5)*width;
printf("c: %f,%f,%f\n",c[0],c[1],c[2]);
/* Assign temporary array*/
nvpbuf = 0; /* Number of particles to tesselate, including
buffer */
nvp = 0; /* Without the buffer */
for (i=0; i<np; i++) {
isitinbuf = 1;
isitinmain = 1;
DL {
rtemp[d] = (double)r[i][d] - (double)c[d];
if (rtemp[d] > 0.5) rtemp[d] --;
if (rtemp[d] < -0.5) rtemp[d] ++;
isitinbuf = isitinbuf && (fabs(rtemp[d]) < totwidth2);
isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
}
if (isitinbuf) nvpbuf++;
if (isitinmain) nvp++;
}
nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard
points */
parts = (coordT *)malloc(3*nvpbuf*sizeof(coordT));
orig = (int *)malloc(nvpbuf*sizeof(int));
if (parts == NULL) {
printf("Unable to allocate parts\n");
fflush(stdout);
}
if (orig == NULL) {
printf("Unable to allocate orig\n");
fflush(stdout);
}
nvp = 0; nvpall = 0; /* nvp = number of particles without buffer */
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
for (i=0; i<np; i++) {
isitinmain = 1;
DL {
rtemp[d] = r[i][d] - c[d];
if (rtemp[d] > 0.5) rtemp[d] --;
if (rtemp[d] < -0.5) rtemp[d] ++;
isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
}
if (isitinmain) {
parts[3*nvp] = rtemp[0];
parts[3*nvp+1] = rtemp[1];
parts[3*nvp+2] = rtemp[2];
orig[nvp] = i;
nvp++;
if (rtemp[0] < xmin) xmin = rtemp[0];
if (rtemp[0] > xmax) xmax = rtemp[0];
if (rtemp[1] < ymin) ymin = rtemp[1];
if (rtemp[1] > ymax) ymax = rtemp[1];
if (rtemp[2] < zmin) zmin = rtemp[2];
if (rtemp[2] > zmax) zmax = rtemp[2];
}
}
printf("nvp = %d\n",nvp);
printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
nvpbuf = nvp;
for (i=0; i<np; i++) {
isitinbuf = 1;
DL {
rtemp[d] = r[i][d] - c[d];
if (rtemp[d] > 0.5) rtemp[d] --;
if (rtemp[d] < -0.5) rtemp[d] ++;
isitinbuf = isitinbuf && (fabs(rtemp[d])<totwidth2);
}
if ((isitinbuf > 0) &&
((fabs(rtemp[0])>width2)||(fabs(rtemp[1])>width2)||(fabs(rtemp[2])>width2))) {
/*printf("%3.3f ",sqrt(rtemp[0]*rtemp[0] + rtemp[1]*rtemp[1] +
rtemp[2]*rtemp[2]));
printf("|%2.2f,%2.2f,%2.2f,%f,%f",r[i][0],r[i][1],r[i][2],width2,totwidth2);*/
parts[3*nvpbuf] = rtemp[0];
parts[3*nvpbuf+1] = rtemp[1];
parts[3*nvpbuf+2] = rtemp[2];
orig[nvpbuf] = i;
nvpbuf++;
if (rtemp[0] < xmin) xmin = rtemp[0];
if (rtemp[0] > xmax) xmax = rtemp[0];
if (rtemp[1] < ymin) ymin = rtemp[1];
if (rtemp[1] > ymax) ymax = rtemp[1];
if (rtemp[2] < zmin) zmin = rtemp[2];
if (rtemp[2] > zmax) zmax = rtemp[2];
}
}
printf("nvpbuf = %d\n",nvpbuf);
printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
nvpall = nvpbuf;
predict = pow(width+2.*bf,3)*(float)np;
printf("There should be ~ %f points; there are %d\n",predict,nvpbuf);
for (i=0;i<np;i++) free(r[i]);
free(r);
/* Add guard points */
for (i=0; i<NGUARD+1; i++) {
for (j=0; j<NGUARD+1; j++) {
/* Bottom */
parts[3*nvpall] = -width2 + (float)i * s;
parts[3*nvpall+1] = -width2 + (float)j * s;
parts[3*nvpall+2] = -width2 - g;
nvpall++;
/* Top */
parts[3*nvpall] = -width2 + (float)i * s;
parts[3*nvpall+1] = -width2 + (float)j * s;
parts[3*nvpall+2] = width2 + g;
nvpall++;
}
}
for (i=0; i<NGUARD+1; i++) { /* Don't want to overdo the corners*/
for (j=0; j<NGUARD+1; j++) {
parts[3*nvpall] = -width2 + (float)i * s;
parts[3*nvpall+1] = -width2 - g;
parts[3*nvpall+2] = -width2 + (float)j * s;
nvpall++;
parts[3*nvpall] = -width2 + (float)i * s;
parts[3*nvpall+1] = width2 + g;
parts[3*nvpall+2] = -width2 + (float)j * s;
nvpall++;
}
}
for (i=0; i<NGUARD+1; i++) {
for (j=0; j<NGUARD+1; j++) {
parts[3*nvpall] = -width2 - g;
parts[3*nvpall+1] = -width2 + (float)i * s;
parts[3*nvpall+2] = -width2 + (float)j * s;
nvpall++;
parts[3*nvpall] = width2 + g;
parts[3*nvpall+1] = -width2 + (float)i * s;
parts[3*nvpall+2] = -width2 + (float)j * s;
nvpall++;
}
}
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
for (i=nvpbuf;i<nvpall;i++) {
if (parts[3*i] < xmin) xmin = parts[3*i];
if (parts[3*i] > xmax) xmax = parts[3*i];
if (parts[3*i+1] < ymin) ymin = parts[3*i+1];
if (parts[3*i+1] > ymax) ymax = parts[3*i+1];
if (parts[3*i+2] < zmin) zmin = parts[3*i+2];
if (parts[3*i+2] > zmax) zmax = parts[3*i+2];
}
printf("Added guard points to total %d points (should be %d)\n",nvpall,
nvpbuf + 6*(NGUARD+1)*(NGUARD+1));
printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
/* Do tesselation*/
printf("File read. Tessellating ...\n"); fflush(stdout);
exitcode = delaunadj(parts, nvp, nvpbuf, nvpall, &adjs);
/* Now calculate volumes*/
printf("Now finding volumes ...\n"); fflush(stdout);
vols = (float *)malloc(nvp*sizeof(float));
for (i=0; i<nvp; i++) { /* Just the original particles
Assign adjacency coordinate array*/
/* Volumes */
for (j = 0; j < adjs[i].nadj; j++)
DL {
deladjs[3*j + d] = parts[3*adjs[i].adj[j]+d] - parts[3*i+d];
if (deladjs[3*j+d] < -0.5) deladjs[3*j+d]++;
if (deladjs[3*j+d] > 0.5) deladjs[3*j+d]--;
}
exitcode = vorvol(deladjs, points, intpoints, adjs[i].nadj, &(vols[i]));
vols[i] *= (float)np;
if (i % 1000 == 0)
printf("%d: %d, %f\n",i,adjs[i].nadj,vols[i]);
}
/* Get the adjacencies back to their original values */
for (i=0; i<nvp; i++)
for (j = 0; j < adjs[i].nadj; j++)
adjs[i].adj[j] = orig[adjs[i].adj[j]];
totalvol = 0.;
for (i=0;i<nvp; i++) {
totalvol += (double)vols[i];
}
printf("Average volume = %g\n",totalvol/(float)nvp);
/* Now the output!
First number of particles */
sprintf(outfile,"part.%s.%02d.%02d.%02d",suffix,b[0],b[1],b[2]);
printf("Output to %s\n\n",outfile);
out = fopen(outfile,"w");
if (out == NULL) {
printf("Unable to open %s\n",outfile);
exit(0);
}
fwrite(&np,1, sizeof(int),out);
fwrite(&nvp,1, sizeof(int),out);
printf("nvp = %d\n",nvp);
/* Tell us where the original particles were */
fwrite(orig,sizeof(int),nvp,out);
/* Volumes*/
fwrite(vols,sizeof(float),nvp,out);
/* Adjacencies */
for (i=0;i<nvp;i++) {
fwrite(&(adjs[i].nadj),1,sizeof(int),out);
if (adjs[i].nadj > 0)
fwrite(adjs[i].adj,adjs[i].nadj,sizeof(int),out);
else printf("0");
}
fclose(out);
return(0);
}

159
zobov/vozinit.c Normal file
View file

@ -0,0 +1,159 @@
#include "qhull_a.h"
#include "voz.h"
#define NUMCPU 2
#define DL for (d=0;d<3;d++)
#define BF 1e30
int posread(char *posfile, float ***p, float fact);
int main(int argc, char *argv[]) {
int i, np;
float **rfloat, rtemp[3];
FILE *pos, *scr;
char *posfile, scrfile[80], systemstr[90], *suffix;
float xmin,xmax,ymin,ymax,zmin,zmax;
int isitinbuf;
char isitinmain, d;
int numdiv;
int p;
int nvp, nvpall, nvpbuf, nvpmin, nvpmax, nvpbufmin, nvpbufmax; /* yes, the insurance */
float width, width2, totwidth, totwidth2, bf, s, g;
float border, boxsize;
float c[3];
int b[3];
if (argc != 6) {
printf("Wrong number of arguments.\n");
printf("arg1: position file\n");
printf("arg2: buffer size (default 0.1)\n");
printf("arg3: box size\n");
printf("arg4: number of divisions (default 2)\n");
printf("arg5: suffix describing this run\n\n");
exit(0);
}
posfile = argv[1];
suffix = argv[2];
if (sscanf(suffix,"%f",&border) != 1) {
printf("That's no border size; try again.\n");
exit(0);
}
suffix = argv[3];
if (sscanf(suffix,"%f",&boxsize) != 1) {
printf("That's no boxsize; try again.\n");
exit(0);
}
suffix = argv[4];
if (sscanf(suffix,"%d",&numdiv) != 1) {
printf("That's no number of divisions; try again.\n");
exit(0);
}
if (numdiv < 2) {
printf("Cannot have a number of divisions less than 2. Resetting to 2:\n");
numdiv = 2;
}
suffix = argv[5];
/* Read the position file */
np = posread(posfile,&rfloat,1./boxsize);
/* Boxsize should be the range in r, yielding a range 0-1 */
width = 1./(float)numdiv;
width2 = 0.5*width;
if (border > 0.) bf = border;
else bf = 0.1;
/* In units of 0-1, the thickness of each subregion's buffer*/
totwidth = width+2.*bf;
totwidth2 = width2 + bf;
s = width/(float)NGUARD;
if ((bf*bf - 2.*s*s) < 0.) {
printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.",
totwidth/sqrt(0.5*bf*bf));
printf("bf = %f\n",bf);
exit(0);
}
g = (bf / 2.)*(1. + sqrt(1 - 2.*s*s/(bf*bf)));
printf("s = %f, bf = %f, g = %f.\n",s,bf,g);
nvpmax = 0; nvpbufmax = 0; nvpmin = np; nvpbufmin = np;
for (b[0] = 0; b[0] < numdiv; b[0]++) {
c[0] = ((float)b[0]+0.5)*width;
for (b[1] = 0; b[1] < numdiv; b[1]++) {
c[1] = ((float)b[1]+0.5)*width;
for (b[2] = 0; b[2] < numdiv; b[2]++) {
c[2] = ((float)b[2]+0.5)*width;
nvp = 0; /* Number of particles excluding buffer */
nvpbuf = 0; /* Number of particles to tesselate, including
buffer */
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
for (i=0; i<np; i++) {
isitinbuf = 1; isitinmain = 1;
for (d=0; d<3; d++) {
rtemp[d] = rfloat[i][d] - c[d];
if (rtemp[d] > 0.5) rtemp[d] --;
if (rtemp[d] < -0.5) rtemp[d] ++;
isitinbuf = isitinbuf && (fabs(rtemp[d]) < totwidth2);
isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
}
if (isitinbuf) {
nvpbuf++;
}
if (isitinmain) {
nvp++;
if (rtemp[0] < xmin) xmin = rtemp[0];
if (rtemp[0] > xmax) xmax = rtemp[0];
if (rtemp[1] < ymin) ymin = rtemp[1];
if (rtemp[1] > ymax) ymax = rtemp[1];
if (rtemp[2] < zmin) zmin = rtemp[2];
if (rtemp[2] > zmax) zmax = rtemp[2];
}
}
if (nvp > nvpmax) nvpmax = nvp;
if (nvpbuf > nvpbufmax) nvpbufmax = nvpbuf;
if (nvp < nvpmin) nvpmin = nvp;
if (nvpbuf < nvpbufmin) nvpbufmin = nvpbuf;
printf("b=(%d,%d,%d), c=(%f,%f,%f), nvp=%d, nvpbuf=%d\n",
b[0],b[1],b[2],c[0],c[1],c[2],nvp,nvpbuf);
}
}
}
printf("Nvp range: %d,%d\n",nvpmin,nvpmax);
printf("Nvpbuf range: %d,%d\n",nvpbufmin,nvpbufmax);
/* Output script file */
sprintf(scrfile,"scr%s",suffix);
printf("Writing script file to %s.\n",scrfile);fflush(stdout);
scr = fopen(scrfile,"w");
if (scr == NULL) {
printf("Problem opening script file.\n");
fflush(stdout);
exit(0);
}
fprintf(scr,"#!/bin/bash -f\n");
p = 0;
for (b[0]=0;b[0]<numdiv; b[0]++) {
for (b[1] = 0; b[1] < numdiv; b[1]++) {
for (b[2] = 0; b[2] < numdiv; b[2]++) {
fprintf(scr,"./voz1b1 %s %f %f,%f,%f %s %d %d %d %d &\n",
posfile,border,boxsize,boxsize,boxsize,suffix,numdiv,b[0],b[1],b[2]);
p++;
if ((p == NUMCPU)) { fprintf(scr, "wait\n"); p = 0; }
}
}
}
fprintf(scr,"wait\n");
fprintf(scr,"./voztie %d %s\n",numdiv,suffix);
fclose(scr);
sprintf(systemstr,"chmod u+x %s",scrfile);
system(systemstr);
return(0);
}

182
zobov/voztie.c Normal file
View file

@ -0,0 +1,182 @@
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "voz.h"
int main(int argc, char *argv[]) {
FILE *part, *adj, *vol;
char partfile[80], *suffix, adjfile[80], volfile[80];
float *vols, volstemp;
PARTADJ *adjs;
int numdiv,np,np2,na;
int i,j,k,p,nout;
int nvp,npnotdone,nvpmax, nvpsum, *orig;
double avgnadj, avgvol;
if (argc != 3) {
printf("Wrong number of arguments.\n");
printf("arg1: number of divisions (default 2)\n");
printf("arg2: suffix describing this run\n\n");
exit(0);
}
if (sscanf(argv[1],"%d",&numdiv) != 1) {
printf("That's no number of divisions; try again.\n");
exit(0);
}
if (numdiv < 2) {
printf("Cannot have a number of divisions less than 2. Resetting to 2:\n");
numdiv = 2;
}
suffix = argv[2];
np = -1; nvpmax = -1; nvpsum = 0;
for (i = 0; i < numdiv; i++) {
for (j = 0; j < numdiv; j++) {
for (k = 0; k < numdiv; k++) {
sprintf(partfile,"part.%s.%02d.%02d.%02d",suffix,i,j,k);
part = fopen(partfile,"r");
if (part == NULL) {
printf("Unable to open file %s.\n\n",partfile);
exit(0);
}
fread(&np2,1,sizeof(int),part);
fread(&nvp,1,sizeof(int),part);
if (np == -1)
np = np2;
else
if (np2 != np) {
printf("Incompatible total particle numbers: %d,%d\n\n",np,np2);
exit(0);
}
if (nvp > nvpmax) nvpmax = nvp;
fclose(part);
}
}
}
printf("We have %d particles to tie together.\n",np); fflush(stdout);
printf("The maximum number of particles in a file is %d.\n",nvpmax);
adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ));
if (adjs == NULL) printf("Couldn't allocate adjs.\n");
vols = (float *)malloc(np*sizeof(float));
if (vols == NULL) printf("Couldn't allocate vols.\n");
orig = (int *)malloc(nvpmax*sizeof(int));
if (orig == NULL) printf("Couldn't allocate orig.\n");
if ((vols == NULL) || (orig == NULL) || (adjs == NULL)) {
printf("Not enough memory to allocate. Exiting.\n");
exit(0);
}
for (p=0;p<np;p++)
vols[p] = -1.;
for (i = 0; i < numdiv; i++) {
for (j = 0; j < numdiv; j++) {
for (k = 0; k < numdiv; k++) {
sprintf(partfile,"part.%s.%02d.%02d.%02d",suffix,i,j,k);
part = fopen(partfile,"r");
if (part == NULL) {
printf("Unable to open file %s.\n\n",partfile);
exit(0);
}
fread(&np2,1,sizeof(int),part);
fread(&nvp,1,sizeof(int),part);
/*printf("nvp = %d\n",nvp);fflush(stdout);*/
nvpsum += nvp;
fread(orig,nvp,sizeof(int),part);
for (p=0;p<nvp;p++) {
fread(&volstemp,1,sizeof(float),part);
if (vols[orig[p]] > -1.)
if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3) {
// printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n",
// orig[p],vols[orig[p]],volstemp);
// exit(0);
}
vols[orig[p]] = volstemp;
}
for (p=0;p<nvp;p++) {
fread(&na,1,sizeof(int),part);
if (na > 0) {
adjs[orig[p]].nadj = na;
adjs[orig[p]].adj = (int *)malloc(na*sizeof(int));
if (adjs[orig[p]].adj == NULL) {
printf("Couldn't allocate adjs[orig[%d]].adj.\n",p);
exit(0);
}
fread(adjs[orig[p]].adj,na,sizeof(int),part);
} else {
printf("0"); fflush(stdout);
}
}
fclose(part);
printf("%d ",k);
}
}
}
printf("\n");
npnotdone = 0; avgnadj = 0.; avgvol = 0.;
for (p=0;p<np;p++) {
if (vols[p] == -1.) npnotdone++;
avgnadj += (double)(adjs[p].nadj);
avgvol += (double)(vols[p]);
}
if (npnotdone > 0)
printf("%d particles not done!\n");
printf("%d particles done more than once.\n",nvpsum-np);
avgnadj /= (double)np;
avgvol /= (double)np;
printf("Average # adjacencies = %lf (%f for Poisson)\n",avgnadj,
48.*3.141593*3.141593/35.+2.);
printf("Average volume = %lf\n",avgvol);
/* Now the output! */
sprintf(adjfile,"adj%s.dat",suffix);
sprintf(volfile,"vol%s.dat",suffix);
printf("Outputting to %s, %s\n\n",adjfile,volfile);
adj = fopen(adjfile,"w");
if (adj == NULL) {
printf("Unable to open %s\n",adjfile);
exit(0);
}
fwrite(&np,1, sizeof(int),adj);
/* Adjacencies: first the numbers of adjacencies,
and the number we're actually going to write per particle */
for (i=0;i<np;i++)
fwrite(&adjs[i].nadj,1,sizeof(int),adj);
/* Now the lists of adjacencies (without double counting) */
for (i=0;i<np;i++)
if (adjs[i].nadj > 0) {
nout = 0;
for (j=0;j<adjs[i].nadj; j++) if (adjs[i].adj[j] > i) nout++;
fwrite(&nout,1,sizeof(int),adj);
for (j=0;j<adjs[i].nadj; j++)
if (adjs[i].adj[j] > i)
fwrite(&(adjs[i].adj[j]),1,sizeof(int),adj);
}
fclose(adj);
/* Volumes */
vol = fopen(volfile,"w");
if (vol == NULL) {
printf("Unable to open %s\n",volfile);
exit(0);
}
fwrite(&np,1, sizeof(int),vol);
fwrite(vols,sizeof(float),np,vol);
fclose(vol);
return(0);
}

241
zobov/vozutil.c Normal file
View file

@ -0,0 +1,241 @@
#include "qhull_a.h"
#include "voz.h"
#define FOREACHvertex2_(vertices) FOREACHsetelement_(vertexT, vertices2,vertex2)
/*char qh_version[] = "user_eg 3.1 2001/10/04"; */
int compar(const void * n1, const void * n2) {
int i1,i2;
i1 = *(int *)n1;
i2 = *(int *)n2;
return 2*(i1 > i2) - 1 + (i1 == i2);
}
/* Finds Delaunay adjacencies of a set of points */
int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs) {
int dim= 3; /* dimension of points */
boolT ismalloc= False; /* True if qhull should free points in qh_freeqhull() or reallocation */
char flags[250]; /* option flags for qhull, see qh_opt.htm */
FILE *outfile= stdout; /* output from qh_produce_output()
use NULL to skip qh_produce_output() */
FILE *errfile= stderr; /* error messages from qhull code */
int exitcode; /* 0 if no error from qhull */
int curlong, totlong; /* memory remaining after qh_memfreeshort */
int i, ver, count;
int numfacets, numsimplicial, numridges, totneighbors, numneighbors,
numcoplanars, numtricoplanars;
setT *vertices, *vertices2, *vertex_points, *coplanar_points;
vertexT *vertex, **vertexp;
vertexT *vertex2, **vertex2p;
int vertex_i, vertex_n;
facetT *facet, **facetp, *neighbor, **neighborp;
pointT *point, **pointp;
int numdiv;
PARTADJ adjst;
adjst.adj = (int *)malloc(MAXVERVER*sizeof(int));
if (adjst.adj == NULL) {
printf("Unable to allocate adjst.adj\n");
exit(0);
}
/* Delaunay triangulation*/
sprintf (flags, "qhull s d QJ");
exitcode= qh_new_qhull (dim, nvpall, points, ismalloc,
flags, outfile, errfile);
if (!exitcode) { /* if no error */
/* 'qh facet_list' contains the convex hull */
/* From qh_printvneighbors */
qh_countfacets(qh facet_list, NULL, 0, &numfacets, &numsimplicial,
&totneighbors, &numridges, &numcoplanars, &numtricoplanars);
qh_vertexneighbors();
vertices= qh_facetvertices (qh facet_list, NULL, 0);
vertex_points= qh_settemp (nvpall);
coplanar_points= qh_settemp (nvpall);
qh_setzero (vertex_points, 0, nvpall);
qh_setzero (coplanar_points, 0, nvpall);
FOREACHvertex_(vertices)
qh_point_add (vertex_points, vertex->point, vertex);
FORALLfacet_(qh facet_list) {
FOREACHpoint_(facet->coplanarset)
qh_point_add (coplanar_points, point, facet);
}
ver = 0;
FOREACHvertex_i_(vertex_points) {
(*adjs)[ver].nadj = 0;
if (vertex) {
/* Count the neighboring vertices, check that all are real
neighbors */
adjst.nadj = 0;
FOREACHneighbor_(vertex) {
if ((*adjs)[ver].nadj > -1) {
if (neighbor->visitid) {
vertices2 = neighbor->vertices;
FOREACHvertex2_(vertices2) {
if (ver != qh_pointid(vertex2->point)) {
adjst.adj[adjst.nadj] = (int)qh_pointid(vertex2->point);
adjst.nadj ++;
if (adjst.nadj > MAXVERVER) {
printf("Increase MAXVERVER to at least %d!\n",adjst.nadj);
exit(0);
}
}
}
} else {
printf(" %d",ver);
(*adjs)[ver].nadj = -1; /* There are unreal vertices here */
}
}
}
} else (*adjs)[ver].nadj = -2;
/* Enumerate the unique adjacencies*/
if (adjst.nadj >= 4) {
qsort((void *)adjst.adj, adjst.nadj, sizeof(int), &compar);
count = 1;
for (i=1; i<adjst.nadj; i++)
if (adjst.adj[i] != adjst.adj[i-1]) {
if (adjst.adj[i] >= nvpbuf) {
printf("Guard point encountered. Increase border and/or nguard.\n");
printf("P:(%f,%f,%f), G: (%f,%f,%f)\n",points[3*ver],points[3*ver+1],points[3*ver+2],
points[3*adjst.adj[i]],points[3*adjst.adj[i]+1],points[3*adjst.adj[i]+2]);
}
count++;
}
(*adjs)[ver].adj = (int *)malloc(count*sizeof(int));
if ((*adjs)[ver].adj == NULL) {
printf("Unable to allocate (*adjs)[ver].adj\n");
exit(0);
}
(*adjs)[ver].adj[0] = adjst.adj[0];
count = 1;
for (i=1; i<adjst.nadj; i++)
if (adjst.adj[i] != adjst.adj[i-1]) {
(*adjs)[ver].adj[count] = adjst.adj[i];
count++;
}
(*adjs)[ver].nadj = count;
} else {
printf("Number of adjacencies %d < 4, particle %d -> %d\n",adjst.nadj,ver,ver);
exit(0);
}
ver++;
if (ver == nvp) break;
}
qh_settempfree (&coplanar_points);
qh_settempfree (&vertex_points);
qh_settempfree (&vertices);
}
qh_freeqhull(!qh_ALL); /* free long memory */
qh_memfreeshort (&curlong, &totlong); /* free short memory and memory allocator */
if (curlong || totlong)
fprintf (errfile, "qhull internal warning (delaunadj): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
free(adjst.adj);
return exitcode;
}
/* Calculates the Voronoi volume from a set of Delaunay adjacencies */
int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol) {
int dim= 3; /* dimension of points */
boolT ismalloc= False; /* True if qhull should free points in qh_freeqhull() or reallocation */
char flags[250]; /* option flags for qhull, see qh_opt.htm */
FILE *outfile= NULL; /* output from qh_produce_output()
use NULL to skip qh_produce_output() */
FILE *errfile= stderr; /* error messages from qhull code */
int exitcode; /* 0 if no error from qhull */
facetT *facet; /* set by FORALLfacets */
int curlong, totlong; /* memory remaining after qh_memfreeshort */
coordT *point, *normp, *coordp, *feasiblep, *deladj;
int i, j, k;
boolT zerodiv;
float runsum;
char region;
/*coordT *points;
pointT *intpoints;*/
/* make point array from adjacency coordinates (add offset)*/
/*points = (coordT *)malloc(4*numpoints*sizeof(coordT));
if (points == NULL) {
printf("Unable to allocate points\n");
exit(0);
}*/
for (i=0; i<numpoints; i++) {
runsum = 0.;
deladj = deladjs + 3*i;
point = points + 4*i;
for (j=0;j<3;j++) {
runsum += deladj[j]*deladj[j];
point[j] = deladj[j];
}
point[3] = -0.5*runsum;
}
sprintf (flags, "qhull H0");
exitcode= qh_new_qhull (4, numpoints, points, ismalloc,
flags, outfile, errfile);
numpoints = 0;
if (!exitcode) { /* if no error */
FORALLfacets {
numpoints++;
}
/* Now we know how many points */
/*intpoints = (pointT *)malloc(dim*numpoints*sizeof(pointT));
if (intpoints == NULL) {
printf("Unable to allocate intpoints\n");
exit(0);
}*/
j = 0;
FORALLfacets {
if (!qh feasible_point) {
fprintf (stdout, "qhull input error (qh_printafacet): option 'Fp' needs qh feasible_point\n");
qh_errexit( qh_ERRinput, NULL, NULL);
}
point= coordp= intpoints + j*3;
j++;
normp= facet->normal;
feasiblep= qh feasible_point;
if (facet->offset < -qh MINdenom) {
for (k= qh hull_dim; k--; )
*(coordp++)= (*(normp++) / - facet->offset) + *(feasiblep++);
}else {
for (k= qh hull_dim; k--; ) {
*(coordp++)= qh_divzero (*(normp++), facet->offset, qh MINdenom_1,
&zerodiv) + *(feasiblep++);
if (zerodiv) {
qh_memfree (point, qh normal_size);
printf("LABELprintinfinite\n");
exit(0);
}
}
}
}
}
qh_freeqhull (!qh_ALL);
qh_memfreeshort (&curlong, &totlong);
/* Now we calculate the volume */
sprintf (flags, "qhull FA");
exitcode= qh_new_qhull (dim, numpoints, intpoints, ismalloc,
flags, outfile, errfile);
qh_getarea(qh facet_list);
*vol = qh totvol;
qh_freeqhull (!qh_ALL);
qh_memfreeshort (&curlong, &totlong);
if (curlong || totlong)
fprintf (errfile, "qhull internal warning (vorvol): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
/*free(points); free(intpoints);*/
return exitcode;
}

25
zobov/zobov_readme.txt Normal file
View file

@ -0,0 +1,25 @@
This is the readme file for ZOBOV (ZOnes Bordering On Voidness),
version 1.0 (Jan 5, 2008) a void-finding algorithm for sets of
particles (e.g. a cosmological N-body simulation) by Mark Neyrinck.
It is an inversion of the halo-finding algorithm VOBOZ (VOronoi BOund
Zones), which was also written by Mark Neyrinck, with the guidance of
Nick Gnedin and Andrew Hamilton.
The gz file which once contained this file should include the
following files:
Makefile jozov.c readme.txt voz1b1.c voztie.c
findrtop.c readfiles.c voz.h vozinit.c vozutil.c
For instructions on installing and running ZOBOV, please see
http://www.ifa.hawaii.edu/~neyrinck/zobov/. To use ZOBOV, you will
need to download the Qhull package, at http://www.qhull.org. For help
in understanding ZOBOV, please see the paper describing it
(arXiv:0712.3049), If these resources are insufficient, feel free to
email Mark.Neyrinck@colorado.edu with questions.
This is free software. It may be freely copied, modified, and
redistributed, as long as ZOBOV and its author are acknowledged.
There is no warranty or other guarantee of fitness for ZOBOV; it is
provided "as is." I welcome bug reports (but do not guarantee their
fixing), which should be sent to Mark.Neyrinck@colorado.edu.