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/_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..be4ed061 100644
--- a/lenstools/extern/azimuth.c
+++ b/lenstools/extern/azimuth.c
@@ -116,4 +116,800 @@ 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;x1<size_small_x;x1++){
+		for(y1=0;y1<size_small_y;y1++){
+			for(z1=0;z1<size_small_z;z1++){
+
+				// k1 is small momentum
+				k1x = signed_long(x1,size_x-x1) * kpixX;
+				k1y = signed_long(y1,size_y-y1) * kpixY;
+				k1z = z1 * kpixZ;
+				k1 = sqrt(k1x*k1x + k1y*k1y + k1z*k1z);
+				p1 = x1*size_y*size_z + y1*size_z + z1;
+
+				//Decide in which bin_small this small momentum pixel falls into
+				for(b_small=0;b_small<smallbins;b_small++){
+					if(k1>ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){
+						for(x2=0;x2<size_x;x2++){
+							for(y2=0;y2<size_y;y2++){
+								for(z2=0;z2<size_z;z2++){
+
+									// k2 is large momentum
+									k2x = signed_long(x2,size_x-x2) * kpixX;
+									k2y = signed_long(y2,size_y-y2) * kpixY;
+									k2z = z2 * kpixZ;
+									k2 = sqrt(k2x*k2x + k2y*k2y + k2z*k2z);
+									p2 = x2*size_y*size_z + y2*size_z + z2;
+									// k3 = k1 + k2
+									k3x = k1x+k2x;
+									k3y = k1y+k2y;
+									k3z = k1z+k2z;
+									k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z);
+
+									if(odd_or_even == 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
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_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 + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_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;y1<size_y;y1++){
+			for(z1=0;z1<size_small_z;z1++){
+
+				// k1 is small momentum
+				k1x = signed_long(x1,size_x-x1) * kpixX;
+				k1y = signed_long(y1,size_y-y1) * kpixY;
+				k1z = z1 * kpixZ;
+				k1 = sqrt(k1x*k1x + k1y*k1y + k1z*k1z);
+				p1 = x1*size_y*size_z + y1*size_z + z1;
+
+				//Decide in which bin_small this small momentum pixel falls into
+				for(b_small=0;b_small<smallbins;b_small++){
+					if(k1>ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){
+						for(x2=0;x2<size_x;x2++){
+							for(y2=0;y2<size_y;y2++){
+								for(z2=0;z2<size_z;z2++){
+
+									// k2 is large momentum
+									k2x = signed_long(x2,size_x-x2) * kpixX;
+									k2y = signed_long(y2,size_y-y2) * kpixY;
+									k2z = z2 * kpixZ;
+									k2 = sqrt(k2x*k2x + k2y*k2y + k2z*k2z);
+									p2 = x2*size_y*size_z + y2*size_z + z2;
+									// k3 = k1 + k2
+									k3x = k1x+k2x;
+									k3y = k1y+k2y;
+									k3z = k1z+k2z;
+									k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z);
+
+									if(odd_or_even == 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
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_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 + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_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;x1<size_x;x1++){
+		for(y1=0;y1<size_small_y;y1++){
+			for(z1=0;z1<size_small_z;z1++){
+
+				// k1 is small momentum
+				k1x = signed_long(x1,size_x-x1) * kpixX;
+				k1y = signed_long(y1,size_y-y1) * kpixY;
+				k1z = z1 * kpixZ;
+				k1 = sqrt(k1x*k1x + k1y*k1y + k1z*k1z);
+				p1 = x1*size_y*size_z + y1*size_z + z1;
+
+				//Decide in which bin_small this small momentum pixel falls into
+				for(b_small=0;b_small<smallbins;b_small++){
+					if(k1>ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){
+						for(x2=0;x2<size_x;x2++){
+							for(y2=0;y2<size_y;y2++){
+								for(z2=0;z2<size_z;z2++){
+
+									// k2 is large momentum
+									k2x = signed_long(x2,size_x-x2) * kpixX;
+									k2y = signed_long(y2,size_y-y2) * kpixY;
+									k2z = z2 * kpixZ;
+									k2 = sqrt(k2x*k2x + k2y*k2y + k2z*k2z);
+									p2 = x2*size_y*size_z + y2*size_z + z2;
+									// k3 = k1 + k2
+									k3x = k1x+k2x;
+									k3y = k1y+k2y;
+									k3z = k1z+k2z;
+									k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z);
+
+									if(odd_or_even == 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
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_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 + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_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;y1<size_y;y1++){
+			for(z1=0;z1<size_small_z;z1++){
+
+				// k1 is small momentum
+				k1x = signed_long(x1,size_x-x1) * kpixX;
+				k1y = signed_long(y1,size_y-y1) * kpixY;
+				k1z = z1 * kpixZ;
+				k1 = sqrt(k1x*k1x + k1y*k1y + k1z*k1z);
+				p1 = x1*size_y*size_z + y1*size_z + z1;
+
+				//Decide in which bin_small this small momentum pixel falls into
+				for(b_small=0;b_small<smallbins;b_small++){
+					if(k1>ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){
+						for(x2=0;x2<size_x;x2++){
+							for(y2=0;y2<size_y;y2++){
+								for(z2=0;z2<size_z;z2++){
+
+									// k2 is large momentum
+									k2x = signed_long(x2,size_x-x2) * kpixX;
+									k2y = signed_long(y2,size_y-y2) * kpixY;
+									k2z = z2 * kpixZ;
+									k2 = sqrt(k2x*k2x + k2y*k2y + k2z*k2z);
+									p2 = x2*size_y*size_z + y2*size_z + z2;
+									// k3 = k1 + k2
+									k3x = k1x+k2x;
+									k3y = k1y+k2y;
+									k3z = k1z+k2z;
+									k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z);
+
+									if(odd_or_even == 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
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_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 + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_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;x1<size_small_x;x1++){
+		for(y1=0;y1<size_small_y;y1++){
+			for(z1=0;z1<size_small_z;z1++){
+
+				// k1 is small momentum
+				k1x = signed_long(x1,size_x-x1) * kpixX;
+				k1y = signed_long(y1,size_y-y1) * kpixY;
+				k1z = z1 * kpixZ;
+				k1 = sqrt(k1x*k1x + k1y*k1y + k1z*k1z);
+				p1 = x1*size_y*size_z + y1*size_z + z1;
+
+				//Decide in which bin_small this small momentum pixel falls into
+				for(b_small=0;b_small<smallbins;b_small++){
+					if(k1>ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){
+						for(x2=0;x2<size_x;x2++){
+							for(y2=0;y2<size_y;y2++){
+								for(z2=0;z2<size_z;z2++){
+
+									// k2 is large momentum
+									k2x = signed_long(x2,size_x-x2) * kpixX;
+									k2y = signed_long(y2,size_y-y2) * kpixY;
+									k2z = z2 * kpixZ;
+									k2 = sqrt(k2x*k2x + k2y*k2y + k2z*k2z);
+									p2 = x2*size_y*size_z + y2*size_z + z2;
+									// k3 = k1 + k2
+									k3x = k1x+k2x;
+									k3y = k1y+k2y;
+									k3z = k1z+k2z;
+									k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z);
+
+									if(odd_or_even == 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
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_edges[b_large] && k2<=klarge_edges[b_large+1]){
+													cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2);
+													for(b_cosine=0;b_cosine<cosinebins;b_cosine++){
+														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]++;
+														}
+													}
+												}
+											}
+										}
+									}
+									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 + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_edges[b_large] && k2<=klarge_edges[b_large+1]){
+													cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2);
+													for(b_cosine=0;b_cosine<cosinebins;b_cosine++){
+														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]++;
+														}
+													}	
+												}
+											}							
+										}
+									}
+								}
+							}
+						}
+					}
+				}
+			}
+		}
+
+		for(y1=size_y - size_small_y+1;y1<size_y;y1++){
+			for(z1=0;z1<size_small_z;z1++){
+
+				// k1 is small momentum
+				k1x = signed_long(x1,size_x-x1) * kpixX;
+				k1y = signed_long(y1,size_y-y1) * kpixY;
+				k1z = z1 * kpixZ;
+				k1 = sqrt(k1x*k1x + k1y*k1y + k1z*k1z);
+				p1 = x1*size_y*size_z + y1*size_z + z1;
+
+				//Decide in which bin_small this small momentum pixel falls into
+				for(b_small=0;b_small<smallbins;b_small++){
+					if(k1>ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){
+						for(x2=0;x2<size_x;x2++){
+							for(y2=0;y2<size_y;y2++){
+								for(z2=0;z2<size_z;z2++){
+
+									// k2 is large momentum
+									k2x = signed_long(x2,size_x-x2) * kpixX;
+									k2y = signed_long(y2,size_y-y2) * kpixY;
+									k2z = z2 * kpixZ;
+									k2 = sqrt(k2x*k2x + k2y*k2y + k2z*k2z);
+									p2 = x2*size_y*size_z + y2*size_z + z2;
+									// k3 = k1 + k2
+									k3x = k1x+k2x;
+									k3y = k1y+k2y;
+									k3z = k1z+k2z;
+									k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z);
+
+									if(odd_or_even == 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
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_edges[b_large] && k2<=klarge_edges[b_large+1]){
+													cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2);
+													for(b_cosine=0;b_cosine<cosinebins;b_cosine++){
+														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]++;
+														}
+													}
+												}
+											}
+										}
+									}
+									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 + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_edges[b_large] && k2<=klarge_edges[b_large+1]){
+													cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2);
+													for(b_cosine=0;b_cosine<cosinebins;b_cosine++){
+														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]++;
+														}
+													}	
+												}
+											}							
+										}
+									}
+								}
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+
+	for(x1=size_x - size_small_x+1;x1<size_x;x1++){
+		for(y1=0;y1<size_small_y;y1++){
+			for(z1=0;z1<size_small_z;z1++){
+
+				// k1 is small momentum
+				k1x = signed_long(x1,size_x-x1) * kpixX;
+				k1y = signed_long(y1,size_y-y1) * kpixY;
+				k1z = z1 * kpixZ;
+				k1 = sqrt(k1x*k1x + k1y*k1y + k1z*k1z);
+				p1 = x1*size_y*size_z + y1*size_z + z1;
+
+				//Decide in which bin_small this small momentum pixel falls into
+				for(b_small=0;b_small<smallbins;b_small++){
+					if(k1>ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){
+						for(x2=0;x2<size_x;x2++){
+							for(y2=0;y2<size_y;y2++){
+								for(z2=0;z2<size_z;z2++){
+
+									// k2 is large momentum
+									k2x = signed_long(x2,size_x-x2) * kpixX;
+									k2y = signed_long(y2,size_y-y2) * kpixY;
+									k2z = z2 * kpixZ;
+									k2 = sqrt(k2x*k2x + k2y*k2y + k2z*k2z);
+									p2 = x2*size_y*size_z + y2*size_z + z2;
+									// k3 = k1 + k2
+									k3x = k1x+k2x;
+									k3y = k1y+k2y;
+									k3z = k1z+k2z;
+									k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z);
+
+									if(odd_or_even == 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
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_edges[b_large] && k2<=klarge_edges[b_large+1]){
+													cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2);
+													for(b_cosine=0;b_cosine<cosinebins;b_cosine++){
+														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]++;
+														}
+													}
+												}
+											}
+										}
+									}
+									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 + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_edges[b_large] && k2<=klarge_edges[b_large+1]){
+													cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2);
+													for(b_cosine=0;b_cosine<cosinebins;b_cosine++){
+														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]++;
+														}
+													}	
+												}
+											}							
+										}
+									}
+								}
+							}
+						}
+					}
+				}
+			}
+		}
+
+		for(y1=size_y - size_small_y+1;y1<size_y;y1++){
+			for(z1=0;z1<size_small_z;z1++){
+
+				// k1 is small momentum
+				k1x = signed_long(x1,size_x-x1) * kpixX;
+				k1y = signed_long(y1,size_y-y1) * kpixY;
+				k1z = z1 * kpixZ;
+				k1 = sqrt(k1x*k1x + k1y*k1y + k1z*k1z);
+				p1 = x1*size_y*size_z + y1*size_z + z1;
+
+				//Decide in which bin_small this small momentum pixel falls into
+				for(b_small=0;b_small<smallbins;b_small++){
+					if(k1>ksmall_edges[b_small] && k1<=ksmall_edges[b_small+1]){
+						for(x2=0;x2<size_x;x2++){
+							for(y2=0;y2<size_y;y2++){
+								for(z2=0;z2<size_z;z2++){
+
+									// k2 is large momentum
+									k2x = signed_long(x2,size_x-x2) * kpixX;
+									k2y = signed_long(y2,size_y-y2) * kpixY;
+									k2z = z2 * kpixZ;
+									k2 = sqrt(k2x*k2x + k2y*k2y + k2z*k2z);
+									p2 = x2*size_y*size_z + y2*size_z + z2;
+									// k3 = k1 + k2
+									k3x = k1x+k2x;
+									k3y = k1y+k2y;
+									k3z = k1z+k2z;
+									k3 = sqrt(k3x*k3x + k3y*k3y + k3z*k3z);
+
+									if(odd_or_even == 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
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_edges[b_large] && k2<=klarge_edges[b_large+1]){
+													cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2);
+													for(b_cosine=0;b_cosine<cosinebins;b_cosine++){
+														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]++;
+														}
+													}
+												}
+											}
+										}
+									}
+									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 + 0.001); // Add or subtract 0.001 to avoid error due to finite precision of real numbers
+											}
+											else{
+												x3 = size_x + (long)(k3x/kpixX - 0.001);
+											}
+											if(k3y>=0){
+												y3 = (long)(k3y/kpixY + 0.001);
+											}
+											else{
+												y3 = size_y + (long)(k3y/kpixY - 0.001);
+											}
+											z3 = (long)(k3z/kpixZ + 0.001);
+											p3 = x3*size_y*size_z + y3*size_z + z3;
+											for(b_large=0;b_large<largebins;b_large++){
+												if(k2>klarge_edges[b_large] && k2<=klarge_edges[b_large+1]){
+													cos_theta = (k1x*k2x + k1y*k2y + k1z*k2z)/(k1*k2);
+													for(b_cosine=0;b_cosine<cosinebins;b_cosine++){
+														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]++;
+														}
+													}	
+												}
+											}							
+										}
+									}
+								}
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+	//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
diff --git a/lenstools/extern/coordinates.c b/lenstools/extern/coordinates.c
index 63f4b871..dc8131a7 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);
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<Npart;n++){
+
+		//Compute the position on the grid in the fastest way
+		i = (positions[3*n] - leftX)/sizeX;
+		j = (positions[3*n + 1] - leftY)/sizeY;
+		k = (positions[3*n + 2] - leftZ)/sizeZ;
+
+		//If the particle lands on the grid, put it in the correct pixel
+		if(i>=0 && i<nx && j>=0 && j<ny && k>=0 && k<nz){
+
+			grid[((int)i)*ny*nz + ((int)j)*nz + (int)k] += 1.0 / volume_grid;
+
+		}
+
+
+	}
+
+
+	return 0;
+
+
+}
+
+//Snap particles on a 3d regularly spaced grid and compute the mass density fluctuation within each 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 n,a,b,c;
+	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<Npart;n++){
+
+		//Compute the position on the grid in the fastest way
+		i = (positions[3*n] - leftX)/sizeX;
+		j = (positions[3*n + 1] - leftY)/sizeY;
+		k = (positions[3*n + 2] - leftZ)/sizeZ;
+
+		//If the particle lands on the grid, put it in the correct pixel
+		if(i>=0 && i<nx && j>=0 && j<ny && k>=0 && k<nz){
+
+			grid[((int)i)*ny*nz + ((int)j)*nz + (int)k] += 1.0 / volume_grid;
+
+		}
+
+
+	}
+
+	for(a=0;a<nx;a++){
+		for(b=0;b<ny;b++){
+			for(c=0;c<nz;c++){
+				grid[a*ny*nz + b*nz + c] = grid[a*ny*nz + b*nz + c]/rho_0 - 1;
+			}
+		}
+	}
+
+
+	return 0;
+
+
+}
 
 //adaptive smoothing
 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){
diff --git a/lenstools/extern/grid.h b/lenstools/extern/grid.h
index ae8ab50c..7dfaec87 100644
--- a/lenstools/extern/grid.h
+++ b/lenstools/extern/grid.h
@@ -4,6 +4,8 @@
 #include <math.h>
 
 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);