diff --git a/Triangle/Data.py b/Triangle/Data.py
index e07716c..7f214d2 100644
--- a/Triangle/Data.py
+++ b/Triangle/Data.py
@@ -4,6 +4,7 @@
 import scipy.sparse
 from scipy.integrate import cumulative_trapezoid
 from scipy.signal import firwin, kaiserord, lfilter
+from collections import UserDict
 import h5py
 
 from Triangle.Constants import *
@@ -64,8 +65,6 @@ def SliceDataDictionary(data_dict, time, time_start, time_end, buffer_time=80, c
 
 
 # ================ methods for general dict =============================
-
-
 def Dict2Array(d, keys=None):
     """
     the generated array must follow the order given by keys
@@ -121,160 +120,170 @@ def MultiplyDicts(dict1, dict2, coef=1.0):
 
 
 # ========================== MOSA/SC dictionary wrapper ============================
-class MOSADict:
-    """
-    NOTE: be careful when manipulating array or dict elements
-    NOTE: use B = A.copy() instead of B = A
-    the keys of MOSADict must be MOSA_labels
-    the values of MOSADict are numpy arrays of the same dim or scalars
-    """
-
-    dict_keys = MOSA_labels
-
-    def __init__(self, data_dict):
-        if isinstance(data_dict, dict):
-            if set(data_dict.keys()) == set(self.dict_keys):
-                self.data = data_dict
-            else:
-                raise ValueError("keys of dictionaries mismatch.")
-        else:
-            raise TypeError("data_dict must be a dict.")
-
-    def __repr__(self):
-        return str(self.data)
-
-    def __getitem__(self, key):
-        return self.data[key]
+class ObjectDict(UserDict):
+    """A dictionary holding values for object."""
+    _allowed_keys = []
 
     def __setitem__(self, key, value):
-        if key in self.dict_keys:
-            self.data[key] = value
+        if key in self._allowed_keys:
+            super().__setitem__(key, value)
         else:
-            raise ValueError(str(key) + " is not a surpported key.")
-
+            raise ValueError(str(key) + " is not a supported key.")
+        
+    def __array__(self, dtype=object):
+        return np.array(self.data, dtype=dtype)
+            
     def __add__(self, other):
-        if isinstance(other, MOSADict):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] + other.data[key]
-            return MOSADict(new_data)
+        new_dict = {}
+        if type(self) == type(other):
+            for key in self.keys():
+                new_dict[key] = self[key] + other[key]
         elif np.isscalar(other):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] + other
-            return MOSADict(new_data)
+            for key in self.keys():
+                new_dict[key] = self[key] + other
         else:
-            raise TypeError("Unsupported operation between instances of 'MOSADict' and '{}'".format(type(other)))
-
+            raise TypeError(f"Unsupported operation between '{type(self)}' and '{type(other)}'.")
+        return self.__class__(new_dict)
+    
     def __radd__(self, other):
         return self.__add__(other)
-
+    
     def __mul__(self, other):
-        if isinstance(other, MOSADict):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] * other.data[key]
-            return MOSADict(new_data)
+        new_dict = {}
+        if type(self) == type(other):
+            for key in self.keys():
+                new_dict[key] = self[key] * other[key]
         elif np.isscalar(other):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] * other
-            return MOSADict(new_data)
+            for key in self.keys():
+                new_dict[key] = self[key] * other
         else:
-            raise TypeError("Unsupported operation between instances of 'MOSADict' and '{}'".format(type(other)))
+            raise TypeError(f"Unsupported operation between '{type(self)}' and '{type(other)}'.")
+        return self.__class__(new_dict)
 
     def __rmul__(self, other):
         return self.__mul__(other)
 
     def __sub__(self, other):
-        if isinstance(other, MOSADict):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] - other.data[key]
-            return MOSADict(new_data)
+        new_dict = {}
+        if type(self) == type(other):
+            for key in self.keys():
+                new_dict[key] = self[key] - other[key]
         elif np.isscalar(other):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] - other
-            return MOSADict(new_data)
+            for key in self.keys():
+                new_dict[key] = self[key] - other
         else:
-            raise TypeError("Unsupported operation between instances of 'MOSADict' and '{}'".format(type(other)))
+            raise TypeError(f"Unsupported operation between '{type(self)}' and '{type(other)}'.")
+        return self.__class__(new_dict)
 
     def __rsub__(self, other):
-        if isinstance(other, MOSADict):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = other.data[key] - self.data[key]
-            return MOSADict(new_data)
+        new_dict = {}
+        if type(self) == type(other):
+            for key in self.keys():
+                new_dict[key] = other[key] - self[key]
         elif np.isscalar(other):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = other - self.data[key]
-            return MOSADict(new_data)
+            for key in self.keys():
+                new_dict[key] = other - self[key]
         else:
-            raise TypeError("Unsupported operation between instances of 'MOSADict' and '{}'".format(type(other)))
-
+            raise TypeError(f"Unsupported operation between '{type(self)}' and '{type(other)}'.")
+        return self.__class__(new_dict)
+        
     def __neg__(self):
-        return MOSADict({key: -self.data[key] for key in self.dict_keys})
+        return self.__class__({key: -self[key] for key in self.keys()})
 
     def __truediv__(self, other):
-        if isinstance(other, MOSADict):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] / other.data[key]
-            return MOSADict(new_data)
+        new_dict = {}
+        if type(self) == type(other):
+            for key in self.keys():
+                new_dict[key] = self[key] / other[key]
         elif np.isscalar(other):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] / other
-            return MOSADict(new_data)
+            for key in self.keys():
+                new_dict[key] = self[key] / other
         else:
-            raise TypeError("Unsupported operation between instances of 'MOSADict' and '{}'".format(type(other)))
+            raise TypeError(f"Unsupported operation between '{type(self)}' and '{type(other)}'.")
+        return self.__class__(new_dict)
 
     def __rtruediv__(self, other):
-        if isinstance(other, MOSADict):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = other.data[key] / self.data[key]
-            return MOSADict(new_data)
+        new_dict = {}
+        if type(self) == type(other):
+            for key in self.keys():
+                new_dict[key] = other[key] / self[key]
         elif np.isscalar(other):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = other / self.data[key]
-            return MOSADict(new_data)
+            for key in self.keys():
+                new_dict[key] = other / self[key]
+        else:
+            raise TypeError(f"Unsupported operation between '{type(self)}' and '{type(other)}'.")
+        return self.__class__(new_dict)
+    
+    def copy(self):
+        if type(self[self._allowed_keys[0]]) is np.ndarray:
+            d = {}
+            for key in self._allowed_keys:
+                d[key] = np.copy(self.data[key])
+            return self.__class__(d)
         else:
-            raise TypeError("Unsupported operation between instances of 'MOSADict' and '{}'".format(type(other)))
+            return self.__class__(super().copy())
 
-    def keys(self):
-        return self.data.keys()
+    def integrate(self, fsample):
+        d = {}
+        for k in self.keys():
+            d[k] = cumulative_trapezoid(np.insert(self[k], 0, 0), dx=1 / fsample)
+        return self.__class__(d)
 
-    def values(self):
-        return self.data.values()
+    def derivative(self, fsample):
+        d = {}
+        for k in self.keys():
+            d[k] = np.gradient(self[k], 1 / fsample)
+        return self.__class__(d)
 
-    def items(self):
-        return self.data.items()
+    def downsampled(self, fsample, downsample, kaiser_filter_coef):
+        d = {}
+        for k in self.keys():
+            d[k] = downsampling(self[k], fsample, downsample, kaiser_filter_coef)
+        return self.__class__(d)
 
-    def copy(self):
-        if type(self.data[MOSA_labels[0]]) is np.ndarray:
-            d = {}
-            for key in MOSA_labels:
-                d[key] = np.copy(self.data[key])
-            return MOSADict(d)
+    def detrended(self, order=2):
+        d = {}
+        for k in self.keys():
+            d[k] = detrending(self[k], order=order)
+        return self.__class__(d)
+
+    def drop_edge_points(self, points1=0, points2=0):
+        d = {}
+        if points2 == 0:
+            for k in self.keys():
+                d[k] = self[k][points1:]
         else:
-            return MOSADict(self.data.copy())
+            for k in self.keys():
+                d[k] = self[k][points1:-points2]
+        return self.__class__(d)
+
+    def get_data_by_idx(self, idx):
+        d = {}
+        for key in self.keys():
+            d[key] = self[key][idx]
+        return self.__class__(d)
+
+class MOSADict(ObjectDict):
+    """
+    NOTE: be careful when manipulating array or dict elements
+    NOTE: use B = A.copy() instead of B = A
+    the keys of MOSADict must be MOSA_labels
+    the values of MOSADict are numpy arrays of the same dim or scalars
+    """
+    _allowed_keys = MOSA_labels
 
     def reverse(self):
         d = {}
-        for k in MOSA_labels:
+        for k in self.keys():
             k_re = k[1] + k[0]
-            d[k] = self.data[k_re]
+            d[k] = self[k_re]
         return MOSADict(d)
 
     def adjacent(self):
         d = {}
-        for k in MOSA_labels:
+        for k in self.keys():
             k_ad = adjacent_MOSA_labels[k]
-            d[k] = self.data[k_ad]
+            d[k] = self[k_ad]
         return MOSADict(d)
 
     def timedelay(self, fsample, delay, doppler=None, order=31, pool=None):
@@ -286,10 +295,10 @@ def timedelay(self, fsample, delay, doppler=None, order=31, pool=None):
             if isinstance(doppler, MOSADict):
                 if pool is None:
                     for k in MOSA_labels:
-                        d[k] = timeshift(self.data[k], -delay[k] * fsample, order=order)
+                        d[k] = timeshift(self[k], -delay[k] * fsample, order=order)
                         # d[k] = timeshift(self.data[k], -delay[k] * fsample, order=order) * (1. - doppler[k])
                 else:
-                    param_arr = [(self.data[k], -delay[k] * fsample, order) for k in MOSA_labels]
+                    param_arr = [(self[k], -delay[k] * fsample, order) for k in MOSA_labels]
                     d_arr = pool.starmap(timeshift, param_arr)
                     for j, k in enumerate(MOSA_labels):
                         d[k] = d_arr[j]
@@ -297,219 +306,27 @@ def timedelay(self, fsample, delay, doppler=None, order=31, pool=None):
             elif doppler is None:
                 if pool is None:
                     for k in MOSA_labels:
-                        d[k] = timeshift(self.data[k], -delay[k] * fsample, order=order)
+                        d[k] = timeshift(self[k], -delay[k] * fsample, order=order)
                 else:
-                    param_arr = [(self.data[k], -delay[k] * fsample, order) for k in MOSA_labels]
+                    param_arr = [(self[k], -delay[k] * fsample, order) for k in MOSA_labels]
                     d_arr = pool.starmap(timeshift, param_arr)
                     for j, k in enumerate(MOSA_labels):
                         d[k] = d_arr[j]
                 return MOSADict(d)
             else:
-                raise TypeError("unsurpported doppler type.")
-        else:
-            raise TypeError("unsurpported delay type.")
-
-    def integrate(self, fsample):
-        d = {}
-        for k in MOSA_labels:
-            d[k] = cumulative_trapezoid(np.insert(self.data[k], 0, 0), dx=1 / fsample)
-        return MOSADict(d)
-
-    def derivative(self, fsample):
-        d = {}
-        for k in MOSA_labels:
-            d[k] = np.gradient(self.data[k], 1 / fsample)
-        return MOSADict(d)
-
-    def downsampled(self, fsample, downsample, kaiser_filter_coef):
-        d = {}
-        for k in MOSA_labels:
-            d[k] = downsampling(self.data[k], fsample, downsample, kaiser_filter_coef)
-        return MOSADict(d)
-
-    def detrended(self, order=2):
-        d = {}
-        for k in MOSA_labels:
-            d[k] = detrending(self.data[k], order=order)
-        return MOSADict(d)
-
-    def drop_edge_points(self, points1=0, points2=0):
-        d = {}
-        if points2 == 0:
-            for k in MOSA_labels:
-                d[k] = self.data[k][points1:]
-        else:
-            for k in MOSA_labels:
-                d[k] = self.data[k][points1:-points2]
-        return MOSADict(d)
-
-    def get_data_by_idx(self, idx):
-        d = {}
-        for key in MOSA_labels:
-            d[key] = self.data[key][idx]
-        return MOSADict(d)
-
-
-class SCDict:
-    """
-    NOTE: be careful when manipulating array or dict elements
-    NOTE: use B = A.copy() instead of B = A
-    the keys of SCDict must be SC_labels
-    the values of SCDict are numpy arrays of the same dim or scalars
-    """
-
-    dict_keys = SC_labels
-
-    def __init__(self, data_dict):
-        if isinstance(data_dict, dict):
-            if set(data_dict.keys()) == set(self.dict_keys):
-                self.data = data_dict
-            else:
-                raise ValueError("keys of dictionaries mismatch.")
-        else:
-            raise TypeError("data_dict must be a dict.")
-
-    def __repr__(self):
-        return str(self.data)
-
-    def __getitem__(self, key):
-        return self.data[key]
-
-    def __setitem__(self, key, value):
-        if key in self.dict_keys:
-            self.data[key] = value
-        else:
-            raise ValueError(str(key) + " is not a surpported key.")
-
-    def __add__(self, other):
-        if isinstance(other, SCDict):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] + other.data[key]
-            return SCDict(new_data)
-        elif np.isscalar(other):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] + other
-            return SCDict(new_data)
-        else:
-            raise TypeError("Unsupported operation between instances of 'SCDict' and '{}'".format(type(other)))
-
-    def __radd__(self, other):
-        return self.__add__(other)
-
-    def __mul__(self, other):
-        if isinstance(other, SCDict):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] * other.data[key]
-            return SCDict(new_data)
-        elif np.isscalar(other):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] * other
-            return SCDict(new_data)
+                raise TypeError("unsupported doppler type.")
         else:
-            raise TypeError("Unsupported operation between instances of 'SCDict' and '{}'".format(type(other)))
-
-    def __rmul__(self, other):
-        return self.__mul__(other)
-
-    def __sub__(self, other):
-        if isinstance(other, SCDict):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] - other.data[key]
-            return SCDict(new_data)
-        elif np.isscalar(other):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] - other
-            return SCDict(new_data)
-        else:
-            raise TypeError("Unsupported operation between instances of 'SCDict' and '{}'".format(type(other)))
-
-    def __rsub__(self, other):
-        if isinstance(other, SCDict):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = other.data[key] - self.data[key]
-            return SCDict(new_data)
-        elif np.isscalar(other):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = other - self.data[key]
-            return SCDict(new_data)
-        else:
-            raise TypeError("Unsupported operation between instances of 'SCDict' and '{}'".format(type(other)))
-
-    def __neg__(self):
-        return SCDict({key: -self.data[key] for key in self.dict_keys})
-
-    def __truediv__(self, other):
-        if isinstance(other, SCDict):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] / other.data[key]
-            return SCDict(new_data)
-        elif np.isscalar(other):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = self.data[key] / other
-            return SCDict(new_data)
-        else:
-            raise TypeError("Unsupported operation between instances of 'SCDict' and '{}'".format(type(other)))
-
-    def __rtruediv__(self, other):
-        if isinstance(other, SCDict):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = other.data[key] / self.data[key]
-            return SCDict(new_data)
-        elif np.isscalar(other):
-            new_data = {}
-            for key in self.dict_keys:
-                new_data[key] = other / self.data[key]
-            return SCDict(new_data)
-        else:
-            raise TypeError("Unsupported operation between instances of 'SCDict' and '{}'".format(type(other)))
-
-    def keys(self):
-        return self.data.keys()
-
-    def values(self):
-        return self.data.values()
-
-    def items(self):
-        return self.data.items()
-
-    def copy(self):
-        if type(self.data[SC_labels[0]]) is np.ndarray:
-            d = {}
-            for key in SC_labels:
-                d[key] = np.copy(self.data[key])
-            return SCDict(d)
-        else:
-            return SCDict(self.data.copy())
+            raise TypeError("unsupported delay type.")
+    
+class SCDict(ObjectDict):
+    _allowed_keys = SC_labels
 
     def toMOSA(self):
         d = {}
         for key in MOSA_labels:
-            d[key] = self.data[key[0]]
+            d[key] = self[key[0]]
         return MOSADict(d)
-
-    def integrate(self, fsample):
-        d = {}
-        for k in SC_labels:
-            d[k] = cumulative_trapezoid(np.insert(self.data[k], 0, 0), dx=1 / fsample)
-        return SCDict(d)
-
-    def derivative(self, fsample):
-        d = {}
-        for k in SC_labels:
-            d[k] = np.gradient(self.data[k], 1 / fsample)
-        return SCDict(d)
-
+    
     def timedelay(self, fsample, delay, order=31, pool=None):
         """
         delay must be SCDicts
@@ -526,31 +343,8 @@ def timedelay(self, fsample, delay, order=31, pool=None):
                     d[k] = d_arr[j]
             return SCDict(d)
         else:
-            raise TypeError("unsurpported delay type.")
-
-    def downsampled(self, fsample, downsample, kaiser_filter_coef):
-        d = {}
-        for k in SC_labels:
-            d[k] = downsampling(self.data[k], fsample, downsample, kaiser_filter_coef)
-        return SCDict(d)
-
-    def detrended(self, order=2):
-        d = {}
-        for k in SC_labels:
-            d[k] = detrending(self.data[k], order=order)
-        return SCDict(d)
-
-    def drop_edge_points(self, points1=0, points2=0):
-        d = {}
-        if points2 == 0:
-            for k in SC_labels:
-                d[k] = self.data[k][points1:]
-        else:
-            for k in SC_labels:
-                d[k] = self.data[k][points1:-points2]
-        return SCDict(d)
-
-
+            raise TypeError("unsupported delay type.")
+        
 def assign_noise_for_MOSAs(arrays_or_psds, fsample, size):
     """
     if arrays_or_psds is numpy array or psd function, it will be filled in to a dictionary (arrays are copied)
@@ -835,8 +629,9 @@ def timeshift(data, shifts, order=31):
 
 def store_dict_to_h5(h5parent, dictitem):
     for k, v in dictitem.items():
-        if isinstance(v, dict):
+        if isinstance(v, (dict, UserDict)):
             group = h5parent.create_group(k)
+            group.attrs["type"] = v.__class__.__name__
             store_dict_to_h5(group, v)
         else:
             if isinstance(v, (str, int, float)):
@@ -847,7 +642,15 @@ def store_dict_to_h5(h5parent, dictitem):
 
 
 def read_dict_from_h5(h5parent):
-    result = {}
+    try:
+        type_name = h5parent.attrs["type"]
+        if type_name == "MOSADict":
+            result = MOSADict({})
+        elif type_name == "SCDict":
+            result = SCDict({})
+        else:   result = {}
+    except KeyError:    result = {}
+    
     for k, v in h5parent.items():
         if isinstance(v, h5py.Group):
             result[k] = read_dict_from_h5(v)
diff --git a/Triangle/Noise.py b/Triangle/Noise.py
index f841c45..3a23c3e 100644
--- a/Triangle/Noise.py
+++ b/Triangle/Noise.py
@@ -786,9 +786,9 @@ def PSD_ACC(self, f, AACC=SACC_nominal):
     def TDI_P_ij(self, P_ij_strings, f):
         """
         Args:
-            P_ij_strings is a MosaDict, each item like: [(1., ["12", "21"]), (-1., ["-13"]), ...], repersenting [D_12 D_21 - A_13 ...]
+            P_ij_strings is a Dict, each item like: [(1., ["12", "21"]), (-1., ["-13"]), ...], representing [D_12 D_21 - A_13 ...]
         Returns:
-            P_ij: MosaDict, each item (Nf)
+            P_ij: Dict, each item (Nf)
         """
         # calculate single delay operators
         tmp = -1.0j * TWOPI * f  # -2pif
diff --git a/Triangle/TDI.py b/Triangle/TDI.py
index a25a6af..9a5273b 100644
--- a/Triangle/TDI.py
+++ b/Triangle/TDI.py
@@ -254,13 +254,11 @@ def CalculatePathCombination(self, plus_paths, minus_paths, eta_dict, doppler=Tr
             logger.debug("Calculating minus path " + minus_paths[i])
             tdi_minus += self.CalculatePath(Path_string=minus_paths[i], eta_dict=eta_dict, doppler=doppler)
         return tdi_plus - tdi_minus
-
-    def CalculateNestedDelay(self, Delay_string, measurement, doppler=True):
+    
+    def CalculateNestedDelayTime(self, Delay_string):
         """
-        calculate nested delay of any measurement
+        calculate the total delay time of a multiple delay operator
         """
-        if Delay_string == "":
-            return measurement
         N_delay = len(Delay_string) - 1
         delays = []
         for i in range(N_delay):
@@ -273,30 +271,21 @@ def CalculateNestedDelay(self, Delay_string, measurement, doppler=True):
                 total_delay_time * self.fsample,
                 order=self.delay_order,
             )
+        return total_delay_time
+
+    def CalculateNestedDelay(self, Delay_string, measurement, doppler=True):
+        """
+        calculate nested delay of any measurement
+        """
+        if Delay_string == "":
+            return measurement
+        total_delay_time = self.CalculateNestedDelayTime(Delay_string)
         if doppler:
             total_delay_doppler = np.gradient(total_delay_time, 1.0 / self.fsample)
             return (1.0 + total_delay_doppler) * timeshift(measurement, total_delay_time * self.fsample, order=self.order)
         else:
             return timeshift(measurement, total_delay_time * self.fsample, order=self.order)
 
-    def CalculateNestedDelayTime(self, Delay_string):
-        """
-        calculate the total delay time of a multiple delay operator
-        """
-        N_delay = len(Delay_string) - 1
-        delays = []
-        for i in range(N_delay):
-            delays.append(Delay_string[i] + Delay_string[i + 1])
-        total_delay_time = np.zeros(self.size)
-        for i in range(N_delay):
-            delay_idx = delays[i]
-            total_delay_time = total_delay_time - timeshift(
-                self.delays[delay_idx],
-                total_delay_time * self.fsample,
-                order=self.delay_order,
-            )
-        return total_delay_time
-
     def CalculateNestedDelayCombination(self, Delay_string_plus, Delay_string_minus, measurement, doppler=True):
         N_plus = len(Delay_string_plus)
         N_minus = len(Delay_string_minus)
@@ -976,6 +965,23 @@ def TDI_noise(self, f, P_ij=None, P_ij_strings=None):
             total_noise += noise_oms_ij[key] * PSDOMS + noise_acc_ij[key] * PSDACC  # (Nf)
         return total_noise
 
+    def TDI_average_response(self, f, P_ij=None, P_ij_strings=None, Nsource=1024):
+        """Calculate average response from Pij or Pij strings. 
+        The response function is averaged over N sources in 
+        random directions.
+        """
+        if P_ij is None:
+            P_ij = self.TDI_P_ij(P_ij_strings=P_ij_strings, f=f)
+
+        response_arr = []
+        for _ in tqdm(range(Nsource)):
+            longitude = np.random.uniform(0, TWOPI)
+            latitude = np.arcsin(np.random.uniform(-1, 1))
+            response_arr.append(self.TDI_response_function(lam=longitude, beta=latitude, f=f, P_ij=P_ij))
+        response_arr = np.array(response_arr)  # (Nsource, Nf)
+        response_avg = np.mean(response_arr, axis=0)  # (Nf)
+        return response_avg
+
     def TDI_sensitivity(self, f, P_ij=None, P_ij_strings=None, Nsource=1024):
         """
         calculate sensitivity from Pij or the Pij strings. the response funciton is averaged over N sources in random directions.
@@ -988,16 +994,9 @@ def TDI_sensitivity(self, f, P_ij=None, P_ij_strings=None, Nsource=1024):
         PSD = self.TDI_noise(f=f, P_ij=P_ij)  # (Nf)
 
         # calculate the average response function
-        Response_arr = []
-        for _ in tqdm(range(Nsource)):
-            longitude = np.random.uniform(0, TWOPI)
-            latitude = np.arcsin(np.random.uniform(-1, 1))
-            Response_arr.append(self.TDI_response_function(lam=longitude, beta=latitude, f=f, P_ij=P_ij))
-        Response_arr = np.array(Response_arr)  # (Nsource, Nf)
-        Response_avg = np.mean(Response_arr, axis=0)  # (Nf)
-
+        response_avg = self.TDI_average_response(f, P_ij, Nsource=Nsource)
         # return sensitivity
-        return PSD / Response_avg
+        return PSD / response_avg
 
     def TDI_noise_CSD(self, f, P_ij=None, Q_ij=None, P_ij_strings=None, Q_ij_strings=None, return_PSD=False):
         """calculate the noise CSD of channel P and Q, i.e. 2/T <P^* Q>, reduce to PSD if P = Q"""