diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py
index 7874d10b13b6b0414323b2d7fd72aa66fac2bd44..e212ce51575ac175f668c67396e2a6ded7249d75 100755
--- a/ObservationSim/MockObject/Galaxy.py
+++ b/ObservationSim/MockObject/Galaxy.py
@@ -87,8 +87,6 @@ class Galaxy(MockObject):
             gal = gal.shear(gal_shear)
             gal = galsim.Convolve(psf, gal)
             
-            # if fd_shear is not None:
-            #     gal = gal.shear(fd_shear)
             objs.append(gal)
         final = galsim.Sum(objs)
         return final
@@ -143,6 +141,11 @@ class Galaxy(MockObject):
         bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
         bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
         bulge = bulge.shear(bulge_shape)
+        
+        if fd_shear:
+            g1 += fd_shear.g1
+            g2 += fd_shear.g2
+        gal_shear = galsim.Shear(g1=g1, g2=g2)
 
         for i in range(len(bandpass_list)):
             bandpass = bandpass_list[i]
@@ -167,14 +170,9 @@ class Galaxy(MockObject):
             # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
 
             psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
-            # disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
-            # disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
-            # disk = disk.shear(disk_shape)
-            # bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
-            # bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
-            # bulge = bulge.shear(bulge_shape)
             
             gal_temp = self.bfrac * bulge + (1.0 - self.bfrac) * disk
+            gal_temp = gal_temp.shear(gal_shear)
 
             gal_temp = gal_temp.withFlux(nphotons)
             if not big_galaxy: # Not apply PSF for very big galaxy
@@ -190,27 +188,6 @@ class Galaxy(MockObject):
             # kfrac = np.random.random()*(1.0 - self.bfrac)
             # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
 
-            # gal = gal.withFlux(nphotons)
-            # gal_shear = galsim.Shear(g1=g1, g2=g2)
-            # gal = gal.shear(gal_shear)
-
-            # if not big_galaxy: # Not apply PSF for very big galaxy
-            #     gal = galsim.Convolve(psf, gal)
-            #     # if fd_shear is not None:
-            #     #     gal = gal.shear(fd_shear)
-
-            # # Use (explicit) stamps to draw
-            # stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True)
-            
-            # xmax = max(xmax, stamp.xmax - stamp.xmin)
-            # ymax = max(ymax, stamp.ymax - stamp.ymin)
-
-            # photons = stamp.photons
-            # photons.x += x_nominal
-            # photons.y += y_nominal
-            # photons_list.append(photons)
-            # del gal
-
         # # [C6 TEST]
         # print('xmax = %d, ymax = %d '%(xmax, ymax))
         # # Output memory usage
@@ -218,13 +195,6 @@ class Galaxy(MockObject):
         # top_stats = snapshot.statistics('lineno')
         # for stat in top_stats[:10]:
         #     print(stat)
-        if fd_shear:
-            g1 += fd_shear.g1
-            g2 += fd_shear.g2
-        gal_shear = galsim.Shear(g1=g1, g2=g2)
-        gal = gal.shear(gal_shear)
-        # if fd_shear is not None:
-        #     gal = gal.shear(fd_shear)
 
         stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True)
         photons = stamp.photons
@@ -232,7 +202,6 @@ class Galaxy(MockObject):
         photons.y += y_nominal
         photons_list.append(photons)
 
-        # stamp = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.1))
         stamp.wcs = real_wcs_local
         stamp.setCenter(x_nominal, y_nominal)
         bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py
index 28c097e64a3f4fb1fa834e62bb3a2513cca1c070..a70f99788a1f5a955072670dcbd17ce14defa4e4 100755
--- a/ObservationSim/ObservationSim.py
+++ b/ObservationSim/ObservationSim.py
@@ -220,7 +220,7 @@ class Observation(object):
                     continue
                 
                 # [TODO] Testing
-                chip_output.Log_info("mag_%s = %.3f"%(filt.filter_type.lower(), obj.param["mag_%s"%filt.filter_type.lower()]))
+                # chip_output.Log_info("mag_%s = %.3f"%(filt.filter_type.lower(), obj.param["mag_%s"%filt.filter_type.lower()]))
 
                 # Exclude very bright/dim objects (for now)
                 if cut_filter.is_too_bright(
diff --git a/ObservationSim/PSF/FieldDistortion.py b/ObservationSim/PSF/FieldDistortion.py
index 1296aa6d02be4c814bfbf96352c513a58e0a9881..7c95fd6b7b3369cb2ae5792569dc293989d02a6c 100644
--- a/ObservationSim/PSF/FieldDistortion.py
+++ b/ObservationSim/PSF/FieldDistortion.py
@@ -102,12 +102,8 @@ class FieldDistortion(object):
         g_abs = np.sqrt(g1k_fd**2 + g2k_fd**2)
         phi = cmath.phase(complex(g1k_fd, g2k_fd))
         # g_abs = 0.7
-        g1k_fd = g_abs * np.cos(phi - 2*img_rot)
-        g2k_fd = g_abs * np.sin(phi - 2*img_rot)
-        # g1k_fd = g_abs * np.cos(phi - 2*img_rot)
-        # g2k_fd = -g_abs * np.sin(phi - 2*img_rot)
-        # g1k_fd = g_abs * np.cos(0. - 2*img_rot)
-        # g2k_fd = g_abs * np.sin(0. - 2*img_rot)
+        g1k_fd = g_abs * np.cos(phi + 2*img_rot)
+        g2k_fd = g_abs * np.sin(phi + 2*img_rot)
 
         fd_shear = galsim.Shear(g1=g1k_fd, g2=g2k_fd)
         return galsim.PositionD(x, y), fd_shear