diff --git a/src/common.c b/src/common.c
index b7ea7b2..7cd8f16 100644
--- a/src/common.c
+++ b/src/common.c
@@ -21,6 +21,10 @@ int TaskSmooth=0;
 float *Dens_sm_local;
 fftwf_complex *Cdens_sm_local;
 
+int TaskFractalD=0;
+float *Dens_smD_local;
+fftwf_complex *Cdens_smD_local;
+
 int TaskTidal=0;
 int TaskTidalDiag=0;
 float **Tid_local;
@@ -57,6 +61,8 @@ void end_all(void)
   fftwf_free(Dens_local);
   if(TaskSmooth)
     fftwf_free(Dens_sm_local);
+  if(TaskFractalD)
+    fftwf_free(Dens_smD_local);
   if(TaskTidal) {
     for(i=0;i<6;i++)
       fftwf_free(Tid_local[i]);
@@ -113,6 +119,12 @@ void domain_decomp(void)
       report_error(1,"Couldn't allocate density field\n");
     Cdens_sm_local=(fftwf_complex *)Dens_sm_local;
   }
+  if(TaskFractalD) {
+    Dens_smD_local=fftwf_alloc_real(2*dsize);
+    if(Dens_smD_local==NULL)
+      report_error(1,"Couldn't allocate density field\n");
+    Cdens_smD_local=(fftwf_complex *)Dens_smD_local;
+  }
   if(TaskTidal) {
     Tid_local=my_malloc(6*sizeof(float *));
     Ctid_local=my_malloc(6*sizeof(fftwf_complex *));
diff --git a/src/common.h b/src/common.h
index dccd7ac..e549763 100644
--- a/src/common.h
+++ b/src/common.h
@@ -49,6 +49,10 @@ extern int TaskSmooth;
 extern float *Dens_sm_local;
 extern fftwf_complex *Cdens_sm_local;
 
+extern int TaskFractalD;
+extern float *Dens_smD_local;
+extern fftwf_complex *Cdens_smD_local;
+
 extern int TaskVel;
 extern float **Vel_local;
 
@@ -120,7 +124,7 @@ void write_output(char *prefix);
 void compute_overdensity(ulint np,ulint np_total,float *pos,float *delta,int interp_order);
 void compute_velocity_and_overdensity(ulint np,ulint np_total,float *pos,float *vel,
 				      float *delta,float **velgrid,int interp_order);
-void smooth_density_fourier(float r_smooth);
+void smooth_density_fourier(float r_smooth,int gaus);
 void get_smoothed_density_real(void);
 void get_tidal_field(void);
 void get_linearized_velocity(gad_header head);
diff --git a/src/grid_tools.c b/src/grid_tools.c
index b8ad6ce..cde0d33 100644
--- a/src/grid_tools.c
+++ b/src/grid_tools.c
@@ -635,7 +635,34 @@ void compute_velocity_and_overdensity(ulint np,ulint np_total,float *pos,float *
   }
 }
 
-void smooth_density_fourier(float r_smooth)
+static double window_smooth(double x2, int gaus, int der)
+{
+  if(gaus) {
+    if(der)
+      return (3-x2)*exp(-0.5*x2);
+    else
+      return exp(-0.5*x2);
+  }
+  else { // Top-hat smoothing.
+    if(x2<0.04) {// Taylor expansion below x=0.2 (1E-6 accurate).
+      if(der)
+        return 3-0.5*x2+0.025*x2*x2;
+      else
+        return 1-0.1*x2+0.00357143*x2*x2;
+    }
+    else {
+      double x=sqrt(x2);
+      if(der)
+        return 3*sin(x)/x;
+      else {
+        double j1=(sin(x)-x*cos(x))/x2;
+        return 3*j1/x;
+      }
+    }
+  }
+}
+
+void smooth_density_fourier(float r_smooth, int gaus)
 {
   int iy;
   fftwf_plan plan_tof,plan_tor;
@@ -669,11 +696,15 @@ void smooth_density_fourier(float r_smooth)
       int ix0= ix<=(Ngrid/2) ? ix : ix-Ngrid;
       for(iz=0;iz<Ngrid/2+1;iz++) {
 	float x2=x_smooth2*(ix0*ix0+iy0*iy0+iz*iz);
-	float sm=(float)(exp(-0.5*x2));
+        float sm=(float)(window_smooth((double)x2, gaus, 0));
 	long index=iz+(Ngrid/2+1)*((long)(ix+Ngrid*iy));
 
 	Cdens_local[index]*=norm;
 	Cdens_sm_local[index]=Cdens_local[index]*sm;
+        if(TaskFractalD) {
+          float smD=(float)(window_smooth((double)x2, gaus, 1));
+          Cdens_smD_local[index]=Cdens_local[index]*smD;
+        }
       }
     }
   }
@@ -689,6 +720,7 @@ void get_smoothed_density_real(void)
 {
   fftwf_plan plan_tor;
 
+  // Smoothed overdensity
 #ifdef _DEBUG
   printf("Node %d Planning\n",NodeThis);
 #endif //_DEBUG
@@ -699,6 +731,20 @@ void get_smoothed_density_real(void)
 #endif //_DEBUG
   fftwf_execute(plan_tor);
   fftwf_destroy_plan(plan_tor);
+
+  if(TaskFractalD) {
+    // Smoothed derivative
+#ifdef _DEBUG
+    printf("Node %d Planning\n",NodeThis);
+#endif //_DEBUG
+    plan_tor=fftwf_mpi_plan_dft_c2r_3d(Ngrid,Ngrid,Ngrid,Cdens_smD_local,Dens_smD_local,
+                                       MPI_COMM_WORLD,FFTW_ESTIMATE | FFTW_MPI_TRANSPOSED_IN);
+#ifdef _DEBUG
+    printf("Node %d Transforming back\n",NodeThis);
+#endif //_DEBUG
+    fftwf_execute(plan_tor);
+    fftwf_destroy_plan(plan_tor);
+  }
 }
 
 static void get_tidal_field_fd(void)
diff --git a/src/io.c b/src/io.c
index 5c7524d..9d028ff 100644
--- a/src/io.c
+++ b/src/io.c
@@ -441,6 +441,25 @@ void write_output(char *prefix)
     fclose(fo);
   }
 
+  if(TaskFractalD) {
+    sprintf(fname,"%s_dens_smD.%04d",prefix,NodeThis);
+    num_grids=1;
+    fo=my_fopen(fname,"wb");
+    my_fwrite(&(num_grids),sizeof(int),1,fo);
+    my_fwrite(&(Ngrid),sizeof(int),1,fo);
+    my_fwrite(&(Nx_here),sizeof(int),1,fo);
+    for(ix=0;ix<Nx_here;ix++) {
+      int iy;
+      int ix_id=Ix0_here+ix;
+      my_fwrite(&(ix_id),sizeof(int),1,fo);
+      for(iy=0;iy<Ngrid;iy++) {
+	lint index0=2*(Ngrid/2+1)*((lint)(iy+ix*Ngrid));
+	my_fwrite(&(Dens_smD_local[index0]),sizeof(float),Ngrid,fo);
+      }
+    }
+    fclose(fo);
+  }
+
   if(TaskLinvel) {
     sprintf(fname,"%s_linvel.%04d",prefix,NodeThis);
     num_grids=3;
diff --git a/src/main.c b/src/main.c
index 3d78107..0ae4e66 100644
--- a/src/main.c
+++ b/src/main.c
@@ -6,6 +6,7 @@ int main(int argc,char **argv)
   char fnameOut[256]="output";
   char interp_method[256]="NGP";
   float r_smooth=0.;
+  int gaus=1;
   gad_header head;
   int input_format=1;
   int ii,interp_order=0;
@@ -21,6 +22,7 @@ int main(int argc,char **argv)
     else if(!strcmp(*c,"-interp")) sprintf(interp_method,"%s",*++c);
     else if(!strcmp(*c,"-diag_tidal")) TaskTidalDiag=1;
     else if(!strcmp(*c,"-use_finite_differences")) UseFD=1;
+    else if(!strcmp(*c,"-tophat_smooth")) gaus=0;
     else if(!strcmp(*c,"-do")) {
       int ich=0;
       char d=(*++c)[ich];
@@ -33,8 +35,10 @@ int main(int argc,char **argv)
 	  TaskLinvel=1;
 	else if(d=='u')
 	  TaskNlvel=1;
+	else if(d=='d')
+	  TaskFractalD=1;
 	else
-	  fprintf(stderr,"Unknown task %c. Supported : v, t, l\n",d);
+	  fprintf(stderr,"Unknown task %c. Supported : v, t, l, u\n",d);
 	d=(*c)[++ich];
       }
     }
@@ -50,9 +54,11 @@ int main(int argc,char **argv)
 	fprintf(stderr,"  -smooth     -> smoothing length (in same units as simulation box)\n");
 	fprintf(stderr,"  -interp     -> particle interpolation scheme: NGP, CIC or TSC\n");
 	fprintf(stderr,"  -do         -> concatenate tasks: velocity (v), tidal field (t),"
-		" linearized velocity (l), non-linear velocity (u)\n");
+		" linearized velocity (l), non-linear velocity (u),"
+                " local fractal dim (d)\n");
 	fprintf(stderr,"  -diag_tidal -> diagonalize tidal tensor\n");
 	fprintf(stderr,"  -use_finite_differences -> use central FDs for derivatives\n");
+	fprintf(stderr,"  -tophat_smooth -> Use top-hat smoothing\n");
 	fprintf(stderr,"  -h          -> this help\n\n");
 	return 0;
       }
@@ -64,7 +70,7 @@ int main(int argc,char **argv)
   }
   if(!TaskTidal)
     TaskTidalDiag=0;
-  if(r_smooth>0 || TaskTidal+TaskLinvel+TaskNlvel)
+  if(r_smooth>0 || TaskTidal+TaskLinvel+TaskNlvel+TaskFractalD)
     TaskSmooth=1;
   if(strcmp(interp_method,"NGP") && strcmp(interp_method,"CIC") && strcmp(interp_method,"TSC"))
     report_error(1,"Unknown interpolation scheme %s. Use NGP, CIC or TSC\n",interp_method);
@@ -75,6 +81,10 @@ int main(int argc,char **argv)
     printf("   Output prefix : %s\n",fnameOut);
     printf("   Ngrid : %d\n",Ngrid);
     printf("   R_smooth : %.2lf \n",r_smooth);
+    if(gaus)
+      printf("     (Gaussian smoothing)\n");
+    else
+      printf("     (Top-hat smoothing)\n");
     printf("   Interpolation scheme : %s\n",interp_method);
     printf("   Tasks to do:");
     if(TaskVel)
@@ -89,6 +99,8 @@ int main(int argc,char **argv)
       printf(" LV");
     if(TaskNlvel)
       printf(" NLV");
+    if(TaskFractalD)
+      printf(" FD");
     printf("\n");
     printf("   Use finite differences : %d\n",UseFD);
     printf("\n");
@@ -139,7 +151,7 @@ int main(int argc,char **argv)
   if(TaskSmooth) {
     if(NodeThis==0)
       printf("* Smoothing density field\n");
-    smooth_density_fourier(r_smooth);
+    smooth_density_fourier(r_smooth, gaus);
     if(NodeThis==0)
       printf("\n");
   }