diff --git a/.gitignore b/.gitignore
index 6257730b3..bbadc0563 100644
--- a/.gitignore
+++ b/.gitignore
@@ -28,3 +28,4 @@ Hirep-git.code-workspace
 **/clangd
 .ycm_extra_conf.py*
 compile_commands.json
+Spectrum/measure_spectrum
diff --git a/Include/global.h b/Include/global.h
index 0d3b2e90b..b7da350a5 100644
--- a/Include/global.h
+++ b/Include/global.h
@@ -143,6 +143,15 @@ GLB_VAR(suNg_field_flt,*u_gauge_flt,=NULL);
 GLB_VAR(suNf_field,*u_gauge_f,=NULL);
 GLB_VAR(suNg_field,*u_gauge_s,=NULL);
 GLB_VAR(suNf_field_flt,*u_gauge_f_flt,=NULL);
+
+/*APE smearing*/
+GLB_VAR(suNg_field,*u_gauge_APE,=NULL);
+GLB_VAR(suNg_scalar_field,*u_scalar_APE,=NULL);
+GLB_VAR(suNg_field_flt,*u_gauge_APE_flt,=NULL);
+GLB_VAR(suNf_field,*u_gauge_APE_f,=NULL);
+GLB_VAR(suNg_field,*u_gauge_APE_s,=NULL);
+GLB_VAR(suNf_field_flt,*u_gauge_APE_f_flt,=NULL);
+
 #if defined(GAUGE_SPN) && defined(REPR_FUNDAMENTAL)
 GLB_VAR(suNffull_field,*cl_term,=NULL);
 GLB_VAR(suNffull_field,*cl_force,=NULL);
@@ -161,6 +170,14 @@ GLB_VAR(int,gauge_field_active,=0); // whether gauge field interactions is activ
 #define pu_gauge_f(ix,mu) ((u_gauge_f->ptr)+coord_to_index(ix,mu))
 #define pu_gauge_f_flt(ix,mu) ((u_gauge_f_flt->ptr)+coord_to_index(ix,mu))
 
+/*APE smearing*/
+#define pu_gauge_APE(ix,mu) ((u_gauge_APE->ptr)+coord_to_index(ix,mu))
+#define pu_scalar_APE(ix) ((u_scalar_APE->ptr)+ix)
+#define pu_gauge_APE_flt(ix,mu) ((u_gauge_APE_flt->ptr)+coord_to_index(ix,mu))
+#define pu_gauge_APE_f(ix,mu) ((u_gauge_APE_f->ptr)+coord_to_index(ix,mu))
+#define pu_gauge_APE_f_flt(ix,mu) ((u_gauge_APE_f_flt->ptr)+coord_to_index(ix,mu))
+#define pu_gauge_tmp(ix,mu) ((u_gauge_tmp->ptr)+coord_to_index(ix,mu))
+
 /* input parameters */
 #include "input_par.h"
 GLB_VAR(input_glb,glb_var,=init_input_glb(glb_var));
diff --git a/Include/observables.h b/Include/observables.h
index d7a9bef87..2395b4971 100644
--- a/Include/observables.h
+++ b/Include/observables.h
@@ -296,4 +296,32 @@ void init_triplet_discon_correlators();
 
 void ff_observables();
 
+//APE smearing
+double plaq_APE(int ix,int mu,int nu);
+void cplaq_APE(complex *ret,int ix,int mu,int nu);
+double avr_spatial_plaquette_APE();
+void full_plaquette_APE();
+void polyakov_APE();
+
+/* source smearing*/
+void smearing_function(spinor_field *source, int tau, double epsilon);
+void smearing_function_with_APE(spinor_field *source, int tau, double epsilon);
+void smearing_function_volume(spinor_field *source, double epsilon);
+void smearing_function_volume_with_APE(spinor_field *source, double epsilon);
+void create_smeared_source(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear);
+void create_smeared_source_with_APE(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear);
+void create_noise_source_equal_eo_smeared_with_APE(spinor_field *source, double epsilon, int Nsmear);
+void create_noise_source_equal_eo_smeared(spinor_field *source, double epsilon, int Nsmear);
+
+/* Sink smearing*/
+void smeared_propagator(spinor_field* psi, int nm , double epsilon);
+void smeared_propagator_with_APE(spinor_field* psi, int nm , double epsilon);
+void smeared_propagator_volume(spinor_field* psi, int nm , double epsilon);
+void smeared_propagator_volume_with_APE(spinor_field* psi, int nm , double epsilon);
+
+/* Measurements*/
+void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source_max, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff);
+void measure_smearing_ss(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff);
+
+
 #endif
diff --git a/Include/representation.h b/Include/representation.h
index 328d8504b..a7b88d129 100644
--- a/Include/representation.h
+++ b/Include/representation.h
@@ -21,4 +21,6 @@ void _group_represent2_flt(suNf_flt* v, suNg_flt *u);
 void represent_gauge_field();
 void represent_gauge_field_measure();
 
+/* APE smearing*/
+void represent_gauge_field_APE();
 #endif
diff --git a/Include/spectrum.h b/Include/spectrum.h
index 7a77cb2e0..d1021515a 100644
--- a/Include/spectrum.h
+++ b/Include/spectrum.h
@@ -19,6 +19,8 @@
 
 void measure_spectrum_semwall(int nm, double* m, int nhits,int conf_num, double precision,storage_switch swc, data_storage_array **ret);
 void measure_spectrum_discon_semwall(int nm, double* m, int nhits,int conf_num, double precision,storage_switch swc, data_storage_array **ret);
+void measure_spectrum_discon_semwall_smeared(int nm, double* m, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double APE_epsilon, int APE_N, int N_diff);
+void measure_spectrum_discon_semwall_smeared_single_inversion(int nm, double* m, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double APE_epsilon, int APE_N, int N_diff);
 void measure_spectrum_discon_gfwall(int nm, double* m, int conf_num, double precision,storage_switch swc, data_storage_array **ret);
 void measure_spectrum_discon_volume(int nm, double* m, int conf_num, double precision, int dil,storage_switch swc, data_storage_array **ret);
 void measure_spectrum_gfwall(int nm, double* m, int conf_num, double precision,storage_switch swc, data_storage_array **ret);
diff --git a/Include/utils.h b/Include/utils.h
index 98d3d9437..fdda4d801 100644
--- a/Include/utils.h
+++ b/Include/utils.h
@@ -44,6 +44,7 @@ void init_plaq_open_BCs(double *plaq_weight, double *rect_weight, double ct, dou
 
 void free_BCs();
 void apply_BCs_on_represented_gauge_field();
+void apply_BCs_on_represented_gauge_field_APE();
 void apply_BCs_on_fundamental_gauge_field();
 void apply_BCs_on_momentum_field(suNg_av_field *force);
 void apply_BCs_on_spinor_field(spinor_field *sp);
@@ -176,4 +177,15 @@ int reserve_wrk_space_with_pointers(suNg_field **g_wrk_out, int **i_up_wrk_out,
 void release_wrk_space(int id_release);
 void free_wrk_space();
 
+/* APE smearing */
+void APE_smearing(double smear_val, int Nsmear);
+void assign_ud2u_APE_f(void);
+#ifdef GAUGE_SPN
+void cooling_SPN(suNg* g_out, suNg* g_in, suNg* staple, int cooling);
+void subgrb(int i1col, int i2col, suNgfull* B11, suNgfull* C11);
+void vmxsu2(int i1, int i2, suNgfull* A, complex B[4]);
+void subgrb_tau(int n1, int n2, suNgfull* B11, suNgfull* C11);
+void vmxsu2_tau(int n1, int n2, suNgfull* A, complex B[4]);
+#endif
+
 #endif
diff --git a/LibHR/Observables/avr_plaquette.c b/LibHR/Observables/avr_plaquette.c
index 0a47117b1..a6e9ba4f1 100644
--- a/LibHR/Observables/avr_plaquette.c
+++ b/LibHR/Observables/avr_plaquette.c
@@ -416,3 +416,214 @@ double complex avr_plaquette_wrk()
 
   return pa;
 }
+
+
+double plaq_APE(int ix,int mu,int nu){
+  /**
+   * @brief Computes the value of the plaquette at the given position in the mu-nu plane and return only the real part.
+   *        Modified from `plaq(int ix, int mu, int nu)` by changing `pu_gauge` to `pu_gauge_APE`
+   * 
+   * @param ix The position of the starting site.
+   * @param mu The first direction of the plaquette.
+   * @param nu The second direction of the plaquette.
+   * 
+   * @return The real part of the computed plaquette value, possibly weighted if PLAQ_WEIGHTS is defined.
+   */
+  int iy,iz;
+  double p;
+  suNg *v1,*v2,*v3,*v4,w1,w2,w3;
+
+  iy=iup(ix,mu);
+  iz=iup(ix,nu);
+
+  v1=pu_gauge_APE(ix,mu);
+  v2=pu_gauge_APE(iy,nu);
+  v3=pu_gauge_APE(iz,mu);
+  v4=pu_gauge_APE(ix,nu);
+
+  _suNg_times_suNg(w1,(*v1),(*v2));
+  _suNg_times_suNg(w2,(*v4),(*v3));
+  _suNg_times_suNg_dagger(w3,w1,w2);
+
+  _suNg_trace_re(p,w3);
+
+#ifdef PLAQ_WEIGHTS
+  if(plaq_weight==NULL) return p;
+  return plaq_weight[ix*16+mu*4+nu]*p;
+#else
+  return p;
+#endif
+}
+
+
+void cplaq_APE(complex *ret,int ix,int mu,int nu){
+  /**
+   * @brief Computes the value of the plaquette at the given position in the mu-nu plane.
+   *        Modified from `cplaq(complex *ret,int ix,int mu,int nu)` by changing `pu_gauge` to `pu_gauge_APE`
+   * 
+   * @param ret A pointer to store the resulting complex plaquette value.
+   * @param ix The position of the starting site.
+   * @param mu The first direction of the plaquette.
+   * @param nu The second direction of the plaquette.
+   * 
+   * @note If PLAQ_WEIGHTS is defined, the result may include a weighting factor.
+   */
+  int iy,iz;
+  suNg *v1,*v2,*v3,*v4,w1,w2,w3;
+  double tmpre = 0.;
+  double tmpim = 0.;
+
+  iy=iup(ix,mu);
+  iz=iup(ix,nu);
+
+  v1=pu_gauge_APE(ix,mu);
+  v2=pu_gauge_APE(iy,nu);
+  v3=pu_gauge_APE(iz,mu);
+  v4=pu_gauge_APE(ix,nu);
+
+  _suNg_times_suNg(w1,(*v1),(*v2));
+  _suNg_times_suNg(w2,(*v4),(*v3));
+  _suNg_times_suNg_dagger(w3,w1,w2);
+
+  _suNg_trace_re(tmpre, w3);
+#ifdef GAUGE_SON
+  _suNg_trace_im(tmpim, w3);
+#else
+  *ret = tmpre + I * tmpim;
+#endif
+
+#ifdef PLAQ_WEIGHTS
+  if(plaq_weight!=NULL) {
+    *ret *= plaq_weight[ix * 16 + mu * 4 + nu];;
+  }
+#endif
+
+}
+
+
+double avr_spatial_plaquette_APE(){
+  /**
+   * @brief Computes the average (over the whole lattice) plaquette value over all spatial directions 
+   *        for the APE-smeared gauge links.
+   *        Modified from `avr_spatial_plaquette()` by changing `plaq()` to `plaq_APE()`
+   * 
+   */
+  double pa=0.;
+
+  _PIECE_FOR(&glattice,ixp) {
+    if(ixp==glattice.inner_master_pieces) {
+      _OMP_PRAGMA( master )
+      /* wait for gauge field to be transfered */
+      complete_gf_sendrecv(u_gauge_APE);
+      _OMP_PRAGMA( barrier )
+    }
+    _SITE_FOR_SUM(&glattice,ixp,ix,pa) {
+      
+      pa+=plaq_APE(ix,2,1);
+      pa+=plaq_APE(ix,3,1);
+      pa+=plaq_APE(ix,3,2);
+    }
+  }
+
+  global_sum(&pa, 1);
+
+#ifdef BC_T_OPEN
+    pa /= 3.0*NG*GLB_VOLUME*(GLB_T-1)/GLB_T;
+#else
+    pa /= 3.0*NG*GLB_VOLUME;
+#endif
+
+  return pa;
+
+}
+
+void full_plaquette_APE(){
+  /**
+   * @brief Computes and prints the average plaquette value over all directions, 
+   *        including the temporal direction, for the APE-smeared gauge links.
+   *        Modified from `ull_plaquette())` by changing `cplaq(...)` to `cplaq_APE(...)`
+   * 
+   * @note APE smearing is not applied in the temporal direction. 
+   *       As a result, plaquette values involving the temporal direction are zero.
+   */
+
+    static double complex pa[6];
+    static double complex r0;
+    static double complex r1;
+    static double complex r2;
+    static double complex r3;
+    static double complex r4;
+    static double complex r5;
+
+    _OMP_PRAGMA(single)
+    {
+      r0 = 0.;
+      r1 = 0.;
+      r2 = 0.;
+      r3 = 0.;
+      r4 = 0.;
+      r5 = 0.;
+    }
+
+    _PIECE_FOR(&glattice,ixp)
+    {
+        if(ixp == glattice.inner_master_pieces)
+        {
+            _OMP_PRAGMA( master )
+            /* wait for gauge field to be transfered */
+            complete_gf_sendrecv(u_gauge_APE);
+            _OMP_PRAGMA( barrier )
+        }
+
+        _SITE_FOR_SUM(&glattice,ixp,ix,r0re,r0im,r1re,r1im,r2re,r2im,r3re,r3im,r4re,r4im,r5re,r5im)
+        {
+            complex tmp;
+
+            cplaq_APE(&tmp,ix,1,0);
+            r0 += tmp;
+
+            cplaq_APE(&tmp,ix,2,0);
+            r1 += tmp;
+
+            cplaq_APE(&tmp,ix,2,1);
+            r2 += tmp;
+
+            cplaq_APE(&tmp,ix,3,0);
+            r3 += tmp;
+
+            cplaq_APE(&tmp,ix,3,1);
+            r4 += tmp;
+
+            cplaq_APE(&tmp,ix,3,2);
+            r5 += tmp;
+
+        }
+    }
+
+    _OMP_PRAGMA(single)
+    {
+      pa[0] = r0;
+      pa[1] = r1;
+      pa[2] = r2;
+      pa[3] = r3;
+      pa[4] = r4;
+      pa[5] = r5;
+    }
+
+    global_sum((double*)pa,12);
+    for(int k = 0; k < 6; k++)
+    {
+#ifdef BC_T_OPEN
+        pa[k] /= NG*GLB_VOLUME*(GLB_T-1)/GLB_T;
+#else
+        pa[k] /= NG*GLB_VOLUME;
+#endif
+    }
+
+    lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 1, 0, creal(pa[0]), cimag(pa[0]));
+    lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 2, 0, creal(pa[1]), cimag(pa[1]));
+    lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 2, 1, creal(pa[2]), cimag(pa[2]));
+    lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 3, 0, creal(pa[3]), cimag(pa[3]));
+    lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 3, 1, creal(pa[4]), cimag(pa[4]));
+    lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 3, 2, creal(pa[5]), cimag(pa[5]));
+}
\ No newline at end of file
diff --git a/LibHR/Observables/measure_mesons.c b/LibHR/Observables/measure_mesons.c
index 378eb0622..d3e2c0108 100644
--- a/LibHR/Observables/measure_mesons.c
+++ b/LibHR/Observables/measure_mesons.c
@@ -623,3 +623,406 @@ void print_mesons(meson_observable *mo, double norm, int conf, int nm, double *m
   print_corr(mo, lt, conf, nm, mass, label, n_mom);
   zero_corrs(mo);
 }
+
+
+/* Wuppertal smearing*/
+
+void smeared_propagator(spinor_field* psi, int nm , double epsilon){
+    /**
+     * @brief Smears the propagator field `psi` with the Wuppertal smearingwith parameter `epsilon` and one iteration, without APE smearing.
+     * The smearing function is defined as Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf.
+     * See also disucssion in the paragraph below Eq. (21) in the same paper.
+     * 
+     * @param psi The propagator field to be smeared, $S$ in the paper.
+     * @param nm The number of masses
+     * @param epsilon The smearing step size, $\epsilon$ in Eq. (9).
+     */
+
+    int i,t, x, y, z, ix, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down;
+    double norm_factor = 1./(1.+6.*epsilon);
+    suNf_propagator sp_OG, sp_smeared, sp_tmp, sp_right, sp_left, sp_front, sp_back, sp_up, sp_down;
+    _propagator_zero(sp_smeared);
+    
+    spinor_field* smeared_psi = alloc_spinor_field_f(4*nm*NF,&glattice);
+    
+    for(i=0; i<nm; i++) {
+        for (t=0; t<T; t++) {
+            for (x=0; x<X; x++) for (y=0; y<Y; y++) for (z=0; z<Z; z++) {
+                
+                ix=ipt(t,x,y,z);
+                ix_right = iup(ix,1);
+                ix_left  = idn(ix,1);
+                ix_front = iup(ix,2);
+                ix_back  = idn(ix,2);
+                ix_up    = iup(ix,3);
+                ix_down  = idn(ix,3);
+                
+                for (int a=0;a<NF;++a){
+                    for (int beta=0;beta<4;beta++){
+                        _propagator_assign(sp_OG, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix),a,beta);
+                        _propagator_assign(sp_right, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix_right),a,beta);
+                        _propagator_assign(sp_left, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix_left),a,beta);
+                        _propagator_assign(sp_front, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix_front),a,beta);
+                        _propagator_assign(sp_back, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix_back),a,beta);
+                        _propagator_assign(sp_up, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix_up),a,beta);
+                        _propagator_assign(sp_down, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix_down),a,beta);
+                    }
+                }
+                
+                _propagator_zero(sp_smeared);
+                _propagator_mul_assign(sp_OG, norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_OG);
+
+                _suNf_prop_multiply(sp_tmp, *pu_gauge_f(ix,1), sp_right);
+                _propagator_mul_assign(sp_tmp, epsilon*norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_tmp);
+    
+                _suNf_inverse_prop_multiply(sp_tmp, *pu_gauge_f(ix_left,1), sp_left);
+                _propagator_mul_assign(sp_tmp, epsilon*norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_tmp);
+                
+                _suNf_prop_multiply(sp_tmp, *pu_gauge_f(ix,2), sp_front);
+                _propagator_mul_assign(sp_tmp, epsilon*norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_tmp);
+                
+                _suNf_inverse_prop_multiply(sp_tmp, *pu_gauge_f(ix_back,2), sp_back);
+                _propagator_mul_assign(sp_tmp, epsilon*norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_tmp);
+                
+                _suNf_prop_multiply(sp_tmp, *pu_gauge_f(ix,3), sp_up);
+                _propagator_mul_assign(sp_tmp, epsilon*norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_tmp);
+                
+                _suNf_inverse_prop_multiply(sp_tmp, *pu_gauge_f(ix_down,3), sp_down);
+                _propagator_mul_assign(sp_tmp, epsilon*norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_tmp);
+                
+                for (int a=0;a<NF;a++) for (int b=0;b<NF;b++){
+                    for (int alpha=0;alpha<4;alpha++) for (int beta=0;beta<4;beta++){
+                        
+                        _FIELD_AT(&smeared_psi[a*4*nm+beta*nm+i],ix)->c[alpha].c[b] = sp_smeared.c[b].c[alpha].c[beta].c[a];
+                    
+                    }
+                }
+               
+            }
+        }
+    }
+    for (int i=0;i<4*nm*NF;i++){
+        spinor_field_copy_f(psi + i, smeared_psi  + i);
+    }
+    for (int i=0;i<4*nm*NF;i++){
+        start_sf_sendrecv(psi + i);
+        complete_sf_sendrecv(psi + i);
+    }
+    free_spinor_field_f(smeared_psi);
+}
+
+void smeared_propagator_with_APE(spinor_field* psi, int nm , double epsilon){
+  /**
+     * @brief Smears the propagator field `psi` with the Wuppertal smearing with parameter `epsilon` and one iteration.
+     * The smearing function is defined as Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf.
+     * See also disucssion in the paragraph below Eq. (21) in the same paper.
+     * 
+     * @param psi The propagator field to be smeared, $S$ in the paper.
+     * @param nm The number of masses
+     * @param epsilon The smearing step size, $\epsilon$ in Eq. (9).
+     * 
+     * The gauge field is smeared with APE smearing before the Wuppertal smearing of the propagator field.
+     */
+
+    
+    int i,t, x, y, z, ix, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down;
+    double norm_factor = 1./(1.+6.*epsilon);
+    suNf_propagator sp_OG, sp_smeared, sp_tmp, sp_right, sp_left, sp_front, sp_back, sp_up, sp_down;
+    _propagator_zero(sp_smeared);
+    
+    spinor_field* smeared_psi = alloc_spinor_field_f(4*nm*NF, &glattice);
+    
+    for(i=0; i<nm; i++) {
+        for (t=0; t<T; t++) {
+            for (x=0; x<X; x++) for (y=0; y<Y; y++) for (z=0; z<Z; z++) {
+                
+                ix=ipt(t,x,y,z);
+                ix_right = iup(ix,1);
+                ix_left  = idn(ix,1);
+                ix_front = iup(ix,2);
+                ix_back  = idn(ix,2);
+                ix_up    = iup(ix,3);
+                ix_down  = idn(ix,3);
+                
+                for (int a=0;a<NF;++a){
+                    for (int beta=0;beta<4;beta++){
+                        _propagator_assign(sp_OG, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix),a,beta);
+                        _propagator_assign(sp_right, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix_right),a,beta);
+                        _propagator_assign(sp_left, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix_left),a,beta);
+                        _propagator_assign(sp_front, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix_front),a,beta);
+                        _propagator_assign(sp_back, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix_back),a,beta);
+                        _propagator_assign(sp_up, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix_up),a,beta);
+                        _propagator_assign(sp_down, *_FIELD_AT(&psi[a*4*nm+beta*nm+i],ix_down),a,beta);
+                    }
+                }
+                
+                _propagator_zero(sp_smeared);
+                _propagator_mul_assign(sp_OG, norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_OG);
+
+                _suNf_prop_multiply(sp_tmp, *pu_gauge_APE_f(ix,1), sp_right);
+                _propagator_mul_assign(sp_tmp, epsilon*norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_tmp);
+    
+                _suNf_inverse_prop_multiply(sp_tmp, *pu_gauge_APE_f(ix_left,1), sp_left);
+                _propagator_mul_assign(sp_tmp, epsilon*norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_tmp);
+                
+                _suNf_prop_multiply(sp_tmp, *pu_gauge_APE_f(ix,2), sp_front);
+                _propagator_mul_assign(sp_tmp, epsilon*norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_tmp);
+                
+                _suNf_inverse_prop_multiply(sp_tmp, *pu_gauge_APE_f(ix_back,2), sp_back);
+                _propagator_mul_assign(sp_tmp, epsilon*norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_tmp);
+                
+                _suNf_prop_multiply(sp_tmp, *pu_gauge_APE_f(ix,3), sp_up);
+                _propagator_mul_assign(sp_tmp, epsilon*norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_tmp);
+                
+                _suNf_inverse_prop_multiply(sp_tmp, *pu_gauge_APE_f(ix_down,3), sp_down);
+                _propagator_mul_assign(sp_tmp, epsilon*norm_factor);
+                _propagator_add(sp_smeared,sp_smeared, sp_tmp);
+                
+                //FILED_prop_assign(smeared_psi, sp_smeared, ix, i, nm);
+                
+                for (int a=0;a<NF;a++) for (int b=0;b<NF;b++){
+                    for (int alpha=0;alpha<4;alpha++) for (int beta=0;beta<4;beta++){
+                        
+                        _FIELD_AT(&smeared_psi[a*4*nm+beta*nm+i],ix)->c[alpha].c[b] = sp_smeared.c[b].c[alpha].c[beta].c[a];
+                    
+                    }
+                }
+            }
+        }
+    }
+    for (int i=0;i<4*nm*NF;i++){
+        spinor_field_copy_f(psi + i, smeared_psi  + i);
+    }
+    for (int i=0;i<4*nm*NF;i++){
+        start_sf_sendrecv(psi + i);
+        complete_sf_sendrecv(psi + i);
+    }
+    free_spinor_field_f(smeared_psi);
+}
+
+
+// FZ (2024): Smeared propagator with a single source
+void smeared_propagator_volume(spinor_field* psi, int nm , double epsilon){
+    int i,t, x, y, z, ix, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down;
+    double norm_factor = 1./(1.+6.*epsilon);
+    suNf_spin_matrix sp_OG, sp_smeared, sp_tmp, sp_right, sp_left, sp_front, sp_back, sp_up, sp_down;
+    _spinmatrix_zero(sp_smeared);
+    
+    spinor_field* smeared_psi = alloc_spinor_field_f(4*nm,&glattice);
+    
+    for(i=0; i<nm; i++) {
+        for (t=0; t<T; t++) {
+            for (x=0; x<X; x++) for (y=0; y<Y; y++) for (z=0; z<Z; z++) {
+                
+                ix=ipt(t,x,y,z);
+                ix_right = iup(ix,1);
+                ix_left  = idn(ix,1);
+                ix_front = iup(ix,2);
+                ix_back  = idn(ix,2);
+                ix_up    = iup(ix,3);
+                ix_down  = idn(ix,3);
+                
+                // The code identified the rows of the spinmatrix as the spinors. I will copy this
+                // convention for this function.
+                for (int beta=0;beta<4;beta++){
+                    _spinmatrix_assign_row(sp_OG, *_FIELD_AT(&psi[beta*nm+i],ix),beta);
+                    _spinmatrix_assign_row(sp_right, *_FIELD_AT(&psi[beta*nm+i],ix_right),beta);
+                    _spinmatrix_assign_row(sp_left, *_FIELD_AT(&psi[beta*nm+i],ix_left),beta);
+                    _spinmatrix_assign_row(sp_front, *_FIELD_AT(&psi[beta*nm+i],ix_front),beta);
+                    _spinmatrix_assign_row(sp_back, *_FIELD_AT(&psi[beta*nm+i],ix_back),beta);
+                    _spinmatrix_assign_row(sp_up, *_FIELD_AT(&psi[beta*nm+i],ix_up),beta);
+                    _spinmatrix_assign_row(sp_down, *_FIELD_AT(&psi[beta*nm+i],ix_down),beta);
+                }
+                
+                _spinmatrix_zero(sp_smeared);
+
+                // No deicated function for a spin matrix esists. Loop over rows in the spinmatrix
+                // which are spinors and use that method.
+                for (int beta=0;beta<4;beta++){
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],norm_factor, sp_OG.c[beta]);
+                }
+
+                // For the multiplication with a colour matrix we need to loop over every entry in the
+                // spin matrix, i.e. two loops are needed. Assignment and addition can be done in the 
+                // outermost loop since deicated functions for spinors exist.
+                for (int beta=0; beta<4; beta++){
+                  for (int alpha=0; alpha<4; alpha++){
+                   _suNf_multiply(sp_tmp.c[beta].c[alpha], *pu_gauge_f(ix,1), sp_right.c[beta].c[alpha]);
+                  }
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],epsilon*norm_factor, sp_tmp.c[beta]);
+                }
+                
+                // Repeat the same approach for the inverse multiplication and smearing directions
+                for (int beta=0; beta<4; beta++){
+                  for (int alpha=0; alpha<4; alpha++){
+                   _suNf_inverse_multiply(sp_tmp.c[beta].c[alpha], *pu_gauge_f(ix_left,1), sp_left.c[beta].c[alpha]);
+                  }
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],epsilon*norm_factor, sp_tmp.c[beta]);
+                }
+
+                for (int beta=0; beta<4; beta++){
+                  for (int alpha=0; alpha<4; alpha++){
+                   _suNf_multiply(sp_tmp.c[beta].c[alpha], *pu_gauge_f(ix,2), sp_front.c[beta].c[alpha]);
+                  }
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],epsilon*norm_factor, sp_tmp.c[beta]);
+                }
+                
+                for (int beta=0; beta<4; beta++){
+                  for (int alpha=0; alpha<4; alpha++){
+                   _suNf_inverse_multiply(sp_tmp.c[beta].c[alpha], *pu_gauge_f(ix_back,2), sp_back.c[beta].c[alpha]);
+                  }
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],epsilon*norm_factor, sp_tmp.c[beta]);
+                }                
+
+                for (int beta=0; beta<4; beta++){
+                  for (int alpha=0; alpha<4; alpha++){
+                   _suNf_multiply(sp_tmp.c[beta].c[alpha], *pu_gauge_f(ix,3), sp_up.c[beta].c[alpha]);
+                  }
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],epsilon*norm_factor, sp_tmp.c[beta]);
+                }
+                
+                for (int beta=0; beta<4; beta++){
+                  for (int alpha=0; alpha<4; alpha++){
+                   _suNf_inverse_multiply(sp_tmp.c[beta].c[alpha], *pu_gauge_f(ix_down,3), sp_down.c[beta].c[alpha]);
+                  }
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],epsilon*norm_factor, sp_tmp.c[beta]);
+                }           
+                
+                for (int b=0;b<NF;b++){
+                    for (int alpha=0;alpha<4;alpha++) for (int beta=0;beta<4;beta++){
+                        _FIELD_AT(&smeared_psi[beta*nm+i],ix)->c[alpha].c[b] = sp_smeared.c[beta].c[alpha].c[b];
+                        _FIELD_AT(&smeared_psi[beta*nm+i],ix)->c[alpha].c[b] = sp_smeared.c[beta].c[alpha].c[b];
+                    }
+                }
+               
+            }
+        }
+    }
+    for (int i=0;i<4*nm;i++){
+        spinor_field_copy_f(psi + i, smeared_psi  + i);
+    }
+    for (int i=0;i<4*nm;i++){
+        start_sf_sendrecv(psi + i);
+        complete_sf_sendrecv(psi + i);
+    }
+    free_spinor_field_f(smeared_psi);
+}
+
+void smeared_propagator_volume_with_APE(spinor_field* psi, int nm , double epsilon){
+    int i,t, x, y, z, ix, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down;
+    double norm_factor = 1./(1.+6.*epsilon);
+    suNf_spin_matrix sp_OG, sp_smeared, sp_tmp, sp_right, sp_left, sp_front, sp_back, sp_up, sp_down;
+    _spinmatrix_zero(sp_smeared);
+    
+    spinor_field* smeared_psi = alloc_spinor_field_f(4*nm,&glattice);
+    
+    for(i=0; i<nm; i++) {
+        for (t=0; t<T; t++) {
+            for (x=0; x<X; x++) for (y=0; y<Y; y++) for (z=0; z<Z; z++) {
+                
+                ix=ipt(t,x,y,z);
+                ix_right = iup(ix,1);
+                ix_left  = idn(ix,1);
+                ix_front = iup(ix,2);
+                ix_back  = idn(ix,2);
+                ix_up    = iup(ix,3);
+                ix_down  = idn(ix,3);
+                
+                // The code identified the rows of the spinmatrix as the spinors. I will copy this
+                // convention for this function.
+                for (int beta=0;beta<4;beta++){
+                    _spinmatrix_assign_row(sp_OG, *_FIELD_AT(&psi[beta*nm+i],ix),beta);
+                    _spinmatrix_assign_row(sp_right, *_FIELD_AT(&psi[beta*nm+i],ix_right),beta);
+                    _spinmatrix_assign_row(sp_left, *_FIELD_AT(&psi[beta*nm+i],ix_left),beta);
+                    _spinmatrix_assign_row(sp_front, *_FIELD_AT(&psi[beta*nm+i],ix_front),beta);
+                    _spinmatrix_assign_row(sp_back, *_FIELD_AT(&psi[beta*nm+i],ix_back),beta);
+                    _spinmatrix_assign_row(sp_up, *_FIELD_AT(&psi[beta*nm+i],ix_up),beta);
+                    _spinmatrix_assign_row(sp_down, *_FIELD_AT(&psi[beta*nm+i],ix_down),beta);
+                }
+                
+                _spinmatrix_zero(sp_smeared);
+
+                // No deicated function for a spin matrix esists. Loop over rows in the spinmatrix
+                // which are spinors and use that method.
+                for (int beta=0;beta<4;beta++){
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],norm_factor, sp_OG.c[beta]);
+                }
+
+                // For the multiplication with a colour matrix we need to loop over every entry in the
+                // spin matrix, i.e. two loops are needed. Assignment and addition can be done in the 
+                // outermost loop since deicated functions for spinors exist.
+                for (int beta=0; beta<4; beta++){
+                  for (int alpha=0; alpha<4; alpha++){
+                   _suNf_multiply(sp_tmp.c[beta].c[alpha], *pu_gauge_APE_f(ix,1), sp_right.c[beta].c[alpha]);
+                  }
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],epsilon*norm_factor, sp_tmp.c[beta]);
+                }
+                
+                // Repeat the same approach for the inverse multiplication and smearing directions
+                for (int beta=0; beta<4; beta++){
+                  for (int alpha=0; alpha<4; alpha++){
+                   _suNf_inverse_multiply(sp_tmp.c[beta].c[alpha], *pu_gauge_APE_f(ix_left,1), sp_left.c[beta].c[alpha]);
+                  }
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],epsilon*norm_factor, sp_tmp.c[beta]);
+                }
+
+                for (int beta=0; beta<4; beta++){
+                  for (int alpha=0; alpha<4; alpha++){
+                   _suNf_multiply(sp_tmp.c[beta].c[alpha], *pu_gauge_APE_f(ix,2), sp_front.c[beta].c[alpha]);
+                  }
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],epsilon*norm_factor, sp_tmp.c[beta]);
+                }
+                
+                for (int beta=0; beta<4; beta++){
+                  for (int alpha=0; alpha<4; alpha++){
+                   _suNf_inverse_multiply(sp_tmp.c[beta].c[alpha], *pu_gauge_APE_f(ix_back,2), sp_back.c[beta].c[alpha]);
+                  }
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],epsilon*norm_factor, sp_tmp.c[beta]);
+                }                
+
+                for (int beta=0; beta<4; beta++){
+                  for (int alpha=0; alpha<4; alpha++){
+                   _suNf_multiply(sp_tmp.c[beta].c[alpha], *pu_gauge_APE_f(ix,3), sp_up.c[beta].c[alpha]);
+                  }
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],epsilon*norm_factor, sp_tmp.c[beta]);
+                }
+                
+                for (int beta=0; beta<4; beta++){
+                  for (int alpha=0; alpha<4; alpha++){
+                   _suNf_inverse_multiply(sp_tmp.c[beta].c[alpha], *pu_gauge_APE_f(ix_down,3), sp_down.c[beta].c[alpha]);
+                  }
+                  _spinor_mul_add_assign_f(sp_smeared.c[beta],epsilon*norm_factor, sp_tmp.c[beta]);
+                }           
+                
+                for (int b=0;b<NF;b++){
+                    for (int alpha=0;alpha<4;alpha++) for (int beta=0;beta<4;beta++){
+                        _FIELD_AT(&smeared_psi[beta*nm+i],ix)->c[alpha].c[b] = sp_smeared.c[beta].c[alpha].c[b];
+                        _FIELD_AT(&smeared_psi[beta*nm+i],ix)->c[alpha].c[b] = sp_smeared.c[beta].c[alpha].c[b];
+                    }
+                }
+               
+            }
+        }
+    }
+    for (int i=0;i<4*nm;i++){
+        spinor_field_copy_f(psi + i, smeared_psi  + i);
+    }
+    for (int i=0;i<4*nm;i++){
+        start_sf_sendrecv(psi + i);
+        complete_sf_sendrecv(psi + i);
+    }
+    free_spinor_field_f(smeared_psi);
+}
diff --git a/LibHR/Observables/meson_measurements.c b/LibHR/Observables/meson_measurements.c
index 52fe5e328..07959f7a4 100644
--- a/LibHR/Observables/meson_measurements.c
+++ b/LibHR/Observables/meson_measurements.c
@@ -761,3 +761,249 @@ void measure_formfactor_fixed(int ti, int tf, int dt, int nm, double *m, int n_m
   free_gfield_f(u_gauge_old);
   free_propagator_eo();
 }
+
+
+
+/* Updates for Wuppertal smearing by HH in 2020*/
+
+void create_smear_source_and_measure(spinor_field* source, spinor_field* prop, int t, int x, int y, int z, int nm, double* m, int n_mom, double epsilon_source, int Nsmear_source, int APE_N,int conf_num, char* label) {
+  /**
+   * @brief Creates a Wuppetal smeared source and measures the meson correlators with point sink.
+   *        The APE smearing on the gauge links are applied if `APE_N` is not 0.
+   *        The smearing function, `smearing_function`, is imbedded in `create_smeared_source_with_APE` and `create_smeared_source`.
+   *        See more details in that function.
+   * 
+   * @param source The source field, $q(x)$.
+   * @param prop The fermion propagator for measuring the correlators.
+   * @param epsilon_source The smearing parameter, $\epsilon$, at the source.
+   * @param Nsmear_source The number of iterations for the Wuppertal smearing at the source.
+   * @param APE_N The number of iteration for the APE smearing, determining using `create_smeared_source_with_APE` or `create_smeared_source`
+   * @param t The time slice where the source is located.
+   * @param x The location of the source.
+   * @param y The location of the source.
+   * @param z The location of the source.
+   * @param nm Number of the measuring fermion mass
+   * @param m fermion masses
+   */
+  for (int  k = 0;k<NF;++k){
+    if (APE_N != 0){
+      create_smeared_source_with_APE(source, t, x, y, z, k, epsilon_source, Nsmear_source);
+    } else{
+      create_smeared_source(source, t, x, y, z, k, epsilon_source, Nsmear_source);
+    }
+    calc_propagator(prop + 4*k, source , 4);//4 for spin components
+    if (n_mom>1){
+      measure_point_mesons_momenta(meson_correlators,prop+4*k, source, nm, t, n_mom);
+    } else{
+      measure_mesons(meson_correlators,prop+4*k, source, nm, t);
+    }
+  }
+      
+  sprintf(label, "source_N%d_sink_N%d",Nsmear_source, 0);
+  print_mesons(meson_correlators,1.,conf_num,nm,m,GLB_T,n_mom,label);
+      
+}
+
+
+void smear_sink_and_measure(spinor_field* prop, spinor_field* source, int nm, double* m, int t, int n_mom, int Nsmear_source, double epsilon_sink, int Nsmear_sink, int N_diff, int APE_N, int conf_num, char* label) {
+    /**
+     * @brief Performs the Wuppertal smearing at the sink and measures the meson correlators, with the smearing steps from 0 to `Nsmear_sink` with a step of `N_diff`.
+     *        The numbers of the smearing iterations at the sink are from 0 to `Nsmear_sink` with a step of `N_diff`.
+     *        The smearing function is imbedded in `smeared_propagator` and `smeared_propagator_with_APE`.
+     * 
+     * @note The APE smearing on the gauge links are applied if `APE_N` is not 0.
+     */
+    for (int N = 0; N<Nsmear_sink; N += N_diff){
+      for(int steps=0;steps<N_diff;steps++){
+            lprintf("SMEAR",0,"sink smearing with APE steps: %d\n ", steps+N+1);
+            if (APE_N != 0){
+              smeared_propagator_with_APE(prop, nm, epsilon_sink);
+            } else {
+              smeared_propagator(prop, nm, epsilon_sink);
+            }
+              
+          }
+          
+          for (int c=0;c<NF;++c){
+            if (n_mom>1){
+              measure_point_mesons_momenta(meson_correlators,prop+4*c, source, nm, t, n_mom);
+            } else {
+              measure_mesons(meson_correlators,prop+4*c, source, nm, t);
+            }
+          }
+          sprintf(label, "source_N%d_sink_N%d",Nsmear_source, N+N_diff);
+          print_mesons(meson_correlators,1.,conf_num,nm,m,GLB_T,n_mom,label);
+      }
+}
+
+
+void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source_max, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff){
+    /**
+     * @brief measures the meson correlators with smeared source and sink.
+     * It measures the meson correlators for all combinations of the source and sink smearing steps,
+     * [0, Nsmear_source_max] and [0, Nsmear_sink], respectively, with a separation of `N_diff`.
+     * 
+     * @note The APE smearing on the gauge links are applied if `APE_N` is not 0.
+     */
+    char label[20];    
+    if (APE_N != 0){
+        APE_smearing(APE_epsilon, APE_N);
+        represent_gauge_field_APE();
+    }
+
+    for (int Nsmear_source=0; Nsmear_source <= Nsmear_source_max; Nsmear_source += N_diff){
+
+      spinor_field* source = alloc_spinor_field_f(4,&glattice);
+      spinor_field* prop = alloc_spinor_field_f(4*nm*NF, &glattice);
+      init_propagator_eo(nm,m,precision);
+
+      create_smear_source_and_measure(source, prop, t, x, y, z, nm, m, n_mom, epsilon_source, Nsmear_source, APE_N, conf_num, label);
+
+      smear_sink_and_measure(prop, source, nm, m, t, n_mom, Nsmear_source, epsilon_sink, Nsmear_sink, N_diff, APE_N, conf_num, label);
+      
+      free_propagator_eo();
+      free_spinor_field_f(source);
+      free_spinor_field_f(prop);
+    }
+}
+
+void measure_smearing_ss(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff){
+  /**
+   * @brief measures the meson correlators with smeared source and sink.
+   * It measures the meson correlators for a specifc source smearing steps, Nsmear_source,
+   * while the sink smearing steps are in the range of [0, Nsmear_sink] with a separation of `N_diff` .
+   * 
+   * @note The APE smearing on the gauge links are applied if `APE_N` is not 0.
+   */
+
+  char label[20];    
+  if (APE_N != 0){
+    APE_smearing(APE_epsilon, APE_N);
+    represent_gauge_field_APE();
+  }
+
+  spinor_field* source = alloc_spinor_field_f(4,&glattice);
+  spinor_field* prop = alloc_spinor_field_f(4*nm*NF, &glattice);
+  init_propagator_eo(nm,m,precision);
+
+  create_smear_source_and_measure(source, prop, t, x, y, z, nm, m, n_mom, epsilon_source, Nsmear_source, APE_N,  conf_num, label);
+
+  smear_sink_and_measure(prop, source, nm, m, t, n_mom, Nsmear_source, epsilon_sink, Nsmear_sink, N_diff, APE_N, conf_num, label);
+        
+  free_propagator_eo();
+  free_spinor_field_f(source);
+  free_spinor_field_f(prop);
+
+}
+
+
+// So far, only implement the measurement for mesons at rest
+// The sources are not localised on any lattice site but spread out
+/* Smearing for disconnected pieces for singlets (FZ Jan. 2024)*/
+void measure_spectrum_discon_semwall_smeared_single_inversion(int nm, double* m, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double APE_epsilon, int APE_N, int N_diff){
+    
+    char label[32];
+    int k,beta;
+    spinor_field* source = alloc_spinor_field_f(4,&glattice);
+    spinor_field* prop = alloc_spinor_field_f(4*nm, &glattice);
+    
+    init_propagator_eo(nm,m,precision);
+    if (APE_N != 0){
+        APE_smearing(APE_epsilon, APE_N);
+        represent_gauge_field_APE();
+    }
+
+    lprintf("SMEAR",0,"source smearing with steps %d up to  %d\n ",N_diff,Nsmear_source);
+
+    for (k=0;k<nhits;++k){
+      // First create a new source
+      create_noise_source_equal_eo(source);
+      
+      // Then calculate the propagator without smearing
+      // (Do we need to swtich the type of the lattice?)
+      for(beta=0;beta<4;beta++) source[beta].type = &glat_even;
+      calc_propagator(prop,source,4); // 4 for spin dilution
+      for(beta=0;beta<4;beta++) source[beta].type = &glattice;
+
+      lprintf("MESON_MEASUREMENTS",0,"Stochastic source #%d\n",k);
+
+      for (int N = 0; N<=Nsmear_source; N += N_diff){
+        // Apply smearing now to the source
+        // Then apply smearing to the propagator
+
+        // For the first measurement we skip smearing
+        if (N != 0){
+          if (APE_N != 0){
+              for (int n=0;n<N_diff;n++){
+                lprintf("SMEAR",0,"smearing step %d\n", n+1);
+                smearing_function_volume_with_APE(source, epsilon_source);
+                smeared_propagator_volume_with_APE(prop, nm, epsilon_source);
+              }
+          }
+          else{
+              for (int n=0;n<N_diff;n++){
+                lprintf("SMEAR",0,"smearing step %d\n", n+1);
+                smearing_function_volume(source, epsilon_source);
+                smeared_propagator_volume(prop, nm, epsilon_source);
+              }
+          }
+        }
+
+        // The time slice is set to tau=0 since this is only needed for connected pieces
+        measure_mesons(discon_correlators,prop,source,nm,0);
+              
+        sprintf(label,"src %d DISCON_SEMWALL smear_N%d",k,N);
+        print_mesons(discon_correlators,1.,conf_num,nm,m,GLB_T,1,label);
+      }
+    } 
+       
+    free_propagator_eo();
+    free_spinor_field_f(source);
+    free_spinor_field_f(prop);
+}
+
+
+// So far, only implement the measurement for mesons at rest
+// The sources are not localised on any lattice site but spread out
+/* Smearing for disconnected pieces for singlets (FZ Jan. 2024)*/
+void measure_spectrum_discon_semwall_smeared(int nm, double* m, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double APE_epsilon, int APE_N, int N_diff){
+    
+    char label[32];
+    int k,beta;
+    spinor_field* source = alloc_spinor_field_f(4,&glattice);
+    spinor_field* prop = alloc_spinor_field_f(4*nm, &glattice);
+    
+    init_propagator_eo(nm,m,precision);
+    if (APE_N != 0){
+        APE_smearing(APE_epsilon, APE_N);
+        represent_gauge_field_APE();
+    }
+
+    lprintf("SMEAR",0,"source smearing with steps %d up to  %d\n ",N_diff,Nsmear_source);
+
+    for (k=0;k<nhits;++k){
+      for (int N = 0; N<=Nsmear_source; N += N_diff){
+    
+        if (APE_N != 0){
+            create_noise_source_equal_eo_smeared_with_APE(source, epsilon_source, N);
+        }
+        else{
+            create_noise_source_equal_eo_smeared(source, epsilon_source, N);
+        }
+        // Do we need to swtich the type of the lattice?
+        for(beta=0;beta<4;beta++) source[beta].type = &glat_even;
+        calc_propagator(prop,source,4); // 4 for spin dilution
+        for(beta=0;beta<4;beta++) source[beta].type = &glattice;
+
+        // The time slice is set to tau=0 since this is only needed for connected pieces
+        measure_mesons(discon_correlators,prop,source,nm,0);
+              
+        sprintf(label,"src %d DISCON_SEMWALL smear_N%d",k,N);
+        print_mesons(discon_correlators,1.,conf_num,nm,m,GLB_T,1,label);
+      }
+    } 
+       
+    free_propagator_eo();
+    free_spinor_field_f(source);
+    free_spinor_field_f(prop);
+}
\ No newline at end of file
diff --git a/LibHR/Observables/polyakov.c b/LibHR/Observables/polyakov.c
index 7c72c93e1..63b488f8a 100644
--- a/LibHR/Observables/polyakov.c
+++ b/LibHR/Observables/polyakov.c
@@ -457,5 +457,258 @@ void polyakov()
 }
 #endif
 
+
+void polyakov_APE() {
+  /**
+   * @brief Computes the value of polyakov loop for the APE-smeared gauge links.
+   *        Modified from `polyakov()` by changing `pu_gauge(...)` to `pu_gauge_APE(...)`
+   * 
+   */
+  int x[4], i4d, i3d, size3d;
+  int loc[4]={T,X,Y,Z};
+  suNg *p;
+  suNg *lp;
+  suNg *bp;
+  suNg tmp;
+  double complex poly;
+  double adjpoly;
+  double dtmp;
+#ifdef WITH_MPI
+  int np[4]={NP_T,NP_X,NP_Y,NP_Z};
+  int mpiret;
+#endif /* WITH_MPI */
+ 
+  for(int mu=0; mu<4; mu++){
+
+    size3d=T*X*Y*Z/loc[mu];
+    p=malloc(sizeof(suNg)*size3d);
+    lp=malloc(sizeof(suNg)*size3d);
+    bp=malloc(sizeof(suNg)*size3d);
+    
+    
+    /* local wilson lines */
+    
+    i3d=0;
+    for(x[(mu+1)%4]=0; x[(mu+1)%4]<loc[(mu+1)%4]; x[(mu+1)%4]++)
+    for(x[(mu+2)%4]=0; x[(mu+2)%4]<loc[(mu+2)%4]; x[(mu+2)%4]++)
+    for(x[(mu+3)%4]=0; x[(mu+3)%4]<loc[(mu+3)%4]; x[(mu+3)%4]++) {
+      _suNg_unit(lp[i3d]);
+      for(x[mu]=0; x[mu]<loc[mu]; x[mu]++) {
+        i4d=ipt(x[0],x[1],x[2],x[3]);
+        _suNg_times_suNg(tmp,lp[i3d],*pu_gauge_APE(i4d,mu));
+        lp[i3d] = tmp;
+      }
+      i3d++;
+    }
+  
+  
+    
+    /* global wilson lines */
+    
+  #ifdef WITH_MPI
+    if(COORD[mu]!=0) {
+      MPI_Status status;
+      mpiret=MPI_Recv(bp, /* buffer */
+          size3d*sizeof(suNg)/sizeof(double), /* lenght in units of doubles */
+          MPI_DOUBLE, /* basic datatype */
+          proc_dn(CID,mu), /* cid of origin */
+          COORD[mu]-1, /* tag of communication */
+          cart_comm, /* use the cartesian communicator */
+          &status);
+  #ifndef NDEBUG
+      if (mpiret != MPI_SUCCESS) {
+        char mesg[MPI_MAX_ERROR_STRING];
+        int mesglen;
+        MPI_Error_string(mpiret,mesg,&mesglen);
+        lprintf("MPI",0,"ERROR: %s\n",mesg);
+        if (status.MPI_ERROR != MPI_SUCCESS) {
+          MPI_Error_string(status.MPI_ERROR,mesg,&mesglen);
+          lprintf("MPI",0,"Req [%d] Source [%d] Tag [%] ERROR: %s\n",
+              0,
+              status.MPI_SOURCE,
+              status.MPI_TAG,
+              mesg);
+        }
+        error(1,1,"WL_Hamiltonian_gauge [wilsonloops.c]","Cannot receive buf_gtf[1]");
+      }
+  #endif /* NDEBUG */
+    }
+  #endif /* WITH_MPI */
+  
+  
+    if(COORD[mu]==0) {
+      i3d=0;
+      for(x[(mu+1)%4]=0; x[(mu+1)%4]<loc[(mu+1)%4]; x[(mu+1)%4]++)
+      for(x[(mu+2)%4]=0; x[(mu+2)%4]<loc[(mu+2)%4]; x[(mu+2)%4]++)
+      for(x[(mu+3)%4]=0; x[(mu+3)%4]<loc[(mu+3)%4]; x[(mu+3)%4]++) {
+        p[i3d]=lp[i3d];
+        i3d++;
+      }
+    } else {
+      i3d=0;
+      for(x[(mu+1)%4]=0; x[(mu+1)%4]<loc[(mu+1)%4]; x[(mu+1)%4]++)
+      for(x[(mu+2)%4]=0; x[(mu+2)%4]<loc[(mu+2)%4]; x[(mu+2)%4]++)
+      for(x[(mu+3)%4]=0; x[(mu+3)%4]<loc[(mu+3)%4]; x[(mu+3)%4]++) {
+        _suNg_times_suNg(p[i3d], bp[i3d], lp[i3d]);
+        i3d++;
+      }
+    }
+    
+    
+  #ifdef WITH_MPI
+    if(COORD[mu]!=np[mu]-1) {
+      mpiret=MPI_Send(p, /* buffer */
+          size3d*sizeof(suNg)/sizeof(double), /* lenght in units of doubles */
+          MPI_DOUBLE, /* basic datatype */
+          proc_up(CID,mu), /* cid of destination */
+          COORD[mu], /* tag of communication */
+          cart_comm /* use the cartesian communicator */
+          );
+  #ifndef NDEBUG
+      if (mpiret != MPI_SUCCESS) {
+        char mesg[MPI_MAX_ERROR_STRING];
+        int mesglen;
+        MPI_Error_string(mpiret,mesg,&mesglen);
+        lprintf("MPI",0,"ERROR: %s\n",mesg);
+        error(1,1,"WL_Hamiltonian_gauge [wilsonloops.c]","Cannot send buf_gtf[0]");
+      }
+  #endif /* NDEBUG */
+    }
+  #endif /* WITH_MPI */
+  
+  
+  
+    /* broadcast wilson lines */
+  
+  #ifdef WITH_MPI
+    int mpiret;
+    if(COORD[mu]==np[mu]-1) {
+      MPI_Request comm_req[np[mu]-1];
+      int destCID=CID;
+      for(int t=0; t<np[mu]-1; t++) {
+        destCID=proc_up(destCID,mu);
+        mpiret=MPI_Isend(p, /* buffer */
+            size3d*sizeof(suNg)/sizeof(double), /* lenght in units of doubles */
+            MPI_DOUBLE, /* basic datatype */
+            destCID, /* cid of destination */
+            t, /* tag of communication */
+            cart_comm, /* use the cartesian communicator */
+            &(comm_req[t]) /* handle to communication request */
+            );
+  #ifndef NDEBUG
+        if (mpiret != MPI_SUCCESS) {
+          char mesg[MPI_MAX_ERROR_STRING];
+          int mesglen;
+          MPI_Error_string(mpiret,mesg,&mesglen);
+          lprintf("MPI",0,"ERROR: %s\n",mesg);
+          error(1,1,"WL_broadcast_polyakov [wilsonloops.c]","Cannot start send polyakov");
+        }
+  #endif /* NDEBUG */
+      }
+      
+      MPI_Status status[np[mu]-1];
+      mpiret=MPI_Waitall(np[mu]-1, comm_req, status);
+  #ifndef NDEBUG
+      if (mpiret != MPI_SUCCESS) {
+        char mesg[MPI_MAX_ERROR_STRING];
+        int mesglen, k;
+        MPI_Error_string(mpiret,mesg,&mesglen);
+        lprintf("MPI",0,"ERROR: %s\n",mesg);
+        for (k=0; k<np[mu]-1; ++k) {
+          if (status[k].MPI_ERROR != MPI_SUCCESS) {
+            MPI_Error_string(status[k].MPI_ERROR,mesg,&mesglen);
+            lprintf("MPI",0,"Req [%d] Source [%d] Tag [%] ERROR: %s\n",
+                k,
+                status[k].MPI_SOURCE,
+                status[k].MPI_TAG,
+                mesg);
+          }
+        }
+        error(1,1,"WL_broadcast_polyakov [wilsonloops.c]","Cannot complete communications");
+      }
+  #endif /* NDEBUG */
+  
+    } else {
+  
+      int sCOORD[4], sCID;
+      sCOORD[0]=COORD[0];sCOORD[1]=COORD[1];sCOORD[2]=COORD[2];sCOORD[3]=COORD[3];
+      sCOORD[mu]=np[mu]-1;
+      mpiret=MPI_Cart_rank(cart_comm, sCOORD, &sCID);
+  #ifndef NDEBUG
+      if (mpiret != MPI_SUCCESS) {
+        char mesg[MPI_MAX_ERROR_STRING];
+        int mesglen;
+        MPI_Error_string(mpiret,mesg,&mesglen);
+        lprintf("MPI",0,"ERROR: %s\n",mesg);
+        error(1,1,"WL_broadcast_polyakov [wilsonloops.c]","Cannot retrieve source CID");
+      }
+  #endif /* NDEBUG */
+      MPI_Status status;
+      mpiret=MPI_Recv(p, /* buffer */
+          size3d*sizeof(suNg)/sizeof(double), /* lenght in units of doubles */
+          MPI_DOUBLE, /* basic datatype */
+          sCID, /* cid of destination */
+          COORD[mu], /* tag of communication */
+          cart_comm, /* use the cartesian communicator */
+          &status);
+  #ifndef NDEBUG
+      if (mpiret != MPI_SUCCESS) {
+        char mesg[MPI_MAX_ERROR_STRING];
+        int mesglen;
+        MPI_Error_string(mpiret,mesg,&mesglen);
+        lprintf("MPI",0,"ERROR: %s\n",mesg);
+        if (status.MPI_ERROR != MPI_SUCCESS) {
+          MPI_Error_string(status.MPI_ERROR,mesg,&mesglen);
+          lprintf("MPI",0,"Req [%d] Source [%d] Tag [%] ERROR: %s\n",
+              0,
+              status.MPI_SOURCE,
+              status.MPI_TAG,
+              mesg);
+        }
+        error(1,1,"WL_broadcast_polyakov [wilsonloops.c]","Cannot receive polyakov");
+      }
+  #endif /* NDEBUG */
+    }
+  #endif /* WITH_MPI */
+  
+  
+    
+    /* trace and average */
+    _complex_0(poly);
+    adjpoly=0.;
+  
+    i3d=0;
+    for(x[(mu+1)%4]=0; x[(mu+1)%4]<loc[(mu+1)%4]; x[(mu+1)%4]++)
+    for(x[(mu+2)%4]=0; x[(mu+2)%4]<loc[(mu+2)%4]; x[(mu+2)%4]++)
+    for(x[(mu+3)%4]=0; x[(mu+3)%4]<loc[(mu+3)%4]; x[(mu+3)%4]++) {
+      _suNg_trace_re(dtmp,p[i3d]);
+/*      lprintf("LOC_POLYAKOV",0,"%d %d %d %d %d %1.8e\n",mu,x[(mu+1)%4],x[(mu+2)%4],x[(mu+3)%4],i3d,dtmp); */
+      poly += dtmp;
+      adjpoly +=dtmp*dtmp;
+#ifndef GAUGE_SON
+      _suNg_trace_im(dtmp,p[i3d]);
+      poly += I * dtmp;
+      adjpoly +=dtmp*dtmp - 1;
+#endif
+      i3d++;
+    }
+  
+    global_sum((double*)&poly,2);
+    global_sum((double*)&adjpoly,1);
+  
+    
+    poly /= NG*(GLB_VOLUME/loc[mu]);
+    adjpoly /= NG*(GLB_VOLUME/loc[mu]);
+    
+    lprintf("FUND_POLYAKOV",0,"Polyakov direction %d = %1.8e %1.8e\n",mu,creal(poly), cimag(poly));
+    lprintf("ADJ_POLYAKOV",0,"Polyakov direction %d = %1.8e\n",mu,adjpoly);
+    
+    free(p);
+    free(lp);
+    free(bp);
+    
+  }
+}
+
 #undef MPIRET
-#undef _MIN
\ No newline at end of file
+#undef _MIN
diff --git a/LibHR/Observables/sources.c b/LibHR/Observables/sources.c
index 7bdef15a8..69acdc08c 100755
--- a/LibHR/Observables/sources.c
+++ b/LibHR/Observables/sources.c
@@ -25,6 +25,7 @@
 #include "gamma_spinor.h"
 #include "spin_matrix.h"
 #include "propagator.h"
+#include "representation.h"
 
 #define PI 3.141592653589793238462643383279502884197
 
@@ -88,6 +89,11 @@ static int random_tau()
 						source[spin](t,x) = \sum_{x%p == 0} \delta_{s spin} 1_color
 
                 z2_volume_source:                  source[spin](t,x) = Z(2) x Z(2) (no dilution)
+    smeared_source:
+              Wuppertal smeared source, see Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf.
+              The involving gauge field is the original one.
+    smeared_source_with_APE:
+              Wuppertal smearing with APE smeared gauge field.
 \***************************************************************************/
 void create_point_source(spinor_field *source, int tau, int color)
 {
@@ -711,3 +717,406 @@ void zero_even_or_odd_site_spinorfield(spinor_field *source, int nspinor, int eo
           }
         }
 }
+
+
+/* Updates for Wuppertal smearing by HH in 2020*/
+
+void smearing_function(spinor_field *source, int tau, double epsilon){
+  /**
+   * @brief Applies the Wuppertal smearing to the source field.
+   * The smearing is done in the spatial directions only, basing on Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf:
+   * $\Phi q(x) = 1/(1+2d\epsilon)[ q(x) + \epsilon\sum_{\mu=1}^{3} U_{\mu}(x)q(x+\mu) ]$
+   * $\Phi$ is the smearing function, $q(x)$ is the source field, $U_{\mu}(x)$ is the link variable in the $\mu$ direction.
+   * $d = 3$ for spatial directions.
+   * 
+   * @param source -> the source field, $q(x)$.
+   * @param epsilon -> the smearing parameter, $\epsilon$.
+   * @param tau is the time slice where the source is located.
+   */
+    
+    int ix, x, y, z, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down;
+    double norm_factor = 1./(1.+6.*epsilon);
+    suNf_spinor spinor_OG, spinor_smeared, spinor_tmp, spinor_right, spinor_left, spinor_front, spinor_back, spinor_up, spinor_down;
+    _spinor_zero_f(spinor_OG);_spinor_zero_f(spinor_tmp);_spinor_zero_f(spinor_right);_spinor_zero_f(spinor_left);
+    _spinor_zero_f(spinor_front);_spinor_zero_f(spinor_back);_spinor_zero_f(spinor_up);_spinor_zero_f(spinor_down);
+        
+    spinor_field* smeared_source = alloc_spinor_field_f(4,&glattice);
+    for (int beta=0;beta<4;++beta){
+        spinor_field_zero_f(&smeared_source[beta]);
+    }
+    
+    if(COORD[0]==tau/T){ // apply the smearing only on the time slice where the source is located.
+        for (x=0; x<X; x++) for (y=0; y<Y; y++) for (z=0; z<Z; z++){
+            
+            // assign the positions: a site (tau, x, y, z) and its neighbors in the spatial directions.
+            ix = ipt(tau - zerocoord[0], x,y,z);
+            ix_right = iup(ix,1);
+            ix_left  = idn(ix,1);
+            ix_front = iup(ix,2);
+            ix_back  = idn(ix,2);
+            ix_up    = iup(ix,3);
+            ix_down  = idn(ix,3);
+    
+    
+    
+            for (int spin=0;spin<4;spin++){
+                spinor_OG    = *_FIELD_AT(&source[spin], ix); // OG: original spinor field at (x,y,z).
+                spinor_right = *_FIELD_AT(&source[spin], ix_right); 
+                spinor_left  = *_FIELD_AT(&source[spin], ix_left);
+                spinor_front = *_FIELD_AT(&source[spin], ix_front);
+                spinor_back  = *_FIELD_AT(&source[spin], ix_back);
+                spinor_up    = *_FIELD_AT(&source[spin], ix_up);
+                spinor_down  = *_FIELD_AT(&source[spin], ix_down);
+                        
+        
+                _spinor_zero_f(spinor_smeared);
+
+                // mutiply the spinor fields with link variable in the spatial directions and add to the smearing field.
+                _suNf_multiply(spinor_tmp.c[spin], *pu_gauge_f(ix,1), spinor_right.c[spin]);
+                _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+        
+                _suNf_inverse_multiply(spinor_tmp.c[spin], *pu_gauge_f(ix_left,1), spinor_left.c[spin]);
+                _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+                    
+                _suNf_multiply(spinor_tmp.c[spin], *pu_gauge_f(ix,2), spinor_front.c[spin]);
+                _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+        
+                _suNf_inverse_multiply(spinor_tmp.c[spin], *pu_gauge_f(ix_back,2), spinor_back.c[spin]);
+                _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+                    
+                _suNf_multiply(spinor_tmp.c[spin], *pu_gauge_f(ix,3), spinor_up.c[spin]);
+                _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+        
+                _suNf_inverse_multiply(spinor_tmp.c[spin], *pu_gauge_f(ix_down,3), spinor_down.c[spin]);
+                _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+                    
+                _vector_mul_add_assign_f(spinor_smeared.c[spin], norm_factor, spinor_OG.c[spin]);
+        
+        _FIELD_AT(&smeared_source[spin], ix)->c[spin]  = spinor_smeared.c[spin];
+               
+        }
+    }
+    
+    // beta is the spinor index.
+    for (int beta=0;beta<4;++beta){
+        spinor_field_copy_f(source + beta, smeared_source +beta);
+    }
+    }
+    
+    for (int beta=0;beta<4;++beta){
+        start_sf_sendrecv(source + beta);
+        complete_sf_sendrecv(source + beta);
+    }
+    free_spinor_field_f(smeared_source);
+}
+
+void smearing_function_with_APE(spinor_field *source, int tau, double epsilon){
+  /**
+   * @brief Applies the Wuppertal smearing to the source field with APE smeared gauge field,
+   *        which is similar to the function `smearing_function` but changes the gauge field to the APE smeared one:
+   *        `pu_gauge_f` -> `pu_gauge_APE_f`
+   */
+    
+    int ix, x, y, z, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down;
+    double norm_factor = 1./(1.+6.*epsilon);
+    suNf_spinor spinor_OG, spinor_smeared, spinor_tmp, spinor_right, spinor_left, spinor_front, spinor_back, spinor_up, spinor_down;
+    _spinor_zero_f(spinor_OG);_spinor_zero_f(spinor_tmp);_spinor_zero_f(spinor_right);_spinor_zero_f(spinor_left);
+    _spinor_zero_f(spinor_front);_spinor_zero_f(spinor_back);_spinor_zero_f(spinor_up);_spinor_zero_f(spinor_down);
+    
+    spinor_field* smeared_source = alloc_spinor_field_f(4, &glattice);
+    for (int beta=0;beta<4;++beta){
+        spinor_field_zero_f(&smeared_source[beta]);
+    }
+    
+    represent_gauge_field_APE();
+    if(COORD[0]==tau/T){
+    for (x=0; x<X; x++) for (y=0; y<Y; y++) for (z=0; z<Z; z++) {
+            
+        //lprintf("SMEAR",0,"smearing function at(%d,%d,%d,%d)\n",tau,x,y,z);
+        ix = ipt(tau - zerocoord[0], x, y, z);
+        ix_right = iup(ix,1);
+        ix_left  = idn(ix,1);
+        ix_front = iup(ix,2);
+        ix_back  = idn(ix,2);
+        ix_up    = iup(ix,3);
+        ix_down  = idn(ix,3);
+            
+        for (int spin=0;spin<4;spin++){
+            spinor_OG    = *_FIELD_AT(&source[spin], ix);
+            spinor_right = *_FIELD_AT(&source[spin], ix_right);
+            spinor_left  = *_FIELD_AT(&source[spin], ix_left);
+            spinor_front = *_FIELD_AT(&source[spin], ix_front);
+            spinor_back  = *_FIELD_AT(&source[spin], ix_back);
+            spinor_up    = *_FIELD_AT(&source[spin], ix_up);
+            spinor_down  = *_FIELD_AT(&source[spin], ix_down);
+                
+            _spinor_zero_f(spinor_smeared);
+            
+            _suNf_multiply(spinor_tmp.c[spin], *pu_gauge_APE_f(ix,1), spinor_right.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+            _suNf_inverse_multiply(spinor_tmp.c[spin], *pu_gauge_APE_f(ix_left,1), spinor_left.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+            
+            _suNf_multiply(spinor_tmp.c[spin], *pu_gauge_APE_f(ix,2), spinor_front.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+            _suNf_inverse_multiply(spinor_tmp.c[spin], *pu_gauge_APE_f(ix_back,2), spinor_back.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+            
+            _suNf_multiply(spinor_tmp.c[spin], *pu_gauge_APE_f(ix,3), spinor_up.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+            _suNf_inverse_multiply(spinor_tmp.c[spin], *pu_gauge_APE_f(ix_down,3), spinor_down.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+            
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], norm_factor, spinor_OG.c[spin]);
+                
+            *_FIELD_AT(&smeared_source[spin], ix)  = spinor_smeared;
+        }
+    }
+    for (int beta=0;beta<4;++beta){
+        spinor_field_copy_f(source + beta, smeared_source +beta);
+    }
+    }
+    
+    for (int beta=0;beta<4;++beta){
+        start_sf_sendrecv(source + beta);
+        complete_sf_sendrecv(source + beta);
+    }
+    free_spinor_field_f(smeared_source);
+}
+
+void create_smeared_source(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear){
+  /**
+   * @brief Creates a smeared source field at a given position (t,x,y,z) with a given color by applying
+   *        the Wuppertal smearing function on a point source based on Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf.
+   *        
+   * @param t The time slice where the source is located.
+   * @param x The location of the source.
+   * @param y The location of the source.
+   * @param z The location of the source.
+   * @param Nsmear is the iteration number at the source
+   * @param epsilon is the step size, $\epsilon$ in the reference.
+   * 
+   * @note The gauge field is the original one.
+   */
+    
+    int beta;
+    
+    for (beta=0;beta<4;++beta){
+        spinor_field_zero_f(&source[beta]);
+    }
+    
+    create_point_source_loc(source, t, x, y, z, color);
+    
+    lprintf("SMEAR",0,"Smeared Source at (%d,%d,%d,%d) with APE smearing \n",t, x,y,z);
+    lprintf("SMEAR",0,"source smearing epsilon = %f iterations: \n", epsilon);
+    
+    for (int n=0;n<Nsmear;n++){
+        lprintf("SMEAR",0,"%d...", n+1);
+        smearing_function(source, t, epsilon);
+    }
+}
+
+void create_smeared_source_with_APE(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear){
+  /**
+   * @brief Creates a smeared source field at a given position (t,x,y,z) with a given color by applying
+   *        the Wuppertal smearing function on a point source based on Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf.
+   *        
+   * @param t The time slice where the source is located.
+   * @param x The location of the source.
+   * @param y The location of the source.
+   * @param z The location of the source.
+   * @param Nsmear is the iteration number at the source
+   * @param epsilon is the step size, $\epsilon$ in the reference.
+   * 
+   * @note The gauge field is the APE smeared one.
+   */
+    int beta;
+    
+    for (beta=0;beta<4;++beta){
+        spinor_field_zero_f(&source[beta]);
+    }
+    
+    create_point_source_loc(source, t, x, y, z, color);
+    
+    lprintf("SMEAR",0,"Smeared Source at (%d,%d,%d,%d) with APE smearing \n",t, x,y,z);
+    lprintf("SMEAR",0,"source smearing epsilon = %f iterations: \n", epsilon);
+    
+    for (int n=0;n<Nsmear;n++){
+        lprintf("SMEAR",0,"%d...", n+1);
+        smearing_function_with_APE(source, t, epsilon);
+    }
+}
+
+// FZ 2024: Add extra smeared sources for disconnected pieces
+void smearing_function_volume(spinor_field *source, double epsilon){
+    
+    int ix, t, x, y, z, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down;
+    double norm_factor = 1./(1.+6.*epsilon);
+    suNf_spinor spinor_OG, spinor_smeared, spinor_tmp, spinor_right, spinor_left, spinor_front, spinor_back, spinor_up, spinor_down;
+    _spinor_zero_f(spinor_OG);_spinor_zero_f(spinor_tmp);_spinor_zero_f(spinor_right);_spinor_zero_f(spinor_left);
+    _spinor_zero_f(spinor_front);_spinor_zero_f(spinor_back);_spinor_zero_f(spinor_up);_spinor_zero_f(spinor_down);
+        
+    spinor_field* smeared_source = alloc_spinor_field_f(4,&glattice);
+    for (int beta=0;beta<4;++beta){
+        spinor_field_zero_f(&smeared_source[beta]);
+    }
+    
+    for (t=0; t<T; t++) for (x=0; x<X; x++) for (y=0; y<Y; y++) for (z=0; z<Z; z++){
+        
+        ix = ipt(t, x,y,z);
+        ix_right = iup(ix,1);
+        ix_left  = idn(ix,1);
+        ix_front = iup(ix,2);
+        ix_back  = idn(ix,2);
+        ix_up    = iup(ix,3);
+        ix_down  = idn(ix,3);
+
+        for (int spin=0;spin<4;spin++){
+            spinor_OG    = *_FIELD_AT(&source[spin], ix);
+            spinor_right = *_FIELD_AT(&source[spin], ix_right);
+            spinor_left  = *_FIELD_AT(&source[spin], ix_left);
+            spinor_front = *_FIELD_AT(&source[spin], ix_front);
+            spinor_back  = *_FIELD_AT(&source[spin], ix_back);
+            spinor_up    = *_FIELD_AT(&source[spin], ix_up);
+            spinor_down  = *_FIELD_AT(&source[spin], ix_down);
+                    
+    
+            _spinor_zero_f(spinor_smeared);
+    
+            _suNf_multiply(spinor_tmp.c[spin], *pu_gauge_f(ix,1), spinor_right.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+    
+            _suNf_inverse_multiply(spinor_tmp.c[spin], *pu_gauge_f(ix_left,1), spinor_left.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+                
+            _suNf_multiply(spinor_tmp.c[spin], *pu_gauge_f(ix,2), spinor_front.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+    
+            _suNf_inverse_multiply(spinor_tmp.c[spin], *pu_gauge_f(ix_back,2), spinor_back.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+                
+            _suNf_multiply(spinor_tmp.c[spin], *pu_gauge_f(ix,3), spinor_up.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+    
+            _suNf_inverse_multiply(spinor_tmp.c[spin], *pu_gauge_f(ix_down,3), spinor_down.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+                
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], norm_factor, spinor_OG.c[spin]);
+    
+            _FIELD_AT(&smeared_source[spin], ix)->c[spin]  = spinor_smeared.c[spin];       
+        }
+    }
+    
+    for (int beta=0;beta<4;++beta){
+        spinor_field_copy_f(source + beta, smeared_source +beta);
+    }
+    
+    
+    for (int beta=0;beta<4;++beta){
+        start_sf_sendrecv(source + beta);
+        complete_sf_sendrecv(source + beta);
+    }
+    free_spinor_field_f(smeared_source);
+}
+
+void smearing_function_volume_with_APE(spinor_field *source, double epsilon){
+    
+    int ix, t, x, y, z, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down;
+    double norm_factor = 1./(1.+6.*epsilon);
+    suNf_spinor spinor_OG, spinor_smeared, spinor_tmp, spinor_right, spinor_left, spinor_front, spinor_back, spinor_up, spinor_down;
+    _spinor_zero_f(spinor_OG);_spinor_zero_f(spinor_tmp);_spinor_zero_f(spinor_right);_spinor_zero_f(spinor_left);
+    _spinor_zero_f(spinor_front);_spinor_zero_f(spinor_back);_spinor_zero_f(spinor_up);_spinor_zero_f(spinor_down);
+    
+    spinor_field* smeared_source = alloc_spinor_field_f(4, &glattice);
+    for (int beta=0;beta<4;++beta){
+        spinor_field_zero_f(&smeared_source[beta]);
+    }
+    
+    represent_gauge_field_APE();
+
+    for (t=0; t<T; t++) for (x=0; x<X; x++) for (y=0; y<Y; y++) for (z=0; z<Z; z++) {
+            
+        ix = ipt(t, x, y, z);
+        ix_right = iup(ix,1);
+        ix_left  = idn(ix,1);
+        ix_front = iup(ix,2);
+        ix_back  = idn(ix,2);
+        ix_up    = iup(ix,3);
+        ix_down  = idn(ix,3);
+            
+        for (int spin=0;spin<4;spin++){
+            spinor_OG    = *_FIELD_AT(&source[spin], ix);
+            spinor_right = *_FIELD_AT(&source[spin], ix_right);
+            spinor_left  = *_FIELD_AT(&source[spin], ix_left);
+            spinor_front = *_FIELD_AT(&source[spin], ix_front);
+            spinor_back  = *_FIELD_AT(&source[spin], ix_back);
+            spinor_up    = *_FIELD_AT(&source[spin], ix_up);
+            spinor_down  = *_FIELD_AT(&source[spin], ix_down);
+                
+            _spinor_zero_f(spinor_smeared);
+            
+            _suNf_multiply(spinor_tmp.c[spin], *pu_gauge_APE_f(ix,1), spinor_right.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+            _suNf_inverse_multiply(spinor_tmp.c[spin], *pu_gauge_APE_f(ix_left,1), spinor_left.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+            
+            _suNf_multiply(spinor_tmp.c[spin], *pu_gauge_APE_f(ix,2), spinor_front.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+            _suNf_inverse_multiply(spinor_tmp.c[spin], *pu_gauge_APE_f(ix_back,2), spinor_back.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+            
+            _suNf_multiply(spinor_tmp.c[spin], *pu_gauge_APE_f(ix,3), spinor_up.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+            _suNf_inverse_multiply(spinor_tmp.c[spin], *pu_gauge_APE_f(ix_down,3), spinor_down.c[spin]);
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], epsilon*norm_factor, spinor_tmp.c[spin]);
+            
+            _vector_mul_add_assign_f(spinor_smeared.c[spin], norm_factor, spinor_OG.c[spin]);
+                
+            *_FIELD_AT(&smeared_source[spin], ix)  = spinor_smeared;
+        }
+    }
+    for (int beta=0;beta<4;++beta){
+        spinor_field_copy_f(source + beta, smeared_source +beta);
+    }
+    
+    for (int beta=0;beta<4;++beta){
+        start_sf_sendrecv(source + beta);
+        complete_sf_sendrecv(source + beta);
+    }
+    free_spinor_field_f(smeared_source);
+}
+
+void create_noise_source_equal_eo_smeared_with_APE(spinor_field *source, double epsilon, int Nsmear){
+    
+    int beta;
+    for (beta=0;beta<4;++beta){
+        spinor_field_zero_f(&source[beta]);
+    }
+    
+    create_noise_source_equal_eo(source);
+    
+    lprintf("SMEAR",0,"Smeared noisy volume source with APE smearing \n");
+    lprintf("SMEAR",0,"source smearing epsilon = %f iterations: \n", epsilon);
+    
+    for (int n=0;n<Nsmear;n++){
+        lprintf("SMEAR",0,"%d...", n+1);
+        smearing_function_volume_with_APE(source, epsilon);
+    }
+}
+void create_noise_source_equal_eo_smeared(spinor_field *source, double epsilon, int Nsmear){
+    
+    int beta;
+    for (beta=0;beta<4;++beta){
+        spinor_field_zero_f(&source[beta]);
+    }
+
+    create_noise_source_equal_eo(source);
+    
+    lprintf("SMEAR",0,"Smeared noisy volume source \n");
+    lprintf("SMEAR",0,"source smearing epsilon = %f iterations: \n", epsilon);
+    
+    for (int n=0;n<Nsmear;n++){
+        lprintf("SMEAR",0,"%d...", n+1);
+        smearing_function_volume(source, epsilon);
+    }
+}
\ No newline at end of file
diff --git a/LibHR/Update/representation.c b/LibHR/Update/representation.c
index 7b60b750d..2d7d01648 100644
--- a/LibHR/Update/representation.c
+++ b/LibHR/Update/representation.c
@@ -473,3 +473,78 @@ void represent_gauge_field()
   compute_clover_term();
 #endif
 }
+
+void represent_gauge_field_APE()
+{
+#ifdef WITH_SMEARING
+    smear_gauge_field();
+#endif
+
+#ifdef ALLOCATE_REPR_GAUGE_FIELD
+  /* loop on local lattice first */
+  _MASTER_FOR(&glattice,ix) {
+//  for(ip=0;ip<glattice.local_master_pieces;ip++)
+//    for(ix=glattice.master_start[ip];ix<=glattice.master_end[ip];ix++)
+      for (int mu=0;mu<4;mu++) {
+#ifdef WITH_SMEARING
+        suNg *u=_4FIELD_AT(u_gauge_APE_s,ix,mu);
+#else
+        suNg *u=pu_gauge_APE(ix,mu);
+#endif
+        suNf *Ru=pu_gauge_APE_f(ix,mu);
+        #ifdef UNROLL_GROUP_REPRESENT
+          _group_represent(*Ru,*u);
+        #else
+          _group_represent2(Ru,u);
+        #endif
+      }
+  }
+  
+  /* wait gauge field transfer */
+  complete_gf_sendrecv(u_gauge_APE);
+
+  /* loop on the rest of master sites */
+_OMP_PRAGMA ( _omp_parallel )
+  for(int ip=glattice.local_master_pieces;ip<glattice.total_gauge_master_pieces;ip++) {
+_OMP_PRAGMA ( _omp_for )
+    for(int ix=glattice.master_start[ip];ix<=glattice.master_end[ip];ix++)
+      for (int mu=0;mu<4;mu++) {
+#ifdef WITH_SMEARING
+        suNg *u=_4FIELD_AT(u_gauge_APE_s,ix,mu);
+#else
+        suNg *u=pu_gauge_APE(ix,mu);
+#endif
+        suNf *Ru=pu_gauge_APE_f(ix,mu);
+        #ifdef UNROLL_GROUP_REPRESENT
+          _group_represent(*Ru,*u);
+        #else
+          _group_represent2(Ru,u);
+        #endif
+      }
+  }
+  
+  apply_BCs_on_represented_gauge_field_APE();
+#else
+  static int first_time=1;
+  /* wait gauge field transfer */
+  complete_gf_sendrecv(u_gauge_APE);
+
+  if(first_time) {
+    first_time=0;
+#ifdef WITH_SMEARING
+      u_gauge_APE_f=(suNf_field *)((void*)u_gauge_APE_s);
+#else
+    u_gauge_APE_f=(suNf_field *)((void*)u_gauge_APE);
+#endif
+    //    apply_BCs_on_represented_gauge_field(); //Already applied when configuration read or initialized
+  }
+  complete_gf_sendrecv(u_gauge_APE);
+  u_gauge_APE_f=(suNf_field *)((void*)u_gauge_APE);
+
+#endif
+  assign_ud2u_APE_f();
+
+#ifdef WITH_CLOVER
+  compute_clover_term();
+#endif
+}
\ No newline at end of file
diff --git a/LibHR/Utils/APE_smearing.c b/LibHR/Utils/APE_smearing.c
new file mode 100644
index 000000000..35fbf932f
--- /dev/null
+++ b/LibHR/Utils/APE_smearing.c
@@ -0,0 +1,135 @@
+//
+//  APE_smearing.c
+//
+//  Code added by HH on 2020.11.2 for APE smearing,
+//  based on the functions written in Fortran provided by Davide Vadacchino.
+//  
+//
+
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+#include "global.h"
+#include "utils.h"
+#include "suN.h"
+#include "memory.h"
+#include "communications.h"
+#include "logger.h"
+#include "complex.h"
+#include "random.h"
+#include "representation.h"
+#include "observables.h"
+#include "update.h"
+
+void APE_smearing(double smear_val, int Nsmear){
+    /**
+     * APE smearing function is used to smooth the gauge field, we take the defination
+     * from Eq. (16) in https://arxiv.org/abs/hep-lat/0409141v2.
+     * 
+     * @param smear_val: Smearing parameter $\alpha$ in the reference
+     * @param Nsmear: Number of smearing iterations.
+     */
+
+    suNg staple, tr1, tr2;
+    int ix, iy, iz, it;
+    int mu, nu, mid, midpmu, midpnu, midmnu, midpmumnu;
+    suNg v, vout, vtmp;
+    
+    double sm1 = 1. - smear_val;  // ( 1 - \alpha )
+    double sm2 = smear_val / 6. ; // ( \alpha \ 6 )
+    
+    suNg_field *u_gauge_tmp = alloc_gfield(&glattice);
+    u_gauge_APE = alloc_gfield(&glattice);
+    
+    #ifdef ALLOCATE_REPR_GAUGE_FIELD
+    u_gauge_APE_f=alloc_gfield_f(&glattice);
+    #endif
+    
+    for (it=0;it<T; it++){
+        for (ix = 0; ix < X; ix++) for (iy = 0; iy < Y; iy++) for (iz = 0; iz < Z; iz++){
+            for (mu = 0; mu < 4; mu++){
+                mid = ipt(it,ix,iy,iz);
+                *pu_gauge_APE(mid, mu) = *pu_gauge(mid, mu);
+            }
+        }
+    }
+    
+    start_gf_sendrecv(u_gauge_APE);
+    complete_gf_sendrecv(u_gauge_APE);
+    represent_gauge_field_APE();
+    
+    lprintf("APE",0,"APE smearing with val=%f \n", smear_val);
+    lprintf("APE",0,"N=0 <p>_s = %1.6f\n",avr_spatial_plaquette_APE());
+    
+    for (int N=0; N<Nsmear; N++){
+        
+        for (it=0;it<T; it++){
+            for (ix = 0; ix < X; ix++) for (iy = 0; iy < Y; iy++) for (iz = 0; iz < Z; iz++){
+                mid = ipt(it, ix, iy, iz);
+                for (mu = 1; mu < 4; mu++){
+                    _suNg_zero(v);
+                    midpmu = iup(mid, mu);
+                
+                    for (int i = 0; i < 2; i++){
+                        nu = (mu + i) % 3 + 1;
+                        midpnu = iup(mid, nu);
+                        midmnu = idn(mid, nu);
+                        midpmumnu = idn(midpmu, nu);
+                    
+                        //Up Staple
+                        _suNg_times_suNg(tr1, *pu_gauge_APE(mid, nu), *pu_gauge_APE(midpnu, mu));
+                        _suNg_times_suNg_dagger(staple, tr1, *pu_gauge_APE(midpmu, nu));
+                        
+                        _suNg_add_assign(v, staple);
+
+                        //Down Staple
+                        _suNg_times_suNg(tr2, *pu_gauge_APE(midmnu, mu), *pu_gauge_APE(midpmumnu, nu));
+                        _suNg_dagger_times_suNg(staple, *pu_gauge_APE(midmnu, nu), tr2);
+
+                        _suNg_add_assign(v, staple);
+                        }
+                
+                    _suNg_mul(tr1, sm1, *pu_gauge_APE(mid, mu));
+                    _suNg_mul(tr2, sm2, v);
+                    
+                    _suNg_zero(vout);
+                    _suNg_add_assign(vout,tr1);
+                    _suNg_add_assign(vout,tr2);
+                    
+                    _suNg_mul(vtmp, 1., vout);
+                    project_to_suNg(&vout); // Project to the group manifold
+                    
+                    #ifdef GAUGE_SPN
+                    cooling_SPN(&v, &vout, &vtmp, 5); // max Re Tr(U V^+)
+                    *pu_gauge_tmp(mid, mu) = v;
+                    
+                    #else
+                    
+                    *pu_gauge_tmp(mid, mu) = vout;
+                    
+                    #endif
+                }
+            }
+        }
+        
+        for (it=0;it<T; it++){
+            for (ix = 0; ix < X; ix++) for (iy = 0; iy < Y; iy++) for (iz = 0; iz < Z; iz++){
+                for (mu = 0; mu < 4; mu++){
+                    mid = ipt(it,ix,iy,iz);
+                    *pu_gauge_APE(mid, mu) = *pu_gauge_tmp(mid, mu);
+                }
+            }
+        }
+        
+        start_gf_sendrecv(u_gauge_APE);
+        complete_gf_sendrecv(u_gauge_APE);
+        represent_gauge_field_APE();
+        
+    }
+    lprintf("APE",0,"N=%d <p>_s = %1.6f\n", Nsmear, avr_spatial_plaquette_APE());
+    polyakov_APE();
+    full_plaquette_APE();
+    lprintf("APE",0,"APE smearing END \n");
+    free_gfield(u_gauge_tmp);
+}
\ No newline at end of file
diff --git a/LibHR/Utils/boundary_conditions.c b/LibHR/Utils/boundary_conditions.c
index 15869eb60..4e11daf68 100644
--- a/LibHR/Utils/boundary_conditions.c
+++ b/LibHR/Utils/boundary_conditions.c
@@ -337,6 +337,33 @@ void apply_BCs_on_clover_term(suNfc_field *cl)
 }
 #endif //defined(GAUGE_SPN) && defined(REPR_FUNDAMENTAL)
 
+static void sp_T_antiperiodic_BCs_APE();
+static void sp_X_antiperiodic_BCs_APE();
+static void sp_Y_antiperiodic_BCs_APE();
+static void sp_Z_antiperiodic_BCs_APE();
+/*static void sp_spatial_theta_BCs(double theta);*/
+static void chiSF_ds_BT_APE(double ds);
+
+void apply_BCs_on_represented_gauge_field_APE() {
+#ifdef BC_T_ANTIPERIODIC
+  sp_T_antiperiodic_BCs_APE();
+#endif
+#ifdef BC_X_ANTIPERIODIC
+  sp_X_antiperiodic_BCs_APE();
+#endif
+#ifdef BC_Y_ANTIPERIODIC
+  sp_Y_antiperiodic_BCs_APE();
+#endif
+#ifdef BC_Z_ANTIPERIODIC
+  sp_Z_antiperiodic_BCs_APE();
+#endif
+#ifdef ROTATED_SF
+#ifndef ALLOCATE_REPR_GAUGE_FIELD
+#error The represented gauge field must be allocated!!!
+#endif
+  chiSF_ds_BT_APE(BCs_pars.chiSF_boundary_improvement_ds);
+#endif
+}
 
 /***************************************************************************/
 /* BOUNDARY CONDITIONS TO BE APPLIED ON THE REPRESENTED GAUGE FIELD        */
@@ -816,6 +843,107 @@ static void gf_open_BCs()
   }
 }
 #endif
+/***************************************************************************/
+/* BOUNDARY CONDITIONS TO BE APPLIED ON THE REPRESENTED APE GAUGE FIELD        */
+/***************************************************************************/
+
+static void sp_T_antiperiodic_BCs_APE() {
+  if(COORD[0]==0) {
+    int index;
+    int ix,iy,iz;
+    suNf *u;
+    for (ix=0;ix<X_EXT;++ix) for (iy=0;iy<Y_EXT;++iy) for (iz=0;iz<Z_EXT;++iz){
+      index=ipt_ext(2*T_BORDER,ix,iy,iz);
+      if(index!=-1) {
+        u=pu_gauge_APE_f(index,0);
+        _suNf_minus(*u,*u);
+      }
+    }
+  }
+}
+
+static void sp_X_antiperiodic_BCs_APE() {
+  if(COORD[1]==0) {
+    int index;
+    int it,iy,iz;
+    suNf *u;
+    for (it=0;it<T_EXT;++it) for (iy=0;iy<Y_EXT;++iy) for (iz=0;iz<Z_EXT;++iz){
+      index=ipt_ext(it,2*X_BORDER,iy,iz);
+      if(index!=-1) {
+        u=pu_gauge_APE_f(index,1);
+        _suNf_minus(*u,*u);
+      }
+    }
+  }
+}
+
+static void sp_Y_antiperiodic_BCs_APE() {
+  if(COORD[2]==0) {
+    int index;
+    int ix,it,iz;
+    suNf *u;
+    for (it=0;it<T_EXT;++it) for (ix=0;ix<X_EXT;++ix) for (iz=0;iz<Z_EXT;++iz){
+      index=ipt_ext(it,ix,2*Y_BORDER,iz);
+      if(index!=-1) {
+        u=pu_gauge_APE_f(index,2);
+        _suNf_minus(*u,*u);
+      }
+    }
+  }
+}
+
+static void sp_Z_antiperiodic_BCs_APE() {
+  if(COORD[3]==0) {
+    int index;
+    int ix,iy,it;
+    suNf *u;
+    for (it=0;it<T_EXT;++it) for (ix=0;ix<X_EXT;++ix) for (iy=0;iy<Y_EXT;++iy){
+      index=ipt_ext(it,ix,iy,2*Z_BORDER);
+      if(index!=-1) {
+        u=pu_gauge_APE_f(index,3);
+        _suNf_minus(*u,*u);
+      }
+    }
+  }
+}
+
+
+static void chiSF_ds_BT_APE(double ds) {
+  if(COORD[0] == 0) {
+    int index;
+    int ix,iy,iz;
+    suNf *u;
+    for (ix=0;ix<X_EXT;++ix) for (iy=0;iy<Y_EXT;++iy) for (iz=0;iz<Z_EXT;++iz){
+      index=ipt_ext(T_BORDER+1,ix,iy,iz);
+      if(index!=-1) {
+        u=pu_gauge_APE_f(index,1);
+        _suNf_mul(*u,ds,*u);
+        u=pu_gauge_APE_f(index,2);
+        _suNf_mul(*u,ds,*u);
+        u=pu_gauge_APE_f(index,3);
+        _suNf_mul(*u,ds,*u);
+      }
+    }
+  }
+  if(COORD[0] == NP_T-1) {
+    int index;
+    int ix,iy,iz;
+    suNf *u;
+    for (ix=0;ix<X_EXT;++ix) for (iy=0;iy<Y_EXT;++iy) for (iz=0;iz<Z_EXT;++iz){
+      index=ipt_ext(T+T_BORDER-1,ix,iy,iz);
+      if(index!=-1) {
+        u=pu_gauge_APE_f(index,1);
+        _suNf_mul(*u,ds,*u);
+        u=pu_gauge_APE_f(index,2);
+        _suNf_mul(*u,ds,*u);
+        u=pu_gauge_APE_f(index,3);
+        _suNf_mul(*u,ds,*u);
+      }
+    }
+  }
+}
+
+
 
 /***************************************************************************/
 /* BOUNDARY CONDITIONS TO BE APPLIED ON THE MOMENTUM FIELDS                */
diff --git a/LibHR/Utils/single_double_utils.c b/LibHR/Utils/single_double_utils.c
index c91ee2de2..eb985e107 100644
--- a/LibHR/Utils/single_double_utils.c
+++ b/LibHR/Utils/single_double_utils.c
@@ -113,3 +113,21 @@ void assign_sd2s(spinor_field_flt *out, spinor_field *in)
     }
   }
 }
+
+/*APE smearing*/
+void assign_ud2u_APE_f(void)
+{
+  if(u_gauge_APE_f_flt!=NULL){
+    double *d;
+    float *f;
+  
+    d = (double*)(u_gauge_APE_f->ptr);
+    f = (float*)(u_gauge_APE_f_flt->ptr);
+
+_OMP_PRAGMA( _omp_parallel )
+_OMP_PRAGMA( _omp_for )
+    for(int i=0; i<4*glattice.gsize_gauge*(sizeof(suNf)/sizeof(double)); i++) {
+      *(f+i) = (float)(*(d+i));
+    }
+  }
+}
\ No newline at end of file
diff --git a/LibHR/Utils/suN_utils.c b/LibHR/Utils/suN_utils.c
index b7ef83e4a..61888f619 100644
--- a/LibHR/Utils/suN_utils.c
+++ b/LibHR/Utils/suN_utils.c
@@ -447,3 +447,214 @@ void covariant_project_to_suNg(suNg *u)
 
   _suNg_times_suNg(*u, tmp1, tmp);
 }
+
+
+void cooling_SPN(suNg* g_out, suNg* g_in, suNg* g_tilde, int cooling){
+  /**
+   * @brief Cooling algorithm for SP(N) gauge fields applies the beta=infinity limit of the Cabibbo-Mariani algorithm,
+   * which gives the maximum value of Re Tr(g_tilde g_in).
+   * 
+   * We follow the procedure described in https://arxiv.org/pdf/hep-lat/0404008, for the SU(N) case.
+   * 
+   * @param g_out: Pointer to the output gauge field after cooling.
+   * @param g_in: Pointer to the input gauge field before cooling, which is U in Eq. (25) of the reference.
+   * @param g_tilde: Pointer to the input gauge field before cooling, which is \tilde{U} in Eq. (25) of the reference.
+   * @param cooling: Number of cooling iterations to perform.
+   * @note `subgrb` and `subgrb_tau` are helper functions to extract the SU(2) sub-matrices, 
+   * while `vmxsu2` and `vmxsu2_tau` are helper functions to multiply the SU(2) sub-matrices.
+   * See detail explanations in Appendex A of the reference: https://arxiv.org/abs/2010.15781.
+   */
+    
+    suNgfull B, U_tilde, S;
+    _suNg_expand( U_tilde, *g_tilde);
+    _suNg_expand( S, *g_in);
+        
+    for (int nbcool=0;nbcool<cooling;nbcool++){
+        _suNgfull_times_suNgfull_dagger(B, S, U_tilde);
+        
+        for (int N1=0;N1<NG/2-1;N1++){
+            for (int N2=N1+1;N2<NG/2;N2++){
+                subgrb(N1, N2, &B, &S);
+                subgrb(NG/2 + N1, NG/2 + N2, &B, &S);
+                subgrb_tau( N1, N2, &B, &S);
+                subgrb(N1, NG/2 + N1, &B, &S);
+                subgrb(N2, NG/2 + N2, &B, &S);
+            }
+        }       
+    }
+    
+    for (int i=0;i<(NG/2)*NG;i++){
+        
+        g_out->c[i] = S.c[i];
+    }
+}
+
+
+void subgrb(int i1col, int i2col, suNgfull* B11, suNgfull* C11){
+  /**
+  * @brief Extracts an SU(2) subgroup from a Sp(N) matrix and applies a transformation.
+  * 
+  * This function targets the SU(2) subgroup defined by the columns `i1col` and `i2col`
+  * of the input Sp(N) matrix `B11`. It constructs an SU(2) of type Eq. (A4) in the paper: https://arxiv.org/abs/2010.15781,
+  * normalizes it, and applies it to both `C11` and `B11` via left multiplication.
+  * 
+  * @param i1col Index of the first column in Sp(N) matrix `B11`.
+  * @param i2col Index of the second column in Sp(N) matrix `B11`.
+  * @param B11 Pointer to the Sp(N) matrix to be updated in-place.
+  * @param C11 Pointer to a Sp(N) matrix also updated by the same SU(2) transformation.
+  */
+    
+    double complex F11, F12, A11[4], ztmp1, ztmp2;    // A11 is the extracted su2 matrix in 1-d array
+    double UMAG;                               // [ 0 1 ]
+    int i1,i2,i3,i4;                           // [ 2 3 ]   --> [0, 1, 2, 3]
+    
+    i1 = i1col + i1col*NG;
+    i2 = i2col + i2col*NG;
+    i3 = i2col + i1col*NG;
+    i4 = i1col + i2col*NG;
+    
+    _complex_add_star(F11, B11->c[i1], B11->c[i2]);
+    _complex_mulr(F11, 0.5, F11);
+    _complex_sub_star(F12, B11->c[i3], B11->c[i4]);
+    _complex_mulr(F12, 0.5, F12);
+    
+    
+    _complex_mul_star(ztmp1, F11, F11);
+    _complex_mul_star(ztmp2, F12, F12);
+    UMAG = sqrt(creal(ztmp1)+creal(ztmp2));
+    UMAG = 1./UMAG;
+    
+    _complex_mulr(F11, UMAG, F11);
+    _complex_mulr(F12, UMAG, F12);
+    
+    _complex_star(A11[0], F11);
+    _complex_star(A11[1], F12);
+    _complex_mulr(A11[2], -1., F12);
+    _complex_mulr(A11[3],  1., F11);
+    
+    vmxsu2(i1col,i2col,C11,A11);
+    vmxsu2(i1col,i2col,B11,A11);
+}
+
+
+void vmxsu2(int i1, int i2, suNgfull* A, double complex B[4]){
+
+  /**
+   * @brief Multiplies columns i1 and i2 of a Sp(N) matrix by a 2×2 SU(2) matrix.
+   *
+   * Applies an SU(2) rotation on columns i1 and i2 of matrix `A`, using matrix `B`
+   * provided in row-major order as a 1D array of 4 complex numbers.
+   * 
+   * @param i1 First column index of Sp(N) matrix A.
+   * @param i2 Second column index of Sp(N) matrix A.
+   * @param A  Pointer to the Sp(N) matrix.
+   * @param B  SU(2) matrix in vector form: {B00, B01, B10, B11}.
+   */
+    
+    double complex C[NG*2];
+    double complex ztmp1,ztmp2;
+    
+    for (int i=0; i<NG; i++){
+        _complex_mul(ztmp1, A->c[i + i1*NG], B[0]);
+        _complex_mul(ztmp2, A->c[i + i2*NG], B[2]);
+        _complex_add(C[i], ztmp1, ztmp2);
+        
+        _complex_mul(ztmp1, A->c[i + i1*NG], B[1]);
+        _complex_mul(ztmp2, A->c[i + i2*NG], B[3]);
+        _complex_add(C[i+NG], ztmp1, ztmp2);
+    }
+    
+    for (int i=0; i<NG; i++){
+        A->c[i + i1*NG] = C[i];
+        A->c[i + i2*NG] = C[i+NG];
+    }
+}
+
+
+void subgrb_tau(int n1, int n2, suNgfull* B11, suNgfull* C11){
+
+  /**
+   * @brief Applies an SU(2) tau-subgroup transformation that mixes upper and lower blocks of Sp(N).
+   *
+   * Constructs a SU(2) matrix of type Eq. (A5) in the paper: https://arxiv.org/abs/2010.15781,
+   * 
+   * Applies the resulting SU(2) matrix to both `B11` and `C11`.
+   *
+   * @param n1 Index in upper block.
+   * @param n2 Index in lower block.
+   * @param B11 Matrix B = S \tilde{U}^\dagger, used to generate the SU(2) transformation.
+   * @param C11 Matrix S, updated with the transformation.
+   */
+    
+    double complex F11, F12, A11[4], ztmp1, ztmp2;    // A11 is the extracted su2 matrix in 1-d array
+    double UMAG;                               // [ 0 1 ]
+    int i1,i2,i3,i4;                           // [ 2 3 ]   --> [0, 1, 2, 3]
+    
+    i1 = n1 + n1*NG;
+    i2 = (NG/2) + n2 + (NG/2 + n2)*NG;
+    i3 = NG/2 + n2 + n1*NG;
+    i4 = (NG/2 + n2)*NG + n1;
+    
+    _complex_add_star(F11, B11->c[i1], B11->c[i2]);
+    _complex_mulr(F11, 0.5, F11);
+    _complex_sub_star(F12, B11->c[i3], B11->c[i4]);
+    _complex_mulr(F12, 0.5, F12);
+    
+    _complex_mul_star(ztmp1, F11, F11);
+    _complex_mul_star(ztmp2, F12, F12);
+    UMAG = sqrt(creal(ztmp1)+creal(ztmp2));
+    UMAG = 1./UMAG;
+    
+    _complex_mulr(F11, UMAG, F11);
+    _complex_mulr(F12, UMAG, F12);
+    
+    _complex_star(A11[0], F11);
+    _complex_star(A11[1], F12);
+    _complex_mulr(A11[2], -1., F12);
+    _complex_mulr(A11[3],  1., F11);
+    
+    vmxsu2_tau(n1, n2, C11, A11);
+    vmxsu2_tau(n1, n2, B11, A11);
+}
+
+
+void vmxsu2_tau(int n1, int n2, suNgfull* A, complex B[4]){
+
+  /**
+   * @brief Applies a (A5)-type SU(2) transformation to a Sp(N) matrix.
+   *
+   * Performs the matrix update defined by a (A5)-type SU(2) matrix.
+   * A temporary copy is used to safely apply the transformation to four sub-blocks:
+   * (n1, n1), (n2, n2), (NG/2+n1, NG/2+n1), and (NG/2+n2, NG/2+n2).
+   * 
+   * @param n1 Index in upper block.
+   * @param n2 Index in lower block.
+   * @param A  Pointer to the Sp(N) matrix to be transformed.
+   * @param B  SU(2) matrix in row-major form.
+   */
+    
+    suNgfull C;
+    double complex ztmp1, ztmp2;
+    
+    _suNgfull_mul(C, 1., *A);
+    
+    for (int i=0; i<NG; i++){
+        
+        _complex_mul(ztmp1, A->c[i + n1*NG], B[0]);
+        _complex_mul(ztmp2, A->c[i +(NG/2+n2)*NG], B[2]);
+        _complex_add(C.c[i + n1*NG], ztmp1, ztmp2);
+        
+        _complex_mul(ztmp1, A->c[i + n2*NG], B[0]);
+        _complex_mul(ztmp2, A->c[i + (NG/2+n1)*NG], B[2]);
+        _complex_add(C.c[i + n2*NG], ztmp1, ztmp2);
+        
+        _complex_mul(ztmp1, A->c[i + n2*NG], B[1]);
+        _complex_mul(ztmp2, A->c[i + (NG/2+n1)*NG], B[3]);
+        _complex_add(C.c[i + (NG/2 + n1)*NG], ztmp1, ztmp2);
+        
+        _complex_mul(ztmp1, A->c[i + n1*NG], B[1]);
+        _complex_mul(ztmp2, A->c[i + (NG/2+n2)*NG], B[3]);
+        _complex_add(C.c[i + (NG/2 + n2)*NG], ztmp1, ztmp2);
+    }
+    _suNgfull_mul(*A, 1., C);
+}
diff --git a/Spectrum/Makefile b/Spectrum/Makefile
index f4b03e703..c324409b1 100644
--- a/Spectrum/Makefile
+++ b/Spectrum/Makefile
@@ -2,7 +2,7 @@ TOPDIR = ..
 MKDIR = $(TOPDIR)/Make
 
 #EXES =  random_cnfg random_spinor mk_mesons_with_z2semwall mk_mesons_with_z2semwall_new  measure_spectrum measure_formfactor #mk_sfcoupling #trunc_mesons mk_mesons
-EXES=measure_spectrum mk_sfcorrelators
+EXES=measure_spectrum
 
 
 
diff --git a/Spectrum/input_file_smeared b/Spectrum/input_file_smeared
new file mode 100644
index 000000000..3e8776e62
--- /dev/null
+++ b/Spectrum/input_file_smeared
@@ -0,0 +1,74 @@
+// Global variables
+GLB_T = 8
+GLB_X = 4
+GLB_Y = 4
+GLB_Z = 4
+NP_T = 2
+NP_X = 1
+NP_Y = 1
+NP_Z = 1
+
+// Replicas
+N_REP = 1
+
+// Random generator
+rlx_level = 1
+rlx_seed = 11228
+rlx_start = new
+rlx_state = ./out/rlx_state
+
+log:default = -1
+log:inverter = -1
+log:forcestat = 0
+
+mes:precision = 1.e-14
+mes:nhits = 1
+mes:nhits_disc = 4
+mes:masses = -0.85
+mes:use_input_mass = 0
+mes:meas_mixed = 0
+mes:nhits_2pt = 1
+mes:momentum = 0
+mes:def_semwall = 0
+mes:def_semwall_nondeg = 0
+mes:def_point = 0
+mes:def_baryon = 0
+mes:def_gfwall = 0
+mes:ext_semwall = 0
+mes:ext_point = 0
+mes:dirichlet_semwall = 0
+mes:dirichlet_point = 0
+mes:dirichlet_gfwall = 0
+mes:discon_semwall = 0
+mes:discon_gfwall = 0
+mes:discon_volume = 0
+mes:dilution = 1
+mes:dirichlet_dt = 2
+mes:configlist = listfile_spectrum
+
+// Measure smeared meson spectrum
+// enable source and sink smearing
+mes:smearing_source_sink = 1
+mes:discon_semwall_smeared = 1
+
+// point source location for source ans sink smearing
+mes:source_t = 0
+mes:source_x = 0
+mes:source_y = 0
+mes:source_z = 0
+
+// source and sink smearing parameters
+mes:smear_epsilon_source = 0.2
+mes:smear_N_source = 50
+mes:smear_epsilon_sink = 0.2
+mes:smear_N_sink = 50
+mes:smear_N_step = 10
+//APE smearing activate if APE_N != 0
+mes:APE_epsilon = 0.4
+mes:APE_N = 10
+
+//Whether to include the four fermion interaction
+ff:on = false
+//Parameters for hopping parameter expansion, 0 in degree_hopping to disable
+mes:degree_hopping = 0
+mes:nhits_hopping = 5
diff --git a/Spectrum/measure_spectrum.c b/Spectrum/measure_spectrum.c
index c8879478c..e7a48b71a 100755
--- a/Spectrum/measure_spectrum.c
+++ b/Spectrum/measure_spectrum.c
@@ -62,6 +62,7 @@ typedef struct _input_mesons
   int fixed_point;
   int fixed_gfwall;
   int discon_semwall;
+  int discon_semwall_smeared;
   int discon_gfwall;
   int discon_volume;
   int def_gfwall;
@@ -74,13 +75,25 @@ typedef struct _input_mesons
   double csw;
   double rho_s;
   double rho_t;
+  int source_t;
+  int source_x;
+  int source_y;
+  int source_z;
+  int smearing_source_sink;
+  double smear_epsilon_source;
+  int smear_N_source;
+  double smear_epsilon_sink;
+  int smear_N_sink;
+  int smear_N_step;
+  int APE_N;
+  double APE_epsilon;
 
   //Currently only implemented for ff
   int nhits_hopping;  //Multiplies the number of hits in the fast part of the hopping parameter expansion
   int degree_hopping; // The degree of the hopping parameter expasion
 
   /* for the reading function */
-  input_record_t read[33];
+  input_record_t read[44];
 } input_mesons;
 
 #define init_input_mesons(varname)                                                                                \
@@ -101,6 +114,7 @@ typedef struct _input_mesons
       {"enable Dirichlet point", "mes:dirichlet_point = %d", INT_T, &(varname).fixed_point},                      \
       {"enable Dirichlet gfwall", "mes:dirichlet_gfwall = %d", INT_T, &(varname).fixed_gfwall},                   \
       {"enable discon semwall", "mes:discon_semwall = %d", INT_T, &(varname).discon_semwall},                     \
+      {"enable discon semwall with smearing", "mes:discon_semwall_smeared = %d",INT_T, &(varname).discon_semwall_smeared},\
       {"enable discon gfwall", "mes:discon_gfwall = %d", INT_T, &(varname).discon_gfwall},                        \
       {"enable discon volume", "mes:discon_volume = %d", INT_T, &(varname).discon_volume},                        \
       {"volume source dilution", "mes:dilution = %d", INT_T, &(varname).dilution},                                \
@@ -117,6 +131,18 @@ typedef struct _input_mesons
       {"hopping expansion degree", "mes:degree_hopping = %d", INT_T, &(varname).degree_hopping},                  \
       {"hopping expansion hits", "mes:nhits_hopping = %d", INT_T, &(varname).nhits_hopping},                      \
       {"Configuration list:", "mes:configlist = %s", STRING_T, &(varname).configlist},                            \
+      {"Smearing Source z", "mes:source_t = %d",INT_T, &(varname).source_t},                                      \
+      {"Smearing Source x", "mes:source_x = %d",INT_T, &(varname).source_x},                                      \
+      {"Smearing Source y", "mes:source_y = %d",INT_T, &(varname).source_y},                                      \
+      {"Smearing Source z", "mes:source_z = %d",INT_T, &(varname).source_z},                                      \
+      {"Smearing a lot at once", "mes:smearing_source_sink = %d",INT_T, &(varname).smearing_source_sink},         \
+      {"smearing source step size", "mes:smear_epsilon_source = %lf", DOUBLE_T, &(varname).smear_epsilon_source}, \
+      {"Number of smeared source steps", "mes:smear_N_source = %d",INT_T, &(varname).smear_N_source},             \
+      {"smearing sink step size", "mes:smear_epsilon_sink = %lf", DOUBLE_T, &(varname).smear_epsilon_sink},       \
+      {"Number of smeared sink steps", "mes:smear_N_sink = %d",INT_T, &(varname).smear_N_sink},                   \
+      {"Number of steps between smearing levels", "mes:smear_N_step = %d",INT_T, &(varname).smear_N_step},        \
+      {"Number of APE smearing steps", "mes:APE_N = %d",INT_T, &(varname).APE_N},                                 \
+      {"APE smearing step size", "mes:APE_epsilon = %lf",DOUBLE_T, &(varname).APE_epsilon},                       \
       {NULL, NULL, INT_T, NULL}                                                                                   \
     }                                                                                                             \
   }
@@ -208,6 +234,12 @@ int main(int argc, char *argv[])
   {
     lprintf("MAIN", 0, "Point sources\n");
   }
+  if (mes_var.APE_N != 0){
+      lprintf("MAIN",0,"APE smearing activate on Gaussian smearing\n");
+  }
+  if (mes_var.smearing_source_sink){
+    lprintf("MAIN",0,"Smear source to several Ns_sink\n");
+  }
   if (mes_var.def_baryon)
   {
     lprintf("MAIN", 0, "Baryon masses\n");
@@ -330,6 +362,9 @@ int main(int argc, char *argv[])
       {
         measure_spectrum_pt(tau, nm, m, mes_var.n_mom, i, mes_var.precision,DONTSTORE, NULL);
       }
+      if (mes_var.smearing_source_sink){
+        measure_smearing_source_sink(mes_var.source_t,mes_var.source_x,mes_var.source_y,mes_var.source_z,nm,m,mes_var.n_mom,mes_var.nhits_2pt,i,mes_var.precision,mes_var.smear_epsilon_source,mes_var.smear_N_source,mes_var.smear_epsilon_sink,mes_var.smear_N_sink, mes_var.APE_epsilon, mes_var.APE_N, mes_var.smear_N_step);
+      }
       if (mes_var.def_baryon)
       {
 #if NG == 3
@@ -370,6 +405,9 @@ int main(int argc, char *argv[])
       {
         measure_spectrum_discon_semwall(nm, m, mes_var.nhits_disc, i, mes_var.precision,DONTSTORE, NULL);
       }
+      if (mes_var.discon_semwall_smeared){
+        measure_spectrum_discon_semwall_smeared(nm,m,mes_var.nhits_disc,i,mes_var.precision,mes_var.smear_epsilon_source,mes_var.smear_N_source, mes_var.APE_epsilon, mes_var.APE_N, mes_var.smear_N_step);
+      }
       if (mes_var.discon_gfwall)
       {
         //       measure_spectrum_discon_gfwall(nm,m,i,mes_var.precision,DONTSTORE, NULL);