diff --git a/seispy/ccpprofile.py b/seispy/ccpprofile.py
index 6a66b426..aff61030 100644
--- a/seispy/ccpprofile.py
+++ b/seispy/ccpprofile.py
@@ -24,6 +24,7 @@ def _prof_range(lat, lon):
 
 def create_center_bin_profile(stations, val=5, method='linear'):
     """Create bins along a profile with given stations
+
     :param stations: Stations along a profile
     :type stations: :class:`seispy.rf2depth_makedata.Station`
     :param val: The interval between two points in km
@@ -109,6 +110,7 @@ def init_profile(lat1, lon1, lat2, lon2, val):
 class CCPProfile():
     def __init__(self, cfg_file=None, log:SetupLog=SetupLog()):
         """CCPProfile class for CCP stacking along a profile
+
         :param cfg_file: Configure file for CCP stacking
         :type cfg_file: str
         :param log: Log class for logging
@@ -118,12 +120,12 @@ def __init__(self, cfg_file=None, log:SetupLog=SetupLog()):
         if cfg_file is None:
             self.cpara = CCPPara()
         elif isinstance(cfg_file, str):
-            self.load_para(cfg_file)
+            self._load_para(cfg_file)
         else:
             raise ValueError('cfg_file must be str format.')
         self.stack_data = []
 
-    def load_para(self, cfg_file):
+    def _load_para(self, cfg_file):
         try:
             self.cpara = ccppara(cfg_file)
         except Exception as e:
@@ -277,12 +279,19 @@ def stack(self):
             self.stack_data.append(boot_stack)   
 
     def save_stack_data(self, format='npz'):
-        """If format is \'npz\', saving stacked data and parameters to local as a npz file. To load the file, please use data = np.load(fname, allow_pickle=True).
-        data['cpara'] is the parameters when CCP stacking.
-        data['stack_data'] is the result of stacked data.
+        """If format is ``npz``, saving stacked data and parameters to local as a npz file. To load the file, please use data = np.load(fname, allow_pickle=True).
+        ``data['cpara']`` is the parameters when CCP stacking.
+        ``data['stack_data']`` is the result of stacked data.
 
-        If format is \'dat\' the stacked data will be save into a txt file with 8 columns, including bin_lat, bin_lon, profile_dis, depth, amp, ci_low, ci_high and count.
-        where bin_lat and bin_lon represent the position of each bin; profile_dis represents the distance in km between each bin and the start point of the profile; depth represents depth of each bin; amp means the stacked amplitude; ci_low and ci_high mean confidence interval with bootstrap method; count represents stacking number of each bin.
+        If format is ``dat`` the stacked data will be save into a txt file with 8 columns,
+        including ``bin_lat``, ``bin_lon``, ``profile_dis``, ``depth``, ``amp``, ``ci_low``, ``ci_high`` and ``count``.
+        
+        - ``bin_lat`` and ``bin_lon`` represent the position of each bin; 
+        - ``profile_dis`` represents the distance in km between each bin and the start point of the profile;
+        - ``depth`` represents depth of each bin;
+        - ``amp`` means the stacked amplitude; 
+        - ``ci_low`` and ``ci_high`` mean confidence interval with bootstrap method;
+        - ``count`` represents stacking number of each bin.
 
         :param format: Format for stacked data
         :type format: str
diff --git a/seispy/distaz.py b/seispy/distaz.py
index ae0479e2..2acb0db0 100644
--- a/seispy/distaz.py
+++ b/seispy/distaz.py
@@ -2,59 +2,15 @@
 import numpy as np
 
 
-def sind(deg):
-    rad = math.radians(deg)
-    return math.sin(rad)
-
-
-def cosd(deg):
-    rad = math.radians(deg)
-    return math.cos(rad)
-
-
-def tand(deg):
-    rad = math.radians(deg)
-    return math.tan(rad)
-
-
-def cotd(deg):
-    rad = math.radians(deg)
-    return math.cos(rad) / math.sin(rad)
-
-
-def asind(x):
-    rad = math.asin(x)
-    return math.degrees(rad)
-
-
-def acosd(x):
-    rad = math.acos(x)
-    return math.degrees(rad)
-
-
-def atand(x):
-    rad = math.atan(x)
-    return math.degrees(rad)
-
-
-def km2deg(kilometers):
-    return kilometers / 111.19
-
-
-def deg2km(degree):
-    return degree * 111.19
-
-
 class distaz:
     """
-    c Subroutine to calculate the Great Circle Arc distance
-    c    between two sets of geographic coordinates
-    c
-    c Equations take from Bullen, pages 154, 155
-    c
-    c T. Owens, September 19, 1991
-    c           Sept. 25 -- fixed az and baz calculations
-    c
+    Subroutine to calculate the Great Circle Arc distance
+        between two sets of geographic coordinates
+    
+    Equations take from Bullen, pages 154, 155
+    
+    T. Owens, September 19, 1991
+              Sept. 25 -- fixed az and baz calculations
     P. Crotwell, Setember 27, 1995
     Converted to c to fix annoying problem of fortran giving wrong
     answers if the input doesn't contain a decimal point.
diff --git a/seispy/eq.py b/seispy/eq.py
index 2043b78f..30b39be2 100644
--- a/seispy/eq.py
+++ b/seispy/eq.py
@@ -25,6 +25,8 @@ def __str__(self):
         return '{}'.format(self.matchkey)
 
 def rotateZNE(st):
+    """Rotate Z, N, E components to Z, N, E components
+    """
     try:
         zne = rotate2zne(
             st[0], st[0].stats.sac.cmpaz, st[0].stats.sac.cmpinc,
@@ -92,6 +94,7 @@ def _cleanstream(self):
     @classmethod
     def from_stream(cls, stream):
         """Create EQ object from obspy stream
+
         :param stream: obspy stream
         :type stream: obspy.Stream
         :return: EQ object
@@ -138,6 +141,13 @@ def channel_correct(self, switchEN=False, reverseE=False, reverseN=False):
             self.st.select(channel='*[N1]')[0].stats.channel = chE
 
     def write(self, path, evt_datetime):
+        """Write raw stream to SAC files
+
+        :param path: path to save SAC files
+        :type path: str
+        :param evt_datetime: event datetime
+        :type evt_datetime: obspy.core.utcdatetime.UTCDateTime
+        """
         for tr in self.st:
             sac = SACTrace.from_obspy_trace(tr)
             sac.b = 0
@@ -158,6 +168,7 @@ def detrend(self):
 
     def filter(self, freqmin=0.05, freqmax=1, order=4):
         """Bandpass filter
+
         :param freqmin: minimum frequency, defaults to 0.05
         :type freqmin: float, optional
         :param freqmax: maximum frequency, defaults to 1
@@ -169,6 +180,7 @@ def filter(self, freqmin=0.05, freqmax=1, order=4):
 
     def get_arrival(self, model, evdp, dis, phase='P'):
         """Get arrival time, ray parameter and incident angle from TauP model
+
         :param model: TauP model
         :type model: TauPyModel
         :param evdp: focal depth
@@ -192,6 +204,7 @@ def get_arrival(self, model, evdp, dis, phase='P'):
 
     def search_inc(self, bazi):
         """Search incident angle for S wave
+
         :param bazi: back azimuth
         :type bazi: float
         :return: incident angle
@@ -211,6 +224,7 @@ def search_inc(self, bazi):
 
     def search_baz(self, bazi, time_b=10, time_e=20, offset=90):
         """Search back azimuth for P wave
+
         :param bazi: back azimuth
         :type bazi: float
         :param time_b: time before P arrival, defaults to 10
@@ -261,6 +275,7 @@ def fix_channel_name(self):
 
     def rotate(self, baz, inc=None, method='NE->RT', search_inc=False, baz_shift=0):
         """Rotate to radial and transverse components
+
         :param baz: back azimuth
         :type baz: float
         :param inc: incident angle, defaults to None
@@ -294,6 +309,7 @@ def rotate(self, baz, inc=None, method='NE->RT', search_inc=False, baz_shift=0):
 
     def snr(self, length=50):
         """Calculate SNR
+
         :param length: length for noise, defaults to 50
         :type length: int, optional
         :return: SNR of E, N, Z components
@@ -317,6 +333,7 @@ def snr(self, length=50):
     
     def get_time_offset(self, event_time=None):
         """Get time offset from SAC header
+
         :param event_time: event time, defaults to None
         :type event_time: obspy.core.utcdatetime.UTCDateTime, optional
         """
@@ -347,6 +364,7 @@ def _get_time(self, time_before, time_after):
 
     def phase_trigger(self, time_before, time_after, prepick=True, stl=5, ltl=10):
         """ Trigger P or S phase
+
         :param time_before: time before P or S arrival
         :type time_before: float
         :param time_after: time after P or S arrival
@@ -445,6 +463,7 @@ def deconvolute(self, shift, time_after, f0=2.0, method='iter', only_r=False,
     
     def decon_p(self, tshift, tcomp=False, **kwargs):
         """Deconvolution for P wave
+
         :param tshift: Time shift before P arrival
         :type tshift: float
         :param tcomp: Whether calculate transverse component, defaults to False
@@ -468,6 +487,7 @@ def decon_p(self, tshift, tcomp=False, **kwargs):
 
     def decon_s(self, tshift, **kwargs):
         """Deconvolution for S wave
+
         :param tshift: Time shift before P arrival
         :type tshift: float
         """
@@ -487,6 +507,7 @@ def saverf(self, path, evtstr=None, shift=0, evla=-12345.,
                evlo=-12345., evdp=-12345., mag=-12345.,
                gauss=0, baz=-12345., gcarc=-12345., only_r=False, **kwargs):
         """Save receiver function to SAC file
+
         :param path: path to save SAC file
         :type path: str
         :param evtstr: event string, defaults to None
@@ -560,6 +581,7 @@ def _s_condition(self, trrf, shift):
 
     def judge_rf(self, gauss, shift, npts, criterion='crust', rmsgate=None):
         """Judge whether receiver function is valid
+        
         :param gauss: Gaussian factor
         :type gauss: float
         :param shift: time shift before P arrival
diff --git a/seispy/hk.py b/seispy/hk.py
index 80f69767..f6aba96f 100644
--- a/seispy/hk.py
+++ b/seispy/hk.py
@@ -136,6 +136,7 @@ def ci(allstack, h, kappa, ev_num):
     """
     Search best H and kappa from stacked matrix.
     Calculate error for H and kappa
+
     :param allstack: stacked HK matrix
     :param h: 1-D array of H
     :param kappa: 1-D array of kappa
diff --git a/seispy/rf2depth_makedata.py b/seispy/rf2depth_makedata.py
index 5e31fb6b..209534cb 100644
--- a/seispy/rf2depth_makedata.py
+++ b/seispy/rf2depth_makedata.py
@@ -74,6 +74,8 @@ class RFDepth():
     def __init__(self, cpara:CCPPara, log:Logger=SetupLog().RF2depthlog,
                  raytracing3d=False, velmod3d=None, modfolder1d=None) -> None:
         """
+        Convert receiver function to depth axis
+        
         :param cpara: CCPPara object
         :type cpara: CCPPara
         :param log: Log object
diff --git a/seispy/rfani.py b/seispy/rfani.py
index f4cbe19c..a8fd6af0 100644
--- a/seispy/rfani.py
+++ b/seispy/rfani.py
@@ -7,14 +7,14 @@
 from seispy.utils import load_cyan_map
 
 
-def joint_stack(energy_r, energy_cc, energy_tc, weight=[0.4, 0.3, 0.3]):
+def _joint_stack(energy_r, energy_cc, energy_tc, weight=[0.4, 0.3, 0.3]):
     energy_r = energy_r / np.max(energy_r)
     energy_cc = energy_cc / np.max(energy_cc)
     energy_tc = energy_tc / np.max(energy_tc)
     return np.exp(np.log(energy_r) * weight[0] + np.log(energy_cc) * weight[1] - np.log(energy_tc) * weight[2])
 
 
-def average_delay(fd, td):
+def _average_delay(fd, td):
     uniq_fd = np.unique(fd)
     if uniq_fd.size > 2:
         raise ValueError('FVD do not converge, {} gotten'.format(uniq_fd))
@@ -55,19 +55,21 @@ def __init__(self, sacdatar, tb, te, tlen=3, val=10, rayp=0.06, model='iasp91'):
         self.nbs = int((self.tb + self.sacdatar.shift) / self.sacdatar.sampling)
         self.nes = int((self.te + self.sacdatar.shift) / self.sacdatar.sampling)
         self.sacdatar.moveoutcorrect(ref_rayp=rayp, velmod=model, replace=True)
-        self.baz_stack(val=val)
-        self.search_peak_amp()
-        self.init_ani_para()
+        self._baz_stack(val=val)
+        self._search_peak_amp()
+        self._init_ani_para()
         self.fvd, self.deltat = np.meshgrid(self.fvd_1d, self.deltat_1d)
 
-    def baz_stack(self, val=10):
+    def _baz_stack(self, val=10):
+        """Stack RF data with back-azimuth.
+        """
         self.stack_range = np.arange(0, 360, val)
         stacked_data = self.sacdatar.bin_stack(lim=[0, 360], val=val)
         self.rfr_baz = stacked_data['data_prime']
         self.rft_baz = stacked_data['datat']
         self.count_baz = stacked_data['count']
 
-    def search_peak_amp(self):
+    def _search_peak_amp(self):
         mean_rf = np.mean(self.rfr_baz, axis=0)
         nmax = extrema(mean_rf[self.nbs:self.nes])+self.nbs
         if nmax.size > 1:
@@ -77,11 +79,13 @@ def search_peak_amp(self):
         self.nb = int(nps - self.tlen / self.sacdatar.sampling)
         self.ne = int(nps + self.tlen / self.sacdatar.sampling)
 
-    def init_ani_para(self):
+    def _init_ani_para(self):
         self.deltat_1d = np.arange(0, 1.55, 0.05)
         self.fvd_1d = np.arange(0, 365, 5)
 
-    def cut_energy_waveform(self, idx, nb, ne):
+    def _cut_energy_waveform(self, idx, nb, ne):
+        """Cut energy waveform with given index, nb and ne.
+        """
         engr = np.zeros([nb.shape[0], nb.shape[1], self.ne-self.nb])
         for i in range(nb.shape[0]):
             for j in range(nb.shape[1]):
@@ -89,6 +93,8 @@ def cut_energy_waveform(self, idx, nb, ne):
         return engr
 
     def radial_energy_max(self):
+        """Calculate the energy of radial component.
+        """
         energy = np.zeros([self.fvd.shape[0], self.fvd.shape[1], self.ne-self.nb])
         # tmp_data = np.zeros(self.ne-self.nb)
         for i, baz in enumerate(self.stack_range):
@@ -96,16 +102,20 @@ def radial_energy_max(self):
             nt_corr = (t_corr / self.sacdatar.sampling).astype(int)
             new_nb = self.nb - nt_corr
             new_ne = self.ne - nt_corr
-            energy += self.cut_energy_waveform(i, new_nb, new_ne)
+            energy += self._cut_energy_waveform(i, new_nb, new_ne)
         energy = np.max(energy ** 2, axis=2)
         energy /= np.max(np.sum(self.rfr_baz[:, self.nb:self.ne], axis=0)**2)
         return energy
 
     def xyz2grd(self, energy):
+        """Interpolate energy to grid.
+        """
         self.fvd, self.deltat = np.meshgrid(self.fvd_1d, self.deltat_1d)
         return griddata(self.ani_points, energy, (self.fvd, self.deltat))
 
     def rotate_to_fast_slow(self):
+        """Rotate RF data to fast and slow direction.
+        """
         self.ani_points = np.empty([0, 2])
         for f in self.fvd_1d:
             for d in self.deltat_1d:
@@ -137,6 +147,15 @@ def rotate_to_fast_slow(self):
         return self.xyz2grd(energy_cc), self.xyz2grd(energy_tc)
 
     def plot_stack_baz(self, enf=60, outpath='./'):
+        """Plot the stack of RF data with back-azimuth.
+        
+        Parameters
+        ----------
+        enf : int, optional
+            Enlarge factor for back-azimuth, by default 60
+        outpath : str, optional
+            Output path for saving the figure, by default './'
+        """
         ml = MultipleLocator(5)
         bound = np.zeros_like(self.sacdatar.time_axis)
         plt.style.use("bmh")
@@ -191,6 +210,19 @@ def plot_stack_baz(self, enf=60, outpath='./'):
         fig.savefig(join(outpath, '{}_baz_stack.png'.format(self.sacdatar.staname)), dpi=400, bbox_inches='tight')
 
     def plot_correct(self, fvd=0, dt=0.44, enf=80, outpath=None):
+        """Plot the RF data with back-azimuth and time after P corrected.
+
+        Parameters
+        ----------
+        fvd : int, optional
+            Fast velocity direction, by default 0
+        dt : float, optional
+            Time delay for correction, by default 0.44
+        enf : int, optional
+            Enlarge factor for back-azimuth, by default 80
+        outpath : str, optional
+            Output path for saving the figure, by default None
+        """
         nt_corr = int((dt/2 / self.sacdatar.sampling))
         # nt_fast = np.arange(self.nb, self.ne) + nt_corr
         # nt_slow = np.arange(self.nb, self.ne) - nt_corr
@@ -230,6 +262,15 @@ def plot_correct(self, fvd=0, dt=0.44, enf=80, outpath=None):
             plt.savefig(join(outpath, 'rf_corrected.png'), dpi=400, bbox_inches='tight')
 
     def search_peak_list(self, energy, opt='max'):
+        """Search the peak of energy.
+
+        Parameters
+        ----------
+        energy : np.ndarray
+            Energy matrix
+        opt : str, optional
+            Option for searching peak, by default 'max'
+        """
         if opt == 'max':
             ind = np.argwhere(energy == np.max(energy))
         elif opt == 'min':
@@ -239,6 +280,15 @@ def search_peak_list(self, energy, opt='max'):
         return self.ani_points[ind][:, 0], self.ani_points[ind][:, 1]
 
     def search_peak(self, energy, opt='max'):
+        """Search the peak of energy.
+
+        Parameters
+        ----------
+        energy : np.ndarray
+            Energy matrix
+        opt : str, optional
+            Option for searching peak, by default 'max'
+        """
         if opt == 'max':
             ind = np.argwhere(energy == np.max(energy))
         elif opt == 'min':
@@ -250,13 +300,20 @@ def search_peak(self, energy, opt='max'):
         for i, j in ind:
             best_fvd.append(self.fvd[i, j])
             best_dt.append(self.deltat[i, j])
-        uniq_fd, uniq_td = average_delay(np.array(best_fvd), np.array(best_dt))
+        uniq_fd, uniq_td = _average_delay(np.array(best_fvd), np.array(best_dt))
         return uniq_fd, uniq_td
 
     def joint_ani(self, weight=[0.5, 0.3, 0.2]):
+        """Joint method for crustal anisotropy estimation.
+
+        Parameters
+        ----------
+        weight : list, optional
+            Weight for three energy matrix, by default [0.5, 0.3, 0.2]
+        """
         self.energy_r = self.radial_energy_max()
         self.energy_cc, self.energy_tc = self.rotate_to_fast_slow()
-        self.energy_joint = joint_stack(self.energy_r, self.energy_cc, self.energy_tc, weight)
+        self.energy_joint = _joint_stack(self.energy_r, self.energy_cc, self.energy_tc, weight)
         self.bf, self.bt = self.search_peak(self.energy_joint, opt='max')
         return self.bf, self.bt
 
diff --git a/seispy/rfcorrect.py b/seispy/rfcorrect.py
index 308bf533..3d3b75ea 100644
--- a/seispy/rfcorrect.py
+++ b/seispy/rfcorrect.py
@@ -55,7 +55,7 @@ def __init__(self, data_path, only_r=False, prime_comp='R'):
             np.loadtxt(evt_lst, dtype=self.dtype, unpack=True, ndmin=1)
         self.rayp = skm2srad(self.rayp)
         self.ev_num = self.evla.shape[0]
-        self.read_sample(data_path)
+        self._read_sample(data_path)
         self.data_prime = np.empty([self.ev_num, self.rflength])
         # self.__dict__['data{}'.format(self.comp.lower())] = np.empty([self.ev_num, self.rflength])
         # eval('self.data_prime = self.data{}'.format(self.comp.lower()))
@@ -72,7 +72,12 @@ def __init__(self, data_path, only_r=False, prime_comp='R'):
                 self.data_prime[_i] = sac.data
         exec('self.data{} = self.data_prime'.format(self.comp.lower()))
 
-    def read_sample(self, data_path):
+    def _read_sample(self, data_path):
+        """Read sample SAC file to get station information
+
+        :param data_path: Path to RF data with SAC format
+        :type data_path: str
+        """
         fname = glob.glob(join(data_path, self.event[0] + '_' + self.phase[0] + '_{}.sac'.format(self.comp)))
         if len(fname) == 0:
             raise FileNotFoundError('No such files with comp of {} in {}'.format(self.comp, data_path))
@@ -90,7 +95,7 @@ def read_sample(self, data_path):
         self.sampling = sample_sac.delta
         self.time_axis = np.arange(self.rflength) * self.sampling - self.shift
 
-    def init_property(self, ev_num):
+    def _init_property(self, ev_num):
         self.ev_num = ev_num
         self.event = np.array(['']*ev_num)
         self.phase = np.array(['']*ev_num)
@@ -145,7 +150,7 @@ def read_stream(cls, stream, rayp, baz, prime_comp='R', stream_t=None):
             baz = np.ones(ev_num)*baz
         if ev_num != rayp.size or ev_num != baz.size:
             raise ValueError('Array length of rayp and baz must be the same as stream')
-        rfsta.init_property(ev_num)
+        rfsta._init_property(ev_num)
         try:
             rfsta.staname = '{}.{}'.format(stream[0].stats.sac.knetwk, stream[0].stats.sac.kstnm)
             rfsta.stla = stream[0].stats.sac.stla
@@ -226,6 +231,7 @@ def _chech_comp(self):
 
     def normalize(self, method='single'):
         """Normalize amplitude of each RFs.
+
         :param method: Method of normalization with ``single`` and ``average`` avaliable.
                      - ``single`` for normalization with max amplitude of current RF.
                      - ``average`` for normalization with average amplitude of current station.
@@ -249,7 +255,6 @@ def normalize(self, method='single'):
     def resample(self, dt):
         """Resample RFs with specified dt
 
-
         Parameters
         ----------
         dt : ``float``
@@ -353,16 +358,51 @@ def psrf_1D_raytracing(self, dep_range=np.arange(0, 150), **kwargs):
         return pplat_s, pplon_s, tps
 
     def psrf_3D_raytracing(self, mod3dpath, dep_range=np.arange(0, 150), srayp=None):
+        """3D back ray tracing to obtained Ps conversion points at discret depths
+
+        :param mod3dpath: Path to 3D velocity model
+        :type mod3dpath: str
+        :param dep_range: Discret conversion depth, defaults to np.arange(0, 150)
+        :type dep_range: numpy.ndarray, optional
+        :param srayp: Ray-parameter lib for Ps phases, If set up to None the rayp of direct is used, defaults to None
+        :type srayp: numpy.lib.npyio.NpzFile, optional
+        :return pplat_s: Latitude of conversion points
+        :return pplon_s: Longitude of conversion points
+        :return tps: Time difference of Ps at each depth
+        :rtype: list
+        """
         self.dep_range = dep_range
         mod3d = Mod3DPerturbation(mod3dpath, dep_range)
         pplat_s, pplon_s, _, _, tps = psrf_3D_raytracing(self, dep_range, mod3d, srayp=srayp)
         return pplat_s, pplon_s, tps
 
     def psrf_3D_moveoutcorrect(self, mod3dpath, **kwargs):
+        """3D moveout correction with 3D velocity model
+
+        :param mod3dpath: Path to 3D velocity model
+        :type mod3dpath: str
+        :param dep_range: Discret conversion depth, defaults to np.arange(0, 150)
+        :type dep_range: numpy.ndarray, optional
+        :param srayp: Ray-parameter lib for Ps phases, If set up to None the rayp of direct is used, defaults to None
+        :type srayp: numpy.lib.npyio.NpzFile, optional
+        :return: 2D array of RFs in depth
+        :rtype: numpy.ndarray
+        """
         warnings.warn('The fuction will be change to RFStation.psrf_3D_timecorrect in the future')
         self.psrf_3D_timecorrect(mod3dpath, **kwargs)
 
     def psrf_3D_timecorrect(self,  mod3dpath, dep_range=np.arange(0, 150), normalize='single', **kwargs):
+        """3D time-to-depth conversion with 3D velocity model
+        
+        :param mod3dpath: Path to 3D velocity model
+        :type mod3dpath: str
+        :param dep_range: Discret conversion depth, defaults to np.arange(0, 150)
+        :type dep_range: numpy.ndarray, optional
+        :param normalize: Normalization method, defaults to 'single', see RFStation.normalize for detail
+        :type normalize: str, optional
+        :return: 2D array of RFs in depth
+        :rtype: numpy.ndarray
+        """
         self.dep_range = dep_range
         mod3d = Mod3DPerturbation(mod3dpath, dep_range)
         pplat_s, pplon_s, pplat_p, pplon_p, raylength_s, raylength_p, tps = psrf_1D_raytracing(self, dep_range, **kwargs)
@@ -399,6 +439,15 @@ def jointani(self, tb, te, tlen=3., stack_baz_val=10, rayp=0.06,
         return best_f, best_t
 
     def slantstack(self, ref_dis=None, rayp_range=None, tau_range=None):
+        """Slant stack for receiver function
+
+        :param ref_dis: reference distance, by default None
+        :type ref_dis: int or float, optional
+        :param rayp_range: range of ray parameter, by default None
+        :type rayp_range: numpy.ndarray, optional
+        :param tau_range: range of tau, by default None
+        :type tau_range: numpy.ndarray, optional
+        """
         self.slant = SlantStack(self.data_prime, self.time_axis, self.dis)
         self.slant.stack(ref_dis, rayp_range, tau_range)
         return self.slant.stack_amp
diff --git a/seispy/slantstack.py b/seispy/slantstack.py
index f333b103..2841e4da 100644
--- a/seispy/slantstack.py
+++ b/seispy/slantstack.py
@@ -26,6 +26,15 @@ def __init__(self, seis, timeaxis, dis) -> None:
         self.syn_drayp = np.array([])
 
     def stack(self, ref_dis=None, rayp_range=None, tau_range=None):
+        """Slant stack for receiver function
+
+        :param ref_dis: reference distance, by default None
+        :type ref_dis: int or float, optional
+        :param rayp_range: range of ray parameter, by default None
+        :type rayp_range: numpy.ndarray, optional
+        :param tau_range: range of tau, by default None
+        :type tau_range: numpy.ndarray, optional
+        """
         if ref_dis is not None and scalar_instance(ref_dis):
             self.ref_dis = ref_dis
         elif ref_dis is None:
@@ -53,6 +62,15 @@ def stack(self, ref_dis=None, rayp_range=None, tau_range=None):
         self.stack_amp /= ev_num
 
     def syn_tps(self, phase_list, velmodel='iasp91', focal_dep=10):
+        """Calculate the theoretical tau and reference rayp for the given phase list
+
+        :param phase_list: phase list for theoretical arrival time
+        :type phase_list: list
+        :param velmodel: velocity model, by default 'iasp91'
+        :type velmodel: str, optional
+        :param focal_dep: focal depth, by default 10
+        :type focal_dep: int or float, optional
+        """
         model = TauPyModel(model=velmodel)
         phase_list.insert(0, 'P')
         arrs = model.get_travel_times(focal_dep, self.ref_dis, phase_list=phase_list)