#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#define RDFLAG O_RDONLY
#define RWFLAG O_RDWR
#define PI 3.141592653589
#define FOD 3.0
int main (int argc, char **argv)
{
int Q, nR;
int i,j;
i = 1;
while (i < argc)
{
if (strcmp(argv[i], "-Q")==0)
Q = atol(argv[++i]);
else if (strcmp(argv[i], "-nR")==0)
nR = atol(argv[++i]);
i++;
}
long N, NN;
N = 1<<Q;
NN = N*N;
short *B0;
short *B1;
B0 = malloc(sizeof(short)*N*N*N);
B1 = malloc(sizeof(short)*N*N*N);
for(i=0; i<N*N*N; i++){
B0[i] = 0;
B1[i] = 0;
}
float xF, yF, zF;
xF = FOD*N;
yF = 0.0;
zF = 0.0;
float nFF;
nFF = (xF*xF + yF*yF + zF*zF);
short px[8] = {-1,-1,-1,-1,1,1,1,1};
short py[8] = {-1,-1,1,1,-1,-1,1,1};
short pz[8] = {-1,1,-1,1,-1,1,-1,1};
short sx[8] = {0,0,0,0,2,2,2,2};
short sy[8] = {0,0,2,2,0,0,2,2};
short sz[8] = {0,2,0,2,0,2,0,2};
double a, aa;
aa = 3.0*1.00394*1.00394;
a = sqrt(aa);
short *qK;
long *qS;
double *qR;
qS = malloc(sizeof(long)*Q);
qR = malloc(sizeof(double)*Q);
qK = malloc(sizeof(short)*Q);
int q, coef;
for(q=0;q<Q;q++){
coef = 1<<(Q-q-1);
qS[q] = coef;
qR[q] = a*coef;
qK[q] = -1;
}
int *xseed;
int *yseed;
int *zseed;
xseed = malloc(sizeof(int)*Q);
yseed = malloc(sizeof(int)*Q);
zseed = malloc(sizeof(int)*Q);
xseed[0] = 0;
yseed[0] = 0;
zseed[0] = 0;
int nP;
nP = 3*N;
int *ind_0;
int *ind_1;
ind_0 = malloc(sizeof(int)*nP);
ind_1 = malloc(sizeof(int)*nP);
int ir0, ir1;
for(j=0;j<nP;j++){
ind_0[j] = -1;
ind_1[j] = -1;
}
float *tau_x;
float *tau_y;
float *tau_z;
tau_x = malloc(sizeof(float)*nR);
tau_y = malloc(sizeof(float)*nR);
tau_z = malloc(sizeof(float)*nR);
printf("before:\n");
float t;
for(j=0; j<nR; j++){
t = -0.5 + (j+0.5)/nR;
tau_x[j] = sqrt(1.0-t*t);
tau_y[j] = t;
tau_z[j] = 0.0;
printf("%d: %f,%f,%f\n", j, tau_x[j],tau_y[j], tau_z[j]);
}
short k;
int nx, ny, nz;
float r, rr, dd, tt, tx, ty, tz, taux, tauy, tauz, tauF, test;
printf("\n");
printf("in loop:\n");
for(j=0; j<nR; j++){
ir0 = 0;
ir1 = 0;
taux = tau_x[j];
tauy = tau_y[j];
tauz = tau_z[j];
printf("%d: %f,%f,%f\n", j, tau_x[j],tau_y[j], tau_z[j]);
tauF = taux*xF + tauy*yF + tauz*zF;
test = nFF-tauF*tauF;
if (test < NN){
q = 0;
while(1){
r = qR[q];
rr = r*r;
dd = rr;
k = qK[q];
while(k<7 && dd>=rr){
k+=1;
nx = xseed[q] + px[k]*qS[q];
ny = yseed[q] + py[k]*qS[q];
nz = zseed[q] + pz[k]*qS[q];
tx = nx-xF;
ty = ny-yF;
tz = nz-zF;
tt = tx*taux + ty*tauy + tz*tauz;
dd = (tx-tt*taux)*(tx-tt*taux) + (ty-tt*tauy)*(ty-tt*tauy) + (tz-tt*tauz)*(tz-tt*tauz);
}
if(dd < rr){
qK[q] = k;
q+=1;
xseed[q] = nx;
yseed[q] = ny;
zseed[q] = nz;
}else{
qK[q] = -1;
q-=1;
}
if (q == (Q-1)){
for(i=0;i<8;i++){
nx = xseed[q] + sx[i];
ny = yseed[q] + sy[i];
nz = zseed[q] + sz[i];
if(nx*nx+ny*ny+nz*nz < NN){
tx = nx-xF;
ty = ny-yF;
tz = nz-zF;
tt = tx*taux + ty*tauy + tz*tauz;
dd = (tx-tt*taux)*(tx-tt*taux) + (ty-tt*tauy)*(ty-tt*tauy) + (tz-tt*tauz)*(tz-tt*tauz);
if(dd<aa){
nx = (N+nx)>>1;
ny = (N+ny)>>1;
nz = (N+nz)>>1;
ind_0[ir0] = nz*NN+ny*N+nx;
ir0+=1;
}
}
}
for(i=0;i<8;i++){
nx = xseed[q] + px[i];
ny = yseed[q] + py[i];
nz = zseed[q] + pz[i];
if(nx*nx+ny*ny+nz*nz<NN){
tx = nx-xF;
ty = ny-yF;
tz = nz-zF;
tt = tx*taux + ty*tauy + tz*tauz;
dd = (tx-tt*taux)*(tx-tt*taux) + (ty-tt*tauy)*(ty-tt*tauy) + (tz-tt*tauz)*(tz-tt*tauz);
if(dd<aa){
nx = (N+nx)>>1;
ny = (N+ny)>>1;
nz = (N+nz)>>1;
ind_1[ir1] = nz*NN+ny*N+nx;
ir1+=1;
}
}
}
qK[q]=-1;
q-=1;
}//if q == Q-1
if(q == -1)
break;
}//while(1)
/*
for (i=0;i<ir0;i++)
B0[ind_0[i]] = 1;
for (i=0;i<ir1;i++)
B1[ind_1[i]] = 1;
int unit;
unit=open("image",RWFLAG|O_CREAT|O_TRUNC|O_APPEND,S_IREAD|S_IWRITE);
if(unit == -1){
printf("Cannot open file %s\n","image");
exit(1);
}
write(unit, B0, sizeof(short)*N*N*N);
write(unit, B1, sizeof(short)*N*N*N);
*/
}
}
printf("\n");
printf("after:\n");
for(j=0; j<nR; j++)
printf("%d: %f,%f,%f\n", j, tau_x[j],tau_y[j], tau_z[j]);
}