From e19e729588131bb483cd6a755b0578a740d8f59d Mon Sep 17 00:00:00 2001 From: jongheekang Date: Sun, 3 May 2015 21:41:55 -0400 Subject: [PATCH 1/6] Modify coordinate.c and h --- lenstools/extern/coordinates.c | 21 +++++++++++++++++++++ lenstools/extern/coordinates.h | 2 ++ 2 files changed, 23 insertions(+) diff --git a/lenstools/extern/coordinates.c b/lenstools/extern/coordinates.c index 63f4b871..cb2b75ab 100644 --- a/lenstools/extern/coordinates.c +++ b/lenstools/extern/coordinates.c @@ -23,4 +23,25 @@ long min_long(long i,long j){ else{ return j; } +} + +int signed_int(int i,int j){ + + if(i<=j){ + return i; + } + else{ + return -j; + } +} + +long signed_long(long i, long j){ + + if(i<=j){ + return i; + } + else{ + return -j; + } + } \ No newline at end of file diff --git a/lenstools/extern/coordinates.h b/lenstools/extern/coordinates.h index e5d21f6a..c179ddf8 100644 --- a/lenstools/extern/coordinates.h +++ b/lenstools/extern/coordinates.h @@ -4,6 +4,8 @@ int min_int(int,int); int max_int(int,int); long min_long(long,long); +int signed_int(int,int); +long signed_long(long,long); static inline long coordinate(long x,long y,long map_size){ return ((y+map_size) % map_size)*map_size + ((x+map_size) % map_size); From caca0e3397019ac9e2a59ea0d3fbc3c38ef7cec1 Mon Sep 17 00:00:00 2001 From: jongheekang Date: Sun, 3 May 2015 21:46:11 -0400 Subject: [PATCH 2/6] Modify grid.c and h and _gadget.c --- lenstools/extern/_gadget.c | 154 +++++++++++++++++++++++++++++++++++++ lenstools/extern/grid.c | 69 +++++++++++++++++ lenstools/extern/grid.h | 2 + 3 files changed, 225 insertions(+) diff --git a/lenstools/extern/_gadget.c b/lenstools/extern/_gadget.c index 89b35cab..6943d0c4 100644 --- a/lenstools/extern/_gadget.c +++ b/lenstools/extern/_gadget.c @@ -20,6 +20,8 @@ static char getPosVel_docstring[] = "Gets the positions or velocities of the par static char getID_docstring[] = "Gets the 4 byte int particles IDs from the Gadget2 snapshot"; static char write_docstring[] = "Writes the particles information to a Gadget snapshot, with a proper header"; static char grid3d_docstring[] = "Put the snapshot particles on a regularly spaced grid"; +static char mass_grid3d_docstring[] = "Put the snapshot particles on a regularly spaced grid and compute mass density"; +static char massfluc_grid3d_docstring[] = "Put the snapshot particles on a regularly spaced grid and compute mass density fluctuation"; static char adaptive_docstring[] = "Put the snapshot particles on a regularly spaced grid using adaptive smoothing"; //Method declarations @@ -28,6 +30,8 @@ static PyObject *_gadget_getPosVel(PyObject *self,PyObject *args); static PyObject *_gadget_getID(PyObject *self,PyObject *args); static PyObject * _gadget_write(PyObject *self,PyObject *args); static PyObject * _gadget_grid3d(PyObject *self,PyObject *args); +static PyObject * _gadget_mass_grid3d(PyObject *self,PyObject *args); +static PyObject * _gadget_massfluc_grid3d(PyObject *self,PyObject *args); static PyObject * _gadget_adaptive(PyObject *self,PyObject *args); //_gadget method definitions @@ -38,6 +42,8 @@ static PyMethodDef module_methods[] = { {"getID",_gadget_getID,METH_VARARGS,getID_docstring}, {"write",_gadget_write,METH_VARARGS,write_docstring}, {"grid3d",_gadget_grid3d,METH_VARARGS,grid3d_docstring}, + {"mass_grid3d",_gadget_mass_grid3d,METH_VARARGS,mass_grid3d_docstring}, + {"massfluc_grid3d",_gadget_massfluc_grid3d,METH_VARARGS,massfluc_grid3d_docstring}, {"adaptive",_gadget_adaptive,METH_VARARGS,adaptive_docstring}, {NULL,NULL,0,NULL} @@ -445,6 +451,154 @@ static PyObject *_gadget_grid3d(PyObject *self,PyObject *args){ } + +//mass_grid3d() implementation +static PyObject *_gadget_mass_grid3d(PyObject *self,PyObject *args){ + + PyObject *positions_obj,*bins_obj; + + //parse input tuple + if(!PyArg_ParseTuple(args,"OO",&positions_obj,&bins_obj)){ + return NULL; + } + + //interpret parsed objects as arrays + PyObject *positions_array = PyArray_FROM_OTF(positions_obj,NPY_FLOAT32,NPY_IN_ARRAY); + PyObject *binsX_array = PyArray_FROM_OTF(PyTuple_GET_ITEM(bins_obj,0),NPY_DOUBLE,NPY_IN_ARRAY); + PyObject *binsY_array = PyArray_FROM_OTF(PyTuple_GET_ITEM(bins_obj,1),NPY_DOUBLE,NPY_IN_ARRAY); + PyObject *binsZ_array = PyArray_FROM_OTF(PyTuple_GET_ITEM(bins_obj,2),NPY_DOUBLE,NPY_IN_ARRAY); + + + //check if anything failed + if(positions_array==NULL || binsX_array==NULL || binsY_array==NULL || binsZ_array==NULL){ + + Py_XDECREF(positions_array); + Py_XDECREF(binsX_array); + Py_XDECREF(binsY_array); + Py_XDECREF(binsZ_array); + + return NULL; + } + + //Get data pointers + float *positions_data = (float *)PyArray_DATA(positions_array); + double *binsX_data = (double *)PyArray_DATA(binsX_array); + double *binsY_data = (double *)PyArray_DATA(binsY_array); + double *binsZ_data = (double *)PyArray_DATA(binsZ_array); + + //Get info about the number of bins + int NumPart = (int)PyArray_DIM(positions_array,0); + int nx = (int)PyArray_DIM(binsX_array,0) - 1; + int ny = (int)PyArray_DIM(binsY_array,0) - 1; + int nz = (int)PyArray_DIM(binsZ_array,0) - 1; + + //Allocate the new array for the grid + npy_intp gridDims[] = {(npy_intp) nx,(npy_intp) ny,(npy_intp) nz}; + PyObject *grid_array = PyArray_ZEROS(3,gridDims,NPY_FLOAT32,0); + + if(grid_array==NULL){ + + Py_DECREF(positions_array); + Py_DECREF(binsX_array); + Py_DECREF(binsY_array); + Py_DECREF(binsZ_array); + + return NULL; + + } + + //Get a data pointer + float *grid_data = (float *)PyArray_DATA(grid_array); + + //Snap the particles on the grid + mass_grid3d(positions_data,NumPart,binsX_data[0],binsY_data[0],binsZ_data[0],binsX_data[1] - binsX_data[0],binsY_data[1] - binsY_data[0],binsZ_data[1] - binsZ_data[0],nx,ny,nz,grid_data); + + //return the grid + Py_DECREF(positions_array); + Py_DECREF(binsX_array); + Py_DECREF(binsY_array); + Py_DECREF(binsZ_array); + + //printf("NumPart = %d, nx = %d, ny = %d, nz = %d, bins_data[0] = %f,bins_data[1] = %f,bins_data[2] = %f \n", NumPart,nx,ny,nz,binsX_data[0], binsX_data[1],binsX_data[2]); + //printf("positions_data[0] = %f, positions_data[1] = %f, positions_data[2] = %f \n",positions_data[0],positions_data[1],positions_data[2]); + return grid_array; + +} + +//massfluc_grid3d() implementation +static PyObject *_gadget_massfluc_grid3d(PyObject *self,PyObject *args){ + + PyObject *positions_obj,*bins_obj; + double rho_0; + + //parse input tuple + if(!PyArg_ParseTuple(args,"OOd",&positions_obj,&bins_obj,&rho_0)){ + return NULL; + } + + //interpret parsed objects as arrays + PyObject *positions_array = PyArray_FROM_OTF(positions_obj,NPY_FLOAT32,NPY_IN_ARRAY); + PyObject *binsX_array = PyArray_FROM_OTF(PyTuple_GET_ITEM(bins_obj,0),NPY_DOUBLE,NPY_IN_ARRAY); + PyObject *binsY_array = PyArray_FROM_OTF(PyTuple_GET_ITEM(bins_obj,1),NPY_DOUBLE,NPY_IN_ARRAY); + PyObject *binsZ_array = PyArray_FROM_OTF(PyTuple_GET_ITEM(bins_obj,2),NPY_DOUBLE,NPY_IN_ARRAY); + + + //check if anything failed + if(positions_array==NULL || binsX_array==NULL || binsY_array==NULL || binsZ_array==NULL){ + + Py_XDECREF(positions_array); + Py_XDECREF(binsX_array); + Py_XDECREF(binsY_array); + Py_XDECREF(binsZ_array); + + return NULL; + } + + //Get data pointers + float *positions_data = (float *)PyArray_DATA(positions_array); + double *binsX_data = (double *)PyArray_DATA(binsX_array); + double *binsY_data = (double *)PyArray_DATA(binsY_array); + double *binsZ_data = (double *)PyArray_DATA(binsZ_array); + + //Get info about the number of bins + int NumPart = (int)PyArray_DIM(positions_array,0); + int nx = (int)PyArray_DIM(binsX_array,0) - 1; + int ny = (int)PyArray_DIM(binsY_array,0) - 1; + int nz = (int)PyArray_DIM(binsZ_array,0) - 1; + + //Allocate the new array for the grid + npy_intp gridDims[] = {(npy_intp) nx,(npy_intp) ny,(npy_intp) nz}; + PyObject *grid_array = PyArray_ZEROS(3,gridDims,NPY_FLOAT32,0); + + if(grid_array==NULL){ + + Py_DECREF(positions_array); + Py_DECREF(binsX_array); + Py_DECREF(binsY_array); + Py_DECREF(binsZ_array); + + return NULL; + + } + + //Get a data pointer + float *grid_data = (float *)PyArray_DATA(grid_array); + + //Snap the particles on the grid + massfluc_grid3d(positions_data,NumPart,binsX_data[0],binsY_data[0],binsZ_data[0],binsX_data[1] - binsX_data[0],binsY_data[1] - binsY_data[0],binsZ_data[1] - binsZ_data[0],nx,ny,nz,rho_0,grid_data); + + //return the grid + Py_DECREF(positions_array); + Py_DECREF(binsX_array); + Py_DECREF(binsY_array); + Py_DECREF(binsZ_array); + + //printf("NumPart = %d, nx = %d, ny = %d, nz = %d, bins_data[0] = %f,bins_data[1] = %f,bins_data[2] = %f \n", NumPart,nx,ny,nz,binsX_data[0], binsX_data[1],binsX_data[2]); + //printf("positions_data[0] = %f, positions_data[1] = %f, positions_data[2] = %f \n",positions_data[0],positions_data[1],positions_data[2]); + return grid_array; + +} + //adaptive() implementation static PyObject * _gadget_adaptive(PyObject *self,PyObject *args){ diff --git a/lenstools/extern/grid.c b/lenstools/extern/grid.c index 580e363d..fe0ccc07 100644 --- a/lenstools/extern/grid.c +++ b/lenstools/extern/grid.c @@ -34,6 +34,75 @@ int grid3d(float *positions,int Npart,double leftX,double leftY,double leftZ,dou } +//Snap particles on a 3d regularly spaced grid and compute the numberdensity (equal to mass density) within each grid +int mass_grid3d(float *positions,int Npart,double leftX,double leftY,double leftZ,double sizeX,double sizeY,double sizeZ,int nx,int ny,int nz,float *grid){ + + int n; + double i,j,k; + double volume_grid = sizeX * sizeY * sizeZ; + + //Cycle through the particles and for each one compute the position on the grid + for(n=0;n=0 && i=0 && j=0 && k=0 && i=0 && j=0 && k int grid3d(float *positions,int Npart,double leftX,double leftY,double leftZ,double sizeX,double sizeY,double sizeZ,int nx,int ny,int nz,float *grid); +int mass_grid3d(float *positions,int Npart,double leftX,double leftY,double leftZ,double sizeX,double sizeY,double sizeZ,int nx,int ny,int nz,float *grid); +int massfluc_grid3d(float *positions,int Npart,double leftX,double leftY,double leftZ,double sizeX,double sizeY,double sizeZ,int nx,int ny,int nz,double rho_0,float *grid); int adaptiveSmoothing(int NumPart,float *positions,double *rp,double *binning0, double *binning1,double center,int direction0,int direction1,int normal,int size0,int size1,int projectAll,double *lensingPlane); static inline double quadraticKernel(double d,double rp){ return (1.0/pow(rp,2)) * pow(1.0 - d/pow(rp,2),2); From db8dbfe52e2916bc3f6e2e910504e904b4d11de7 Mon Sep 17 00:00:00 2001 From: jongheekang Date: Sun, 3 May 2015 21:53:40 -0400 Subject: [PATCH 3/6] Modify azimuth.c and h and _topology.c --- lenstools/extern/_topology.c | 248 +++++++++++ lenstools/extern/azimuth.c | 816 +++++++++++++++++++++++++++++++++++ lenstools/extern/azimuth.h | 2 + 3 files changed, 1066 insertions(+) diff --git a/lenstools/extern/_topology.c b/lenstools/extern/_topology.c index 8155f3e2..02f7492d 100644 --- a/lenstools/extern/_topology.c +++ b/lenstools/extern/_topology.c @@ -25,6 +25,8 @@ static char hessian_docstring[] = "Compute the hessian of a 2D image"; static char minkowski_docstring[] = "Measure the three Minkowski functionals of a 2D image"; static char rfft2_azimuthal_docstring[] = "Measure azimuthal average of Fourier transforms of 2D image"; static char rfft3_azimuthal_docstring[] = "Measure azimuthal average of Fourier transforms of 3D scalar field"; +static char rfft3_azimuthal_bispectrum_docstring[] = "Measure angular averaged bispectrum of Fourier transforms of 3D scalar field"; +static char rfft3_bispectrum_docstring[] = "Measure full bispectrum of Fourier transforms of 3D scalar field"; //method declarations static PyObject *_topology_peakCount(PyObject *self,PyObject *args); @@ -34,6 +36,8 @@ static PyObject *_topology_hessian(PyObject *self,PyObject *args); static PyObject *_topology_minkowski(PyObject *self,PyObject *args); static PyObject *_topology_rfft2_azimuthal(PyObject *self,PyObject *args); static PyObject *_topology_rfft3_azimuthal(PyObject *self,PyObject *args); +static PyObject *_topology_rfft3_azimuthal_bispectrum(PyObject *self,PyObject *args); +static PyObject *_topology_rfft3_bispectrum(PyObject *self,PyObject *args); //_topology method definitions @@ -46,6 +50,8 @@ static PyMethodDef module_methods[] = { {"minkowski",_topology_minkowski,METH_VARARGS,minkowski_docstring}, {"rfft2_azimuthal",_topology_rfft2_azimuthal,METH_VARARGS,rfft2_azimuthal_docstring}, {"rfft3_azimuthal",_topology_rfft3_azimuthal,METH_VARARGS,rfft3_azimuthal_docstring}, + {"rfft3_azimuthal_bispectrum",_topology_rfft3_azimuthal_bispectrum,METH_VARARGS,rfft3_azimuthal_bispectrum_docstring}, + {"rfft3_bispectrum",_topology_rfft3_bispectrum,METH_VARARGS,rfft3_bispectrum_docstring}, {NULL,NULL,0,NULL} } ; @@ -972,3 +978,245 @@ static PyObject *_topology_rfft3_azimuthal(PyObject *self,PyObject *args){ return output; } +// rfft3_azimuthal_bispectrum() implementation +static PyObject *_topology_rfft3_azimuthal_bispectrum(PyObject *self,PyObject *args){ + + /*These are the inputs: the Fourier transforms of the two maps, the size of the pixel in k space, the k bin extremes at which calculate the azimuthal averages*/ + PyObject *ft_map1_obj,*ft_map2_obj,*ft_map3_obj, *ksmall_obj, *klarge_obj; + double kpixX,kpixY,kpixZ; + int size_small_x, size_small_y, size_small_z; + + /*Parse input tuple*/ + if(!PyArg_ParseTuple(args,"OOOdddiiiOO",&ft_map1_obj,&ft_map2_obj,&ft_map3_obj, &kpixX,&kpixY,&kpixZ,&size_small_x,&size_small_y,&size_small_z, &ksmall_obj, &klarge_obj)){ + return NULL; + } + + /*Interpret the parsed objects as numpy arrays*/ + PyObject *ft_map1_array = PyArray_FROM_OTF(ft_map1_obj,NPY_COMPLEX128,NPY_IN_ARRAY); + PyObject *ft_map2_array = PyArray_FROM_OTF(ft_map2_obj,NPY_COMPLEX128,NPY_IN_ARRAY); + PyObject *ft_map3_array = PyArray_FROM_OTF(ft_map3_obj,NPY_COMPLEX128,NPY_IN_ARRAY); + PyObject *ksmall_array = PyArray_FROM_OTF(ksmall_obj,NPY_DOUBLE,NPY_IN_ARRAY); + PyObject *klarge_array = PyArray_FROM_OTF(klarge_obj,NPY_DOUBLE,NPY_IN_ARRAY); + + /*Check if anything failed*/ + if(ft_map1_array==NULL || ksmall_array==NULL || ft_map2_array==NULL || ft_map3_array==NULL || klarge_array==NULL){ + + Py_XDECREF(ft_map1_array); + Py_XDECREF(ft_map2_array); + Py_XDECREF(ft_map3_array); + Py_XDECREF(ksmall_array); + Py_XDECREF(klarge_array); + + return NULL; + } + + /*Get the size of the map fourier transform*/ + int Nside_x = (int)PyArray_DIM(ft_map1_array,0); + int Nside_y = (int)PyArray_DIM(ft_map1_array,1); + int Nside_z = (int)PyArray_DIM(ft_map1_array,2); + + /*Get the number of k bin edges*/ + int Nlarge = (int)PyArray_DIM(klarge_array,0); + int Nsmall = (int)PyArray_DIM(ksmall_array,0); + + /*Build the arrays that will contain the output:bispectrum and hits*/ + npy_intp dims[] = {(npy_intp) Nsmall - 1,(npy_intp) Nlarge-1}; + PyObject *bispectrum_array = PyArray_ZEROS(2,dims,NPY_DOUBLE,0); + PyObject *hits_array = PyArray_ZEROS(2,dims,NPY_LONG,0); + + /*Check for failure*/ + if(bispectrum_array==NULL || hits_array==NULL){ + + Py_DECREF(ft_map1_array); + Py_DECREF(ft_map2_array); + Py_DECREF(ft_map3_array); + Py_DECREF(ksmall_array); + Py_DECREF(klarge_array); + Py_XDECREF(bispectrum_array); + Py_XDECREF(hits_array); + + return NULL; + } + + /*Call the C backend azimuthal average function*/ + if(azimuthal_bispectrum_rfft3((double _Complex *)PyArray_DATA(ft_map1_array),(double _Complex *)PyArray_DATA(ft_map2_array),(double _Complex *)PyArray_DATA(ft_map3_array),Nside_x,Nside_y,Nside_z,size_small_x,size_small_y,size_small_z,kpixX,kpixY,kpixZ,Nsmall,Nlarge,(double *)PyArray_DATA(ksmall_array),(double *)PyArray_DATA(klarge_array),(double *)PyArray_DATA(bispectrum_array),(long *)PyArray_DATA(hits_array))){ + + Py_DECREF(ft_map1_array); + Py_DECREF(ft_map2_array); + Py_DECREF(ft_map3_array); + Py_DECREF(ksmall_array); + Py_DECREF(klarge_array); + Py_XDECREF(bispectrum_array); + Py_XDECREF(hits_array); + + return NULL; + + } + + /*If the call succeeded, build the output tuple and quit*/ + PyObject *output = PyTuple_New(2); + if(PyTuple_SetItem(output,0,hits_array)){ + + Py_DECREF(ft_map1_array); + Py_DECREF(ft_map2_array); + Py_DECREF(ft_map3_array); + Py_DECREF(ksmall_array); + Py_DECREF(klarge_array); + Py_XDECREF(bispectrum_array); + Py_XDECREF(hits_array); + + return NULL; + + } + + if(PyTuple_SetItem(output,1,bispectrum_array)){ + + Py_DECREF(ft_map1_array); + Py_DECREF(ft_map2_array); + Py_DECREF(ft_map3_array); + Py_DECREF(ksmall_array); + Py_DECREF(klarge_array); + Py_XDECREF(bispectrum_array); + Py_XDECREF(hits_array); + + return NULL; + + } + + //Cleanup and return + Py_DECREF(ft_map1_array); + Py_DECREF(ft_map2_array); + Py_DECREF(ft_map3_array); + Py_DECREF(ksmall_array); + Py_DECREF(klarge_array); + + return output; +} + + +// rfft3_bispectrum() implementation +static PyObject *_topology_rfft3_bispectrum(PyObject *self,PyObject *args){ + + /*These are the inputs: the Fourier transforms of the two maps, the size of the pixel in k space, the k bin extremes at which calculate the azimuthal averages*/ + PyObject *ft_map1_obj,*ft_map2_obj,*ft_map3_obj, *ksmall_obj, *klarge_obj, *cosine_obj; + double kpixX,kpixY,kpixZ; + int size_small_x, size_small_y, size_small_z; + + /*Parse input tuple*/ + if(!PyArg_ParseTuple(args,"OOOdddiiiOOO",&ft_map1_obj,&ft_map2_obj,&ft_map3_obj, &kpixX,&kpixY,&kpixZ, &size_small_x, &size_small_y, &size_small_z,&ksmall_obj, &klarge_obj, &cosine_obj)){ + return NULL; + } + + /*Interpret the parsed objects as numpy arrays*/ + PyObject *ft_map1_array = PyArray_FROM_OTF(ft_map1_obj,NPY_COMPLEX128,NPY_IN_ARRAY); + PyObject *ft_map2_array = PyArray_FROM_OTF(ft_map2_obj,NPY_COMPLEX128,NPY_IN_ARRAY); + PyObject *ft_map3_array = PyArray_FROM_OTF(ft_map3_obj,NPY_COMPLEX128,NPY_IN_ARRAY); + PyObject *ksmall_array = PyArray_FROM_OTF(ksmall_obj,NPY_DOUBLE,NPY_IN_ARRAY); + PyObject *klarge_array = PyArray_FROM_OTF(klarge_obj,NPY_DOUBLE,NPY_IN_ARRAY); + PyObject *cosine_array = PyArray_FROM_OTF(cosine_obj,NPY_DOUBLE,NPY_IN_ARRAY); + + /*Check if anything failed*/ + if(ft_map1_array==NULL || ksmall_array==NULL || ft_map2_array==NULL || ft_map3_array==NULL || klarge_array==NULL || cosine_array==NULL){ + + Py_XDECREF(ft_map1_array); + Py_XDECREF(ft_map2_array); + Py_XDECREF(ft_map3_array); + Py_XDECREF(ksmall_array); + Py_XDECREF(klarge_array); + Py_XDECREF(cosine_array); + + return NULL; + } + + /*Get the size of the map fourier transform*/ + int Nside_x = (int)PyArray_DIM(ft_map1_array,0); + int Nside_y = (int)PyArray_DIM(ft_map1_array,1); + int Nside_z = (int)PyArray_DIM(ft_map1_array,2); + + /*Get the number of k bin edges*/ + int Nlarge = (int)PyArray_DIM(klarge_array,0); + int Nsmall = (int)PyArray_DIM(ksmall_array,0); + int Ncosine = (int)PyArray_DIM(cosine_array,0); + + /*Build the arrays that will contain the output:bispectrum and hits*/ + npy_intp dims[] = {(npy_intp) Nsmall - 1,(npy_intp) Nlarge-1, (npy_intp) Ncosine-1}; + PyObject *bispectrum_array = PyArray_ZEROS(3,dims,NPY_DOUBLE,0); + PyObject *hits_array = PyArray_ZEROS(3,dims,NPY_LONG,0); + + /*Check for failure*/ + if(bispectrum_array==NULL || hits_array==NULL){ + + Py_DECREF(ft_map1_array); + Py_DECREF(ft_map2_array); + Py_DECREF(ft_map3_array); + Py_DECREF(ksmall_array); + Py_DECREF(klarge_array); + Py_DECREF(cosine_array); + Py_XDECREF(bispectrum_array); + Py_XDECREF(hits_array); + //Py_XDECREF(hits_test); + + return NULL; + } + + /*Call the C backend bispectrum function*/ + if(bispectrum_rfft3((double _Complex *)PyArray_DATA(ft_map1_array),(double _Complex *)PyArray_DATA(ft_map2_array),(double _Complex *)PyArray_DATA(ft_map3_array),Nside_x,Nside_y,Nside_z,size_small_x,size_small_y,size_small_z,kpixX,kpixY,kpixZ,Nsmall,Nlarge,Ncosine,(double *)PyArray_DATA(ksmall_array),(double *)PyArray_DATA(klarge_array),(double *)PyArray_DATA(cosine_array),(double *)PyArray_DATA(bispectrum_array),(long *)PyArray_DATA(hits_array))){ + + Py_DECREF(ft_map1_array); + Py_DECREF(ft_map2_array); + Py_DECREF(ft_map3_array); + Py_DECREF(ksmall_array); + Py_DECREF(klarge_array); + Py_DECREF(cosine_array); + Py_XDECREF(bispectrum_array); + Py_XDECREF(hits_array); + + return NULL; + + } + + /*If the call succeeded, build the output tuple and quit*/ + PyObject *output = PyTuple_New(2); + if(PyTuple_SetItem(output,0,hits_array)){ + + Py_DECREF(ft_map1_array); + Py_DECREF(ft_map2_array); + Py_DECREF(ft_map3_array); + Py_DECREF(ksmall_array); + Py_DECREF(klarge_array); + Py_DECREF(cosine_array); + Py_XDECREF(bispectrum_array); + Py_XDECREF(hits_array); + + return NULL; + + } + + if(PyTuple_SetItem(output,1,bispectrum_array)){ + + Py_DECREF(ft_map1_array); + Py_DECREF(ft_map2_array); + Py_DECREF(ft_map3_array); + Py_DECREF(ksmall_array); + Py_DECREF(klarge_array); + Py_DECREF(cosine_array); + Py_XDECREF(bispectrum_array); + Py_XDECREF(hits_array); + + return NULL; + + } + + //Cleanup and return + Py_DECREF(ft_map1_array); + Py_DECREF(ft_map2_array); + Py_DECREF(ft_map3_array); + Py_DECREF(ksmall_array); + Py_DECREF(klarge_array); + Py_DECREF(cosine_array); + + return output; +} + + + diff --git a/lenstools/extern/azimuth.c b/lenstools/extern/azimuth.c index 9982acd9..bb96f1a4 100644 --- a/lenstools/extern/azimuth.c +++ b/lenstools/extern/azimuth.c @@ -116,4 +116,820 @@ int azimuthal_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,long size_ return 0; +} + + +/* Compute angular averaged bispectrum of 3D real Fourier transforms of scalar fields (i.e. density)*/ +int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _Complex *ft_map3,long size_x,long size_y,long size_z,int size_small_x, int size_small_y, int size_small_z, double kpixX,double kpixY,double kpixZ,int Nsmall,int Nlarge, double *ksmall_edges, double *klarge_edges ,double *bispectrum_k,long *hits){ + // ft_map1, ft_map2, ft_map3 : input spectra + // size_x, size_y, size_z : the grid size of ft_map + // size_small_x, size_small_y, size_small_z : the maximum grid points that belong to ksmall_edges array; needed to reduced calculation time + // kpixX, kpixY, kpixZ = 2 Pi / box_size : minimum momentum + // Nsmall, Nlarge : the length of ksmall_edges or klarge_edges correspondingly + // ksmall_edges, klarge_edges : small or large momentum at which bispectrum are computed + // bispectrum_k, hits : results + + + //Counters + long x1,y1,z1,x2,y2,z2,x3,y3,z3,p1,p2,p3; + int b_large,b_small; + double k1,k1x,k1y,k1z,k2,k2x,k2y,k2z,k3,k3x,k3y,k3z; + int odd_or_even; + + + //Binning + int largebins = Nlarge - 1; + int smallbins = Nsmall - 1; + + if(size_x % 2 == 0){ + odd_or_even = 1; //even + } + else{ + odd_or_even = 0; //odd + } + + //Loop over all the pixels in the fourier map + for(x1=0;x1ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){ + for(x2=0;x2=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + bispectrum_k[b_small*largebins + b_large] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins + b_large]++; + } + } + } + } + else{ + if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + bispectrum_k[b_small*largebins + b_large] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins + b_large]++; + } + } + } + } + } + } + } + } + } + } + } + + for(y1=size_y-size_small_y+1;y1ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){ + for(x2=0;x2=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + bispectrum_k[b_small*largebins + b_large] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins + b_large]++; + } + } + } + } + else{ + if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + bispectrum_k[b_small*largebins + b_large] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins + b_large]++; + } + } + } + } + } + } + } + } + } + } + } + } + + for(x1=size_x - size_small_x + 1;x1ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){ + for(x2=0;x2=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + bispectrum_k[b_small*largebins + b_large] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins + b_large]++; + } + } + } + } + else{ + if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + bispectrum_k[b_small*largebins + b_large] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins + b_large]++; + } + } + } + } + } + } + } + } + } + } + } + + for(y1=size_y-size_small_y+1;y1ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){ + for(x2=0;x2=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + bispectrum_k[b_small*largebins + b_large] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins + b_large]++; + } + } + } + } + else{ + if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + bispectrum_k[b_small*largebins + b_large] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins + b_large]++; + } + } + } + } + } + } + } + } + } + } + } + } + //Return + return 0; +} + +/* Compute bispectrum of 3D real Fourier transforms of scalar fields (i.e. density)*/ +int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _Complex *ft_map3,long size_x,long size_y, long size_z, int size_small_x, int size_small_y, int size_small_z, double kpixX,double kpixY,double kpixZ,int Nsmall, int Nlarge, int Ncosine, double *ksmall_edges, double *klarge_edges, double *cosines, double *bispectrum_k,long *hits){ + // ft_map1, ft_map2, ft_map3 : input spectra + // size_x, size_y, size_z : the grid size of ft_map + // size_small_x, size_small_y, size_small_z : the maximum grid points that belong to ksmall_edges array; needed to reduced calculation time + // kpixX, kpixY, kpixZ = 2 Pi / box_size : minimum momentum + // Nsmall, Nlarge, Ncosine : the length of ksmall_edges, klarge_edges or cosines correspondingly + // ksmall_edges, klarge_edges : small or large momentum at which bispectrum are computed + // cosines : cosine of angle between momenta at which bispectrum are computed + // bispectrum_k, hits : results + + + //Counters + long x1,y1,z1,x2,y2,z2,x3,y3,z3,p1,p2,p3; + int b_large,b_small; + long b_cosine; + double k1,k1x,k1y,k1z,k2,k2x,k2y,k2z,k3,k3x,k3y,k3z,cos_theta; + int odd_or_even; + + //Binning + int largebins = Nlarge - 1; + int smallbins = Nsmall - 1; + int cosinebins = Ncosine - 1; + + if(size_x % 2 == 0){ + odd_or_even = 1; //even + } + else{ + odd_or_even = 0; //odd + } + //Loop over all the pixels in the fourier map + for(x1=0;x1ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){ + for(x2=0;x2=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2); + for(b_cosine=0;b_cosinecosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ + bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; + //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); + //hits_test[b_cosine]++; + } + } + } + } + } + } + else{ + if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2); + for(b_cosine=0;b_cosinecosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ + bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; + //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); + //hits_test[b_cosine]++; + } + } + } + } + } + } + } + } + } + } + } + } + } + + for(y1=size_y - size_small_y+1;y1ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){ + for(x2=0;x2=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2); + for(b_cosine=0;b_cosinecosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ + bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; + //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); + //hits_test[b_cosine]++; + } + } + } + } + } + } + else{ + if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2); + for(b_cosine=0;b_cosinecosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ + bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; + //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); + //hits_test[b_cosine]++; + } + } + } + } + } + } + } + } + } + } + } + } + } + } + + for(x1=size_x - size_small_x+1;x1ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){ + for(x2=0;x2=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2); + for(b_cosine=0;b_cosinecosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ + bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; + //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); + //hits_test[b_cosine]++; + } + } + } + } + } + } + else{ + if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2); + for(b_cosine=0;b_cosinecosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ + bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; + //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); + //hits_test[b_cosine]++; + } + } + } + } + } + } + } + } + } + } + } + } + } + + for(y1=size_y - size_small_y+1;y1ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){ + for(x2=0;x2=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2); + for(b_cosine=0;b_cosinecosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ + bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; + //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); + //hits_test[b_cosine]++; + } + } + } + } + } + } + else{ + if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=0){ + x3 = (long)(k3x/kpixX); + } + else{ + x3 = size_x + (long)(k3x/kpixX); + } + if(k3y>=0){ + y3 = (long)(k3y/kpixY); + } + else{ + y3 = size_y + (long)(k3y/kpixY); + } + z3 = (long)(k3z/kpixZ); + p3 = x3*size_y*size_z + y3*size_z + z3; + for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ + cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2); + for(b_cosine=0;b_cosinecosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ + bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); + hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; + //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); + //hits_test[b_cosine]++; + } + } + } + } + } + } + } + } + } + } + } + } + } + } + //Return + return 0; } \ No newline at end of file diff --git a/lenstools/extern/azimuth.h b/lenstools/extern/azimuth.h index 7062bf92..872628fa 100644 --- a/lenstools/extern/azimuth.h +++ b/lenstools/extern/azimuth.h @@ -5,5 +5,7 @@ int azimuthal_rfft2(double _Complex *ft_map1,double _Complex *ft_map2,long size_x,long size_y,double map_angle_degrees,int Nvalues,double *lvalues,double *power_l); int azimuthal_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,long size_x,long size_y,long size_z,double kpixX,double kpixY,double kpixZ,int Nvalues,double *kvalues,double *power_k,long *hits); +int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _Complex *ft_map3,long size_x,long size_y,long size_z,int size_small_x, int size_small_y, int size_small_z, double kpixX,double kpixY,double kpixZ,int Nsmall,int Nlarge, double *ksmall_edges, double *klarge_edges ,double *bispectrum_k,long *hits); +int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _Complex *ft_map3,long size_x,long size_y,long size_z,int size_small_x, int size_small_y, int size_small_z,double kpixX,double kpixY,double kpixZ,int Nsmall, int Nlarge, int Ncosine, double *ksmall_edges, double *klarge_edges, double *cosines, double *bispectrum_k,long *hits); #endif \ No newline at end of file From a2444f7e0d4d5b2e122dbb77e6804eb327d1ddee Mon Sep 17 00:00:00 2001 From: jongheekang Date: Mon, 4 May 2015 15:45:15 -0400 Subject: [PATCH 4/6] Fix the indices issue of bispectrum function --- lenstools/extern/azimuth.c | 180 +++++++++++++++++-------------------- 1 file changed, 80 insertions(+), 100 deletions(-) diff --git a/lenstools/extern/azimuth.c b/lenstools/extern/azimuth.c index bb96f1a4..8f662708 100644 --- a/lenstools/extern/azimuth.c +++ b/lenstools/extern/azimuth.c @@ -182,18 +182,18 @@ int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2 if(odd_or_even == 1){ if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ @@ -206,18 +206,18 @@ int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2 else{ if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ @@ -267,18 +267,18 @@ int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2 if(odd_or_even == 1){ if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ @@ -291,18 +291,18 @@ int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2 else{ if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ @@ -354,18 +354,18 @@ int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2 if(odd_or_even == 1){ if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ @@ -378,18 +378,18 @@ int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2 else{ if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ @@ -439,18 +439,18 @@ int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2 if(odd_or_even == 1){ if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ @@ -463,18 +463,18 @@ int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2 else{ if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ @@ -560,20 +560,19 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C if(odd_or_even == 1){ if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; - for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2); @@ -581,8 +580,6 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C if(cos_theta>cosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; - //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); - //hits_test[b_cosine]++; } } } @@ -592,18 +589,18 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C else{ if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ @@ -612,8 +609,6 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C if(cos_theta>cosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; - //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); - //hits_test[b_cosine]++; } } } @@ -660,20 +655,19 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C if(odd_or_even == 1){ if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; - for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2); @@ -681,8 +675,6 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C if(cos_theta>cosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; - //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); - //hits_test[b_cosine]++; } } } @@ -692,18 +684,18 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C else{ if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ @@ -712,8 +704,6 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C if(cos_theta>cosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; - //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); - //hits_test[b_cosine]++; } } } @@ -762,20 +752,19 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C if(odd_or_even == 1){ if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; - for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2); @@ -783,8 +772,6 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C if(cos_theta>cosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; - //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); - //hits_test[b_cosine]++; } } } @@ -794,18 +781,18 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C else{ if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ @@ -814,8 +801,6 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C if(cos_theta>cosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; - //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); - //hits_test[b_cosine]++; } } } @@ -862,20 +847,19 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C if(odd_or_even == 1){ if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; - for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2); @@ -883,8 +867,6 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C if(cos_theta>cosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; - //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); - //hits_test[b_cosine]++; } } } @@ -894,18 +876,18 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C else{ if(k3x>=-kpixX*(size_x-1)/2 && k3x <= kpixX*(size_x-1)/2 && k3y>=-kpixY*(size_y-1)/2 && k3y<=kpixY*(size_y-1)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ - x3 = (long)(k3x/kpixX); + x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } else{ - x3 = size_x + (long)(k3x/kpixX); + x3 = size_x + (long)(k3x/kpixX - 0.001); } if(k3y>=0){ - y3 = (long)(k3y/kpixY); + y3 = (long)(k3y/kpixY + 0.001); } else{ - y3 = size_y + (long)(k3y/kpixY); + y3 = size_y + (long)(k3y/kpixY - 0.001); } - z3 = (long)(k3z/kpixZ); + z3 = (long)(k3z/kpixZ + 0.001); p3 = x3*size_y*size_z + y3*size_z + z3; for(b_large=0;b_largeklarge_edges[b_large] && k2<=klarge_edges[b_large+1]){ @@ -914,8 +896,6 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C if(cos_theta>cosines[b_cosine] && cos_theta<=cosines[b_cosine+1]){ bispectrum_k[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine] += creal(ft_map1[p1]*ft_map2[p2]*conj(ft_map3[p3])); hits[b_small*largebins*cosinebins + b_large * cosinebins + b_cosine]++; - //printf("x3 = %ld, y3 = %ld, z3 = %ld, b_small = %d, b_large = %d, b_cosine = %ld, cos_theta = %lf\n",x3,y3,z3,b_small,b_large,b_cosine,cos_theta); - //hits_test[b_cosine]++; } } } From c4812663d8c52c8cdd1937a371b4a6db9f3e3470 Mon Sep 17 00:00:00 2001 From: jongheekang Date: Wed, 6 May 2015 13:21:54 -0400 Subject: [PATCH 5/6] Modify signed_long and int functions --- lenstools/extern/coordinates.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lenstools/extern/coordinates.c b/lenstools/extern/coordinates.c index cb2b75ab..dc8131a7 100644 --- a/lenstools/extern/coordinates.c +++ b/lenstools/extern/coordinates.c @@ -27,7 +27,7 @@ long min_long(long i,long j){ int signed_int(int i,int j){ - if(i<=j){ + if(i Date: Wed, 6 May 2015 13:26:24 -0400 Subject: [PATCH 6/6] Modify azimuth.c(Modify frequency value conditions in if statement when the size is even) --- lenstools/extern/azimuth.c | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/lenstools/extern/azimuth.c b/lenstools/extern/azimuth.c index 8f662708..be4ed061 100644 --- a/lenstools/extern/azimuth.c +++ b/lenstools/extern/azimuth.c @@ -180,7 +180,7 @@ int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2 k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z); if(odd_or_even == 1){ - if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=-kpixX*size_x/2 && k3x <= kpixX*(size_x-2)/2 && k3y>=-kpixY*size_y/2 && k3y<=kpixY*(size_y-2)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } @@ -265,7 +265,7 @@ int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2 k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z); if(odd_or_even == 1){ - if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=-kpixX*size_x/2 && k3x <= kpixX*(size_x-2)/2 && k3y>=-kpixY*size_y/2 && k3y<=kpixY*(size_y-2)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } @@ -352,7 +352,7 @@ int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2 k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z); if(odd_or_even == 1){ - if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=-kpixX*size_x/2 && k3x <= kpixX*(size_x-2)/2 && k3y>=-kpixY*size_y/2 && k3y<=kpixY*(size_y-2)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } @@ -437,7 +437,7 @@ int azimuthal_bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2 k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z); if(odd_or_even == 1){ - if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=-kpixX*size_x/2 && k3x <= kpixX*(size_x-2)/2 && k3y>=-kpixY*size_y/2 && k3y<=kpixY*(size_y-2)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } @@ -558,7 +558,7 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z); if(odd_or_even == 1){ - if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=-kpixX*size_x/2 && k3x <= kpixX*(size_x-2)/2 && k3y>=-kpixY*size_y/2 && k3y<=kpixY*(size_y-2)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } @@ -653,7 +653,7 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z); if(odd_or_even == 1){ - if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=-kpixX*size_x/2 && k3x <= kpixX*(size_x-2)/2 && k3y>=-kpixY*size_y/2 && k3y<=kpixY*(size_y-2)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } @@ -750,7 +750,7 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z); if(odd_or_even == 1){ - if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=-kpixX*size_x/2 && k3x <= kpixX*(size_x-2)/2 && k3y>=-kpixY*size_y/2 && k3y<=kpixY*(size_y-2)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers } @@ -845,7 +845,7 @@ int bispectrum_rfft3(double _Complex *ft_map1,double _Complex *ft_map2,double _C k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z); if(odd_or_even == 1){ - if(k3x>=-kpixX*(size_x-2)/2 && k3x <= kpixX*size_x/2 && k3y>=-kpixY*(size_y-2)/2 && k3y<=kpixY*size_y/2 && k3z<=kpixZ*(size_z-1)){ + if(k3x>=-kpixX*size_x/2 && k3x <= kpixX*(size_x-2)/2 && k3y>=-kpixY*size_y/2 && k3y<=kpixY*(size_y-2)/2 && k3z<=kpixZ*(size_z-1)){ if(k3x>=0){ x3 = (long)(k3x/kpixX + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers }