From 20f0d7b9f34336c5cfb75ede4f21815694fea65d Mon Sep 17 00:00:00 2001
From: JouSF-1409 <81711503+JouSF-1409@users.noreply.github.com>
Date: Tue, 30 Jan 2024 23:24:50 +0800
Subject: [PATCH 1/7] =?UTF-8?q?=E4=B8=BArayplib=E6=B7=BB=E5=8A=A0=E4=BA=86?=
 =?UTF-8?q?=E5=91=BD=E4=BB=A4=E8=A1=8C=E5=B7=A5=E5=85=B7=EF=BC=8C=20rf2dep?=
 =?UTF-8?q?th=E7=9A=84=E4=B8=80=E7=82=B9=E4=BC=98=E5=8C=96=EF=BC=8C?=
 =?UTF-8?q?=E4=BD=BF=E7=94=A8namedtuple=E6=94=B9=E5=96=84=E4=BA=86Station?=
 =?UTF-8?q?=E7=B1=BB=E7=9A=84=E4=BD=BF=E7=94=A8=20depmodel.py=20=E4=B8=AD?=
 =?UTF-8?q?=E7=9A=84ccp=5Fmodel=E5=A4=AA=E7=B9=81=E7=90=90=EF=BC=8C?=
 =?UTF-8?q?=E6=9C=AA=E6=9D=A5=E4=BC=9A=E6=8B=86=E5=88=86=E9=87=8D=E5=86=99?=
 =?UTF-8?q?=E3=80=82?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 seispy/ccp/rayplib.py       |  95 ++++++++++++++++++++++
 seispy/rf2depth_makedata.py | 157 ++++++++++++++++++++++++++++--------
 2 files changed, 218 insertions(+), 34 deletions(-)
 create mode 100644 seispy/ccp/rayplib.py

diff --git a/seispy/ccp/rayplib.py b/seispy/ccp/rayplib.py
new file mode 100644
index 00000000..f8f40b16
--- /dev/null
+++ b/seispy/ccp/rayplib.py
@@ -0,0 +1,95 @@
+"""
+use obspy.taup instead of taup in java
+"""
+
+import numpy as np
+
+from obspy.taup import TauPyModel
+
+def rayp_input(phase_main:str,
+               source_dist_in_degree:str,
+               conver_layers:str,
+               source_depth_in_km:str,
+               out_put_path:str):
+    """
+    anaylize input informations
+    only allow P, SKS, S, ScS
+    """
+    if phase_main not in ['P', 'S', 'SKS', 'ScS']:
+        raise ValueError('not a valid phase')
+
+    try:
+        # set default step value for 2 input
+        source_dist_in_degree +='/1';source_dist_in_degree = [float(_i) for _i in source_dist_in_degree.split('/')[:4]]
+        conver_layers +='/1';conver_layers = [float(_i) for _i in conver_layers.split('/')[:4]]
+        source_depth_in_km +='/2';source_depth_in_km = [float(_i) for _i in source_depth_in_km.split('/')]
+        if conver_layers[0] < 0:
+            conver_layers[0] = 0
+        conver_layers = np.arange(conver_layers[0], conver_layers[1], conver_layers[2])
+        source_dist_in_degree = np.arange(source_dist_in_degree[0], source_dist_in_degree[1], source_dist_in_degree[2])
+        source_depth_in_km = np.arange(source_depth_in_km[0],source_depth_in_km[1], source_depth_in_km[2])
+    except:
+        raise ValueError('invalid inputs in dist and layers')
+    # here we limit layer_max to 30km
+    if conver_layers[-1] < 30:
+        raise ValueError('theres no need to do such a thing')
+
+
+    rayp = cal_taup('iasp91', phase_main,
+                    conver_layers, source_dist_in_degree, source_depth_in_km)
+
+    np.savez(out_put_path, dis = source_dist_in_degree, dep = source_depth_in_km, layers = conver_layers, rayp = rayp)
+
+
+
+
+def cal_taup(model:str, phase:str,
+             layer:np.ndarray, dist:np.ndarray, depth:np.ndarray):
+    """
+    cal rayp from obspy.taup
+    wont cal layers above 11km, change min_con
+
+    phase_list: list of phase names, like ["P30s", "P50s"]
+    layer : conversion layers
+    dist : numpy.ndarray, dist in deg
+    depth : numpy.ndarray, source_depth in km
+
+    """
+    min_con = 11
+    if phase[-1] == 'P':
+        conv = 's'
+    else:
+        conv = 'p'
+    ## init taup and matrix
+    rayp_lib = np.zeros((dist.shape[0], depth.shape[0], layer.shape[0]))
+    model = TauPyModel(model=model)
+    cal_layers_slice = np.where(layer > min_con); head = cal_layers_slice[0][0]
+    phase_list = ["{}{}{}".format(phase,_d, conv) for _d in layer[cal_layers_slice]]
+
+    ### main
+    for _i, _deg in enumerate(dist):
+        for _j, _dep in enumerate(depth):
+            #print("now at degree {}, dep{}".format(_deg, _dep))
+            for _k, _phase in enumerate(phase_list):
+                arrivals=model.get_travel_times(source_depth_in_km=_dep,
+                                            distance_in_degree=_deg,
+                                            phase_list=[_phase])
+                if len(arrivals) == 0:
+                    continue
+                rayp_lib[_i,_j,head+_k] = arrivals[0].ray_param_sec_degree
+    if head > 0:
+        for _layer in range(head):
+            rayp_lib[:,:,_layer] = rayp_lib[:,:,head]
+    return rayp_lib
+
+
+if __name__ == "__main__":
+    # paras
+    phase = "S"
+    dist = "30/90/5"
+    layers = "0/100/2"
+    depth = '0/100/5'
+    out_put = "G:/Ps_rayp.npz"
+
+    # run
+    rayp_input(phase, dist, layers, depth, out_put)
diff --git a/seispy/rf2depth_makedata.py b/seispy/rf2depth_makedata.py
index 2a45b4e5..e56b5ef0 100644
--- a/seispy/rf2depth_makedata.py
+++ b/seispy/rf2depth_makedata.py
@@ -1,4 +1,18 @@
-from seispy.core.depmodel import _load_mod
+"""
+class RF2depth : process cal RF2depth
+class sta_part, sta_full, _RFInd
+
+"""
+from os.path import join, exists
+import argparse
+import sys
+import glob
+from collections import namedtuple
+from logging import Logger
+
+import numpy as np
+
+from seispy.core.depmodel import _load_mod, DepModel
 from seispy.rfcorrect import RFStation, psrf2depth, psrf_1D_raytracing,\
     psrf_3D_migration, time2depth, psrf_3D_raytracing
 from seispy.core.pertmod import Mod3DPerturbation
@@ -6,10 +20,30 @@
 from seispy.ccppara import ccppara, CCPPara
 from seispy.setuplog import SetupLog
 from seispy.geo import latlon_from, rad2deg
-from os.path import join, exists
-import argparse
-import sys
-import glob
+
+# sta_part is not used
+sta_part = namedtuple('sta_part', ['station', 'stla', 'stlo'])
+sta_full = namedtuple('sta_full', ['station', 'stla', 'stlo', 'stel'])
+_RfInDepth = namedtuple('RfInDepth',
+                       ['station', 'stalat', 'stalon',
+                        'depthrange', 'bazi', 'rayp',
+                        'moveout_correct', 'piercelat', 'piercelon', 'stopindex'])
+class RFInDepth(_RfInDepth):
+    """
+    this class is inheried from a namedtuple, and can be called as a dict
+    mainly to restrict behaviour while output  and
+    provide a better speed while
+
+    warnings: this class is still in develpment
+    """
+    def __getitem__(self, item):
+        """
+        make RFInDepth can be called like dict
+        >>> ll = RFInDepth(station='sta', stalat='stla', stalon='stlo', depthrange='dep', bazi='baz',rayp='rayp', moveout_correct=0, piercelat=0, piercelon=0 )
+        >>> ll['stalat'], ll.stalat
+        ('stla', 'stla')
+        """
+        return getattr(self, item)
 
 
 class Station(object):
@@ -28,20 +62,43 @@ def __init__(self, sta_lst:str):
             self.station, self.stla, self.stlo = np.loadtxt(sta_lst, dtype=dtype, unpack=True, ndmin=1)
             self.stel = np.zeros(self.stla.size)
         self.sta_num = self.stla.shape[0]
+    def __getitem__(self,index):
+        """
+        allows  [] do indexing
+        """
+        if hasattr(self,'stel'):
+            return sta_full(self.station[index],
+                            self.stla[index], self.stlo[index], self.stel[index])
+        else:
+            return sta_full(self.station[index],
+                            self.stla[index], self.stlo[index], 0)
+    def __iter__(self):
+        """
+        allow for...in to iter
+        """
+        for index in range(self.stla.shape[0]):
+            if hasattr(self, 'stel'):
+                yield sta_full(self.station[index],
+                               self.stla[index], self.stlo[index], self.stel[index])
+            else:
+                yield sta_full(self.station[index], self.stla[index], self.stlo[index], 0)
+    def __len__(self):
+        return self.station.shape[0]
+
 
 
 class RFDepth():
     """Convert receiver function to depth axis
     """
-    def __init__(self, cpara:CCPPara, log=SetupLog(), 
-                 raytracing3d=False, velmod3d=None, modfolder1d=None) -> None:
+    def __init__(self, cpara:CCPPara, log:Logger=SetupLog().RF2depthlog,
+                 RAYTRACING3D=False, velmod3d=None, modfolder1d=None) -> None:
         """
         :param cpara: CCPPara object
         :type cpara: CCPPara
         :param log: Log object
-        :type log: SetupLog
-        :param raytracing3d: If True, use 3D ray tracing to calculate the travel time
-        :type raytracing3d: bool
+        :type log: logging.Logger
+        :param RAYTRACING3D: If True, use 3D ray tracing to calculate the travel time
+        :type RAYTRACING3D: bool
         :param velmod3d: Path to 3D velocity model in npz file
         :type velmod3d: str
         :param modfolder1d: Folder path to 1D velocity model files with staname.vel as the file name
@@ -51,22 +108,22 @@ def __init__(self, cpara:CCPPara, log=SetupLog(),
         self.cpara = cpara
         self.modfolder1d = modfolder1d
         self.log = log
-        self.raytracing3d = raytracing3d
+        self.RAYTRACING3D = RAYTRACING3D
         if velmod3d is not None:
             if isinstance(velmod3d, str):
                 self.mod3d = Mod3DPerturbation(velmod3d, cpara.depth_axis, velmod=cpara.velmod)
             else:
-                log.RF2depthlog.error('Path to 3d velocity model should be in str')
+                log.error('Path to 3d velocity model should be in str')
                 sys.exit(1)
         elif modfolder1d is not None:
             if isinstance(modfolder1d, str):
                 if exists(modfolder1d):
                     self.ismod1d = True
                 else:
-                    log.RF2depthlog.error('No such folder of {}'.format(modfolder1d))
+                    log.error('No such folder of {}'.format(modfolder1d))
                     sys.exit(1)
             else:
-                log.RF2depthlog.error('Folder to 1d velocity model files should be in str')
+                log.error('Folder to 1d velocity model files should be in str')
                 sys.exit(1)
         else:
             self.ismod1d = True
@@ -74,6 +131,7 @@ def __init__(self, cpara:CCPPara, log=SetupLog(),
             self.srayp = np.load(cpara.rayp_lib)
         else:
             self.srayp = None
+
         self.sta_info = Station(cpara.stalist)
         self.rfdepth = []
         self._test_comp()
@@ -86,7 +144,7 @@ def _test_comp(self):
                 self.prime_comp = comp
                 break
         if not self.prime_comp:
-            self.log.RF2depthlog.error('No such any RF files in \'R\',' 
+            self.log.error('No such any RF files in \'R\',' 
                                       '\'Q\', \'L\', and \'Z\' components')
             sys.exit(1)
 
@@ -96,29 +154,38 @@ def makedata(self, psphase=1):
         :param psphase: 1 for Ps, 2 for PpPs, 3 for PpSs
         :type psphase: int
         """
-        for i in range(self.sta_info.stla.shape[0]):
-            rfpath = join(self.cpara.rfpath, self.sta_info.station[i])
+        for _i, _sta in enumerate(self.sta_info):
+            rfpath = join(self.cpara.rfpath, _sta.station)
             stadatar = RFStation(rfpath, only_r=True, prime_comp=self.prime_comp)
-            stadatar.stel = self.sta_info.stel[i]
-            stadatar.stla = self.sta_info.stla[i]
-            stadatar.stlo = self.sta_info.stlo[i]
+            stadatar.stel = _sta.stel
+            stadatar.stla = _sta.stla
+            stadatar.stlo = _sta.stlo
             if stadatar.prime_phase == 'P':
                 sphere = True
             else:
                 sphere = False
-            self.log.RF2depthlog.info('the {}th/{} station with {} events'.format(i + 1, self.sta_info.stla.shape[0], stadatar.ev_num))
+            self.log.info('the {}th/{} station with {} events'.format(_i + 1, self.sta_info.shape[0], stadatar.ev_num))
+
+
+            #### 1d model for each station
             if self.ismod1d:
                 if self.modfolder1d is not None:
-                    velmod = _load_mod(self.modfolder1d, self.sta_info.station[i])
+                    velmod = _load_mod(self.modfolder1d, _sta.station)
+
                 else:
                     velmod = self.cpara.velmod
+
                 ps_rfdepth, end_index, x_s, _ = psrf2depth(stadatar, self.cpara.depth_axis,
-                            velmod=velmod, srayp=self.cpara.rayp_lib, sphere=sphere, phase=psphase)
+                            velmod=velmod, srayp=self.cpara.rayp_lib,
+                           sphere=sphere, phase=psphase)
+
                 piercelat, piercelon = np.zeros_like(x_s, dtype=np.float64), np.zeros_like(x_s, dtype=np.float64)
+
                 for j in range(stadatar.ev_num):
                     piercelat[j], piercelon[j] = latlon_from(self.sta_info.stla[i], self.sta_info.stlo[i],
                                                             stadatar.bazi[j], rad2deg(x_s[j]))
             else:
+                ### 3d model interp
                 if self.raytracing3d:
                     pplat_s, pplon_s, pplat_p, pplon_p, newtpds = psrf_3D_raytracing(stadatar, self.cpara.depth_axis, self.mod3d, srayp=self.srayp, sphere=sphere)
                 else:
@@ -130,59 +197,81 @@ def makedata(self, psphase=1):
                     piercelat, piercelon = pplat_s, pplon_s
                 else:
                     piercelat, piercelon = pplat_p, pplon_p
+
                 ps_rfdepth, end_index = time2depth(stadatar, self.cpara.depth_axis, newtpds)
-            rfdep = self._write_rfdep(stadatar, ps_rfdepth, piercelat, piercelon, end_index)
-            self.rfdepth.append(rfdep)
+
+            rfdep = self._append_rfdep(self.cpara,stadatar, ps_rfdepth, piercelat, piercelon, end_index)
         np.save(self.cpara.depthdat, self.rfdepth)
 
-    def _write_rfdep(self, stadata, amp, pplat, pplon, end_index):
+    def _write_rfdep(self,cpara, stadata, amp, pplat, pplon, end_index):
         rfdep = {}
         rfdep['station'] = stadata.staname
         rfdep['stalat'] = stadata.stla
         rfdep['stalon'] = stadata.stlo
-        rfdep['depthrange'] = self.cpara.depth_axis
+        rfdep['depthrange'] = cpara.depth_axis
         rfdep['bazi'] = stadata.bazi
         rfdep['rayp'] = stadata.rayp
         rfdep['moveout_correct'] = amp
         rfdep['piercelat'] = pplat
         rfdep['piercelon'] = pplon
         rfdep['stopindex'] = end_index
-        return rfdep
+        self.rfdepth.append(rfdep)
+
 
 
 def rf2depth():
-    """Convert receiver function to depth axis
+    """
+    CLI for Convert receiver function to depth axis
+    There's  4 branch provided to do RF 2 depth conversion
+
+    1. only -d :do moveout correction
+    2. only -r : do raytracing but no moveout correction
+    3. -d and -r : do moveout correction and raytracing
+    4. -m : use {staname}.vel file for RF2depth conversion
+
+
     """
     parser = argparse.ArgumentParser(description="Convert Ps RF to depth axis")
-    parser.add_argument('-d', help='Path to 3d vel model in npz file for moveout correcting',
+    parser.add_argument('-d',
+                        help='Path to 3d vel model in npz file for moveout correcting',
                         metavar='3d_velmodel_path', type=str, default='')
-    parser.add_argument('-m', help='Folder path to 1d vel model files with staname.vel as the file name',
+    parser.add_argument('-m',
+                        help='Folder path to 1d vel model files with staname.vel as the file name',
                         metavar='1d_velmodel_folder', type=str, default='')
-    parser.add_argument('-r', help='Path to 3d vel model in npz file for 3D ray tracing',
+    parser.add_argument('-r',
+                        help='Path to 3d vel model in npz file for 3D ray tracing',
                         metavar='3d_velmodel_path', type=str, default='')
-    parser.add_argument('cfg_file', type=str, help='Path to configure file')
+    parser.add_argument('cfg_file',
+                        help='Path to configure file',
+                        metavar='ccp.cfg', type=str)
     arg = parser.parse_args()
     cpara = ccppara(arg.cfg_file)
     if arg.d != '' and arg.r != '':
+        #### print help issue
         raise ValueError('Specify only 1 argument in \'-d\' and \'-r\'')
     elif arg.d != '' and arg.r == '' and arg.m == '':
+        #### do 3d moveout correction but 1d rf2depth
         raytracing3d = False
         velmod3d = arg.d
         modfolder1d = None
     elif arg.d == '' and arg.r != '' and arg.m == '':
+        #### do 3d raytraying
         raytracing3d = True
         velmod3d = arg.d
         modfolder1d = None
     elif arg.d == '' and arg.r == '' and arg.m != '':
+        #### use multiple 1d velmod for time2depth convertion
         raytracing3d = False
         velmod3d = None
         modfolder1d = arg.m
     else:
+        ### all last, use default 1D vel model do time2depth conversion
         raytracing3d = False
         velmod3d = None
         modfolder1d = None
     rfd = RFDepth(
-        cpara, raytracing3d=raytracing3d,
+        cpara,
+        RAYTRACING3D=raytracing3d,
         velmod3d=velmod3d,
         modfolder1d=modfolder1d,
     )

From ec5237f60faf54fe5aabea05e38aa4c5ca377733 Mon Sep 17 00:00:00 2001
From: JouSF-1409 <81711503+JouSF-1409@users.noreply.github.com>
Date: Wed, 31 Jan 2024 08:02:57 +0800
Subject: [PATCH 2/7] =?UTF-8?q?=E4=B8=BArayplib=E6=B7=BB=E5=8A=A0=E4=BA=86?=
 =?UTF-8?q?=E5=91=BD=E4=BB=A4=E8=A1=8C=E5=B7=A5=E5=85=B7=EF=BC=8C=20rf2dep?=
 =?UTF-8?q?th=E7=9A=84=E4=B8=80=E7=82=B9=E4=BC=98=E5=8C=96=EF=BC=8C?=
 =?UTF-8?q?=E4=BD=BF=E7=94=A8namedtuple=E6=94=B9=E5=96=84=E4=BA=86Station?=
 =?UTF-8?q?=E7=B1=BB=E7=9A=84=E4=BD=BF=E7=94=A8=20depmodel.py=20=E4=B8=AD?=
 =?UTF-8?q?=E7=9A=84ccp=5Fmodel=E5=A4=AA=E7=B9=81=E7=90=90=EF=BC=8C?=
 =?UTF-8?q?=E6=9C=AA=E6=9D=A5=E4=BC=9A=E6=8B=86=E5=88=86=E9=87=8D=E5=86=99?=
 =?UTF-8?q?=E3=80=82?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 seispy/ccp/rf2depth.py | 8 ++++++++
 1 file changed, 8 insertions(+)
 create mode 100644 seispy/ccp/rf2depth.py

diff --git a/seispy/ccp/rf2depth.py b/seispy/ccp/rf2depth.py
new file mode 100644
index 00000000..0a3d0a08
--- /dev/null
+++ b/seispy/ccp/rf2depth.py
@@ -0,0 +1,8 @@
+from collections import namedtuple
+
+import numpy as np
+
+
+if __name__ == '__main__':
+    import doctest
+    doctree = doctest.testmod()
\ No newline at end of file

From 63d86fcd540c080b7bbd240c53a5f23df0a8ddb9 Mon Sep 17 00:00:00 2001
From: JouSF-1409 <81711503+JouSF-1409@users.noreply.github.com>
Date: Wed, 31 Jan 2024 08:03:56 +0800
Subject: [PATCH 3/7] Update rf2depth_makedata.py

---
 seispy/rf2depth_makedata.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seispy/rf2depth_makedata.py b/seispy/rf2depth_makedata.py
index e56b5ef0..ff8296d9 100644
--- a/seispy/rf2depth_makedata.py
+++ b/seispy/rf2depth_makedata.py
@@ -164,7 +164,7 @@ def makedata(self, psphase=1):
                 sphere = True
             else:
                 sphere = False
-            self.log.info('the {}th/{} station with {} events'.format(_i + 1, self.sta_info.shape[0], stadatar.ev_num))
+            self.log.info('the {}th/{} station with {} events'.format(_i + 1, len(self.sta_info), stadatar.ev_num))
 
 
             #### 1d model for each station

From e315ee75a22f64ba93869284dfb1adaebf6a9111 Mon Sep 17 00:00:00 2001
From: JouSF-1409 <81711503+JouSF-1409@users.noreply.github.com>
Date: Wed, 31 Jan 2024 08:36:46 +0800
Subject: [PATCH 4/7] =?UTF-8?q?=E4=B8=BArayplib=E6=B7=BB=E5=8A=A0=E4=BA=86?=
 =?UTF-8?q?=E5=91=BD=E4=BB=A4=E8=A1=8C=E5=B7=A5=E5=85=B7=EF=BC=8C=20rf2dep?=
 =?UTF-8?q?th=E7=9A=84=E4=B8=80=E7=82=B9=E4=BC=98=E5=8C=96=EF=BC=8C?=
 =?UTF-8?q?=E4=BD=BF=E7=94=A8namedtuple=E6=94=B9=E5=96=84=E4=BA=86Station?=
 =?UTF-8?q?=E7=B1=BB=E7=9A=84=E4=BD=BF=E7=94=A8=20depmodel.py=20=E4=B8=AD?=
 =?UTF-8?q?=E7=9A=84ccp=5Fmodel=E5=A4=AA=E7=B9=81=E7=90=90=EF=BC=8C?=
 =?UTF-8?q?=E6=9C=AA=E6=9D=A5=E4=BC=9A=E6=8B=86=E5=88=86=E9=87=8D=E5=86=99?=
 =?UTF-8?q?=E3=80=82?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 seispy/rf2depth_makedata.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seispy/rf2depth_makedata.py b/seispy/rf2depth_makedata.py
index ff8296d9..94c2fea6 100644
--- a/seispy/rf2depth_makedata.py
+++ b/seispy/rf2depth_makedata.py
@@ -182,7 +182,7 @@ def makedata(self, psphase=1):
                 piercelat, piercelon = np.zeros_like(x_s, dtype=np.float64), np.zeros_like(x_s, dtype=np.float64)
 
                 for j in range(stadatar.ev_num):
-                    piercelat[j], piercelon[j] = latlon_from(self.sta_info.stla[i], self.sta_info.stlo[i],
+                    piercelat[j], piercelon[j] = latlon_from(_sta.stla, _sta.stlo,
                                                             stadatar.bazi[j], rad2deg(x_s[j]))
             else:
                 ### 3d model interp

From 8043706b273af038ee202e80f88f9331b857d76b Mon Sep 17 00:00:00 2001
From: JouSF-1409 <81711503+JouSF-1409@users.noreply.github.com>
Date: Wed, 31 Jan 2024 08:46:34 +0800
Subject: [PATCH 5/7] =?UTF-8?q?=E4=B8=BArayplib=E6=B7=BB=E5=8A=A0=E4=BA=86?=
 =?UTF-8?q?=E5=91=BD=E4=BB=A4=E8=A1=8C=E5=B7=A5=E5=85=B7=EF=BC=8C=20rf2dep?=
 =?UTF-8?q?th=E7=9A=84=E4=B8=80=E7=82=B9=E4=BC=98=E5=8C=96=EF=BC=8C?=
 =?UTF-8?q?=E4=BD=BF=E7=94=A8namedtuple=E6=94=B9=E5=96=84=E4=BA=86Station?=
 =?UTF-8?q?=E7=B1=BB=E7=9A=84=E4=BD=BF=E7=94=A8=20depmodel.py=20=E4=B8=AD?=
 =?UTF-8?q?=E7=9A=84ccp=5Fmodel=E5=A4=AA=E7=B9=81=E7=90=90=EF=BC=8C?=
 =?UTF-8?q?=E6=9C=AA=E6=9D=A5=E4=BC=9A=E6=8B=86=E5=88=86=E9=87=8D=E5=86=99?=
 =?UTF-8?q?=E3=80=82?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 seispy/rf2depth_makedata.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seispy/rf2depth_makedata.py b/seispy/rf2depth_makedata.py
index 94c2fea6..87d7f92f 100644
--- a/seispy/rf2depth_makedata.py
+++ b/seispy/rf2depth_makedata.py
@@ -200,7 +200,7 @@ def makedata(self, psphase=1):
 
                 ps_rfdepth, end_index = time2depth(stadatar, self.cpara.depth_axis, newtpds)
 
-            rfdep = self._append_rfdep(self.cpara,stadatar, ps_rfdepth, piercelat, piercelon, end_index)
+            rfdep = self._write_rfdep(self.cpara,stadatar, ps_rfdepth, piercelat, piercelon, end_index)
         np.save(self.cpara.depthdat, self.rfdepth)
 
     def _write_rfdep(self,cpara, stadata, amp, pplat, pplon, end_index):

From 262727a9dd67ddf179732d53cfe5a8a9483233ec Mon Sep 17 00:00:00 2001
From: JouSF-1409 <81711503+JouSF-1409@users.noreply.github.com>
Date: Wed, 31 Jan 2024 08:46:34 +0800
Subject: [PATCH 6/7] =?UTF-8?q?=E4=B8=BArayplib=E6=B7=BB=E5=8A=A0=E4=BA=86?=
 =?UTF-8?q?=E5=91=BD=E4=BB=A4=E8=A1=8C=E5=B7=A5=E5=85=B7=EF=BC=8C=20rf2dep?=
 =?UTF-8?q?th=E7=9A=84=E4=B8=80=E7=82=B9=E4=BC=98=E5=8C=96=EF=BC=8C?=
 =?UTF-8?q?=E4=BD=BF=E7=94=A8namedtuple=E6=94=B9=E5=96=84=E4=BA=86Station?=
 =?UTF-8?q?=E7=B1=BB=E7=9A=84=E4=BD=BF=E7=94=A8=20depmodel.py=20=E4=B8=AD?=
 =?UTF-8?q?=E7=9A=84ccp=5Fmodel=E5=A4=AA=E7=B9=81=E7=90=90=EF=BC=8C?=
 =?UTF-8?q?=E6=9C=AA=E6=9D=A5=E4=BC=9A=E6=8B=86=E5=88=86=E9=87=8D=E5=86=99?=
 =?UTF-8?q?=E3=80=82?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 seispy/rf2depth_makedata.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seispy/rf2depth_makedata.py b/seispy/rf2depth_makedata.py
index e56b5ef0..1685ce78 100644
--- a/seispy/rf2depth_makedata.py
+++ b/seispy/rf2depth_makedata.py
@@ -200,7 +200,7 @@ def makedata(self, psphase=1):
 
                 ps_rfdepth, end_index = time2depth(stadatar, self.cpara.depth_axis, newtpds)
 
-            rfdep = self._append_rfdep(self.cpara,stadatar, ps_rfdepth, piercelat, piercelon, end_index)
+            rfdep = self._write_rfdep(self.cpara,stadatar, ps_rfdepth, piercelat, piercelon, end_index)
         np.save(self.cpara.depthdat, self.rfdepth)
 
     def _write_rfdep(self,cpara, stadata, amp, pplat, pplon, end_index):

From 5d9baaedb5b4c06be4340619bd2f232177661721 Mon Sep 17 00:00:00 2001
From: Mijian Xu <gomijianxu@gmail.com>
Date: Thu, 1 Feb 2024 00:57:34 +0800
Subject: [PATCH 7/7] add more comments

---
 seispy/catalog.py           |  26 ++++++
 seispy/ccp/rf2depth.py      |   8 --
 seispy/ccppara.py           |   6 ++
 seispy/ccpprofile.py        |  40 +++++---
 seispy/decon.py             |  14 +++
 seispy/distaz.py            |   2 +-
 seispy/eq.py                | 149 +++++++++++++++++++++++++++--
 seispy/geo.py               | 180 ++++++++++++++++++++++++++++++++----
 seispy/harmonics.py         |   2 +
 seispy/para.py              |   6 ++
 seispy/plotR.py             |   6 ++
 seispy/rf2depth_makedata.py |  46 +++------
 12 files changed, 404 insertions(+), 81 deletions(-)
 delete mode 100644 seispy/ccp/rf2depth.py

diff --git a/seispy/catalog.py b/seispy/catalog.py
index 800328ea..0f968e43 100644
--- a/seispy/catalog.py
+++ b/seispy/catalog.py
@@ -5,12 +5,32 @@
 import pandas as pd
 
 def download_catalog(fname, server='IRIS', format='seispy', **kwargs):
+    """
+    Download catalog to local file
+    :param fname: file name
+    :type fname: str
+    :param server: service provider, defaults to IRIS
+    :type server: str
+    :param format: file format, defaults to seispy
+    :type format: str
+    :param kwargs: parameters for Query.get_events
+    :type kwargs: dict
+    """
     query = Query(server=server)
     query.get_events(**kwargs)
     write_catalog(query, fname, format)
 
 
 def write_catalog(query, fname, format):
+    """
+    Write catalog to file
+    :param query: Query object
+    :type query: seispy.io.Query
+    :param fname: file name
+    :type fname: str
+    :param format: file format
+    :type format: str
+    """
     if format == 'seispy':
         with open(fname, 'w') as f:
             for i, row in query.events.iterrows():
@@ -23,6 +43,12 @@ def write_catalog(query, fname, format):
 
 
 def read_catalog_file(fname):
+    """
+    Read catalog file
+    :param fname: catalog file name
+    :type fname: str
+    :return: pandas.DataFrame
+    """
     col = ['date', 'evla', 'evlo', 'evdp', 'mag']
     colcata = ['year', 'month', 'day', 'hour', 'minute', 'second']
     data_cata = pd.read_table(fname, header=None, sep='\s+')
diff --git a/seispy/ccp/rf2depth.py b/seispy/ccp/rf2depth.py
deleted file mode 100644
index 0a3d0a08..00000000
--- a/seispy/ccp/rf2depth.py
+++ /dev/null
@@ -1,8 +0,0 @@
-from collections import namedtuple
-
-import numpy as np
-
-
-if __name__ == '__main__':
-    import doctest
-    doctree = doctest.testmod()
\ No newline at end of file
diff --git a/seispy/ccppara.py b/seispy/ccppara.py
index ad1fc927..ca3d8699 100644
--- a/seispy/ccppara.py
+++ b/seispy/ccppara.py
@@ -61,6 +61,12 @@ def shape(self, value):
 
 
 def ccppara(cfg_file):
+    """ Read configure file for CCP stacking
+    :param cfg_file: Path to configure file
+    :type cfg_file: str
+    :return: CCPPara object
+    :rtype: CCPPara
+    """
     cpara = CCPPara()
     cf = configparser.ConfigParser()
     try:
diff --git a/seispy/ccpprofile.py b/seispy/ccpprofile.py
index 012ebb0e..6a66b426 100644
--- a/seispy/ccpprofile.py
+++ b/seispy/ccpprofile.py
@@ -15,7 +15,7 @@
 import sys
 
 
-def prof_range(lat, lon):
+def _prof_range(lat, lon):
     dis = [0]
     for i in range(lat.size-1):
         dis.append(distaz(lat[i], lon[i], lat[i+1], lon[i+1]).degreesToKilometers())
@@ -23,9 +23,19 @@ 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
+    :type val: float
+    :param method: Method for interpolation
+    :type method: str
+    :return: The location of bins (bin_loca), and length between each bin and the start point (profile_range)
+    :rtype: (numpy.array, numpy.array)
+    """
     if not isinstance(stations, Station):
         raise TypeError('Stations should be seispy.rf2depth_makedata.Station')
-    dis_sta = prof_range(stations.stla, stations.stlo)
+    dis_sta = _prof_range(stations.stla, stations.stlo)
     dis_inter = np.append(np.arange(0, dis_sta[-1], val), dis_sta[-1])
     r, theta, phi = geo2sph(np.zeros(stations.stla.size), stations.stla, stations.stlo)
     # t_po = np.arange(stations.stla.size)
@@ -33,11 +43,10 @@ def create_center_bin_profile(stations, val=5, method='linear'):
     theta_i = interp1d(dis_sta, theta, kind=method, bounds_error=False, fill_value='extrapolate')(dis_inter)
     phi_i = interp1d(dis_sta, phi, kind=method, bounds_error=False, fill_value='extrapolate')(dis_inter)
     _, lat, lon = sph2geo(r, theta_i, phi_i)
-    # dis = prof_range(lat, lon)
     return lat, lon, dis_inter
 
 
-def line_proj(lat1, lon1, lat2, lon2):
+def _line_proj(lat1, lon1, lat2, lon2):
     daz = distaz(lat1, lon1, lat2, lon2)
     az1_begin = (daz.baz - 90) % 360
     az2_begin = (daz.baz + 90) % 360
@@ -46,7 +55,7 @@ def line_proj(lat1, lon1, lat2, lon2):
     return az1_begin, az2_begin, az1_end, az2_end, daz
 
 
-def fix_filename(filename, typ='dat'):
+def _fix_filename(filename, typ='dat'):
     dname = dirname(filename)
     if not exists(dname) and dname != '':
         raise FileExistsError('internal error')
@@ -61,7 +70,7 @@ def fix_filename(filename, typ='dat'):
         return filename
 
 
-def bin_shape(cpara):
+def _bin_shape(cpara):
     if cpara.bin_radius is None:
         depmod = DepModel(cpara.stack_range)
         fzone = km2deg(np.sqrt(0.5*cpara.domperiod*depmod.vs*cpara.stack_range))
@@ -98,11 +107,14 @@ def init_profile(lat1, lon1, lat2, lon2, val):
 
 
 class CCPProfile():
-    def __init__(self, cfg_file=None, log=None):
-        if log is None:
-            self.logger = SetupLog()
-        else:
-            self.logger = log
+    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
+        :type log: :class:`seispy.setuplog.SetupLog`
+        """
+        self.logger = log
         if cfg_file is None:
             self.cpara = CCPPara()
         elif isinstance(cfg_file, str):
@@ -149,7 +161,7 @@ def initial_profile(self):
             self.bin_loca = np.vstack((lat, lon)).T
         else:
             self.bin_loca, self.profile_range = init_profile(*self.cpara.line, self.cpara.slide_val)
-        self.fzone = bin_shape(self.cpara)
+        self.fzone = _bin_shape(self.cpara)
         self._get_sta()
         self._select_sta()
     
@@ -199,7 +211,7 @@ def _write_sta(self):
                 self.stalst[idx,0], self.stalst[idx, 1]))
 
     def _proj_sta(self, width):
-        az1_begin, az2_begin, az1_end, az2_end, daz = line_proj(*self.cpara.line)
+        az1_begin, az2_begin, az1_end, az2_end, daz = _line_proj(*self.cpara.line)
         az_sta_begin = distaz(self.stalst[:, 0], self.stalst[:, 1], self.cpara.line[0], self.cpara.line[1]).az
         az_sta_end = distaz(self.stalst[:, 0], self.stalst[:, 1], self.cpara.line[2], self.cpara.line[3]).az
         if 0 <= daz.baz < 90 or 270 <= daz.baz < 360:
@@ -275,7 +287,7 @@ def save_stack_data(self, format='npz'):
         :param format: Format for stacked data
         :type format: str
         """
-        self.cpara.stackfile = fix_filename(self.cpara.stackfile, format)
+        self.cpara.stackfile = _fix_filename(self.cpara.stackfile, format)
         self.logger.CCPlog.info('Saving stacked data to {}'.format(self.cpara.stackfile))
         if not isinstance(self.cpara.stackfile, str):
             self.logger.CCPlog.error('fname should be in \'str\'')
diff --git a/seispy/decon.py b/seispy/decon.py
index 6b6c7df8..22fd5810 100644
--- a/seispy/decon.py
+++ b/seispy/decon.py
@@ -274,6 +274,20 @@ def deconwater(uin, win, dt, tshift=10., wlevel=0.05, f0=2.0, normalize=False, p
 
 
 def deconvolute(uin, win, dt, method='iter', **kwargs):
+    """ Deconvolute receiver function from waveforms.
+    :param uin: R or Q component for the response function
+    :type uin: np.ndarray
+    :param win: Z or L component for the source function
+    :type win: np.ndarray
+    :param dt: sample interval in second
+    :type dt: float
+    :param method: Method for deconvolution, defaults to 'iter'
+    :type method: str, optional
+    :param kwargs: Parameters for deconvolution
+    :type kwargs: dict
+    :return: RF, rms, [iter]
+    :rtype: np.ndarray, float, int
+    """
     if method.lower() == 'iter':
         return deconit(uin, win, dt, **kwargs)
     elif method.lower() == 'water':
diff --git a/seispy/distaz.py b/seispy/distaz.py
index fdf54720..ae0479e2 100644
--- a/seispy/distaz.py
+++ b/seispy/distaz.py
@@ -72,7 +72,7 @@ class distaz:
     and I vaguely remember a perl port. Long live distaz! 
 
     Mijian Xu, Jan 01, 2016
-    Add np.ndarray to availible input varible.
+    Add np.ndarray to available input.
     """
 
     def __init__(self, lat1, lon1, lat2, lon2):
diff --git a/seispy/eq.py b/seispy/eq.py
index 1de7602b..2043b78f 100644
--- a/seispy/eq.py
+++ b/seispy/eq.py
@@ -63,6 +63,9 @@ def __init__(self, pathname, datestr, suffix='SAC'):
         self.it = 0
         self.trigger_shift = 0
         self.inc_correction = 0
+
+    def __str__(self):
+        return('Event data class {0}'.format(self.datestr))
     
     def _check_comp(self):
         if len(self.st) < 3:
@@ -74,18 +77,26 @@ def _check_comp(self):
             pass
 
     def readstream(self):
+        """Read SAC files to stream
+        """
         self.rf = obspy.Stream()
         self.st = obspy.read(self.filestr)
         self._check_comp()
         self.st.sort()
         self.set_comp()
 
-    def cleanstream(self):
+    def _cleanstream(self):
         self.st = None
         self.rf = None
 
     @classmethod
     def from_stream(cls, stream):
+        """Create EQ object from obspy stream
+        :param stream: obspy stream
+        :type stream: obspy.Stream
+        :return: EQ object
+        :rtype: EQ
+        """
         eq = cls('', '')
         eq.st = stream
         eq.datestr = eq.st[0].stats.starttime.strftime('%Y.%j.%H.%M.%S')
@@ -95,6 +106,8 @@ def from_stream(cls, stream):
         return eq
 
     def set_comp(self):
+        """Set component name
+        """
         if self.st.select(channel='*[E2]'):
             self.comp = 'enz'
         elif self.st.select(channel='*R'):
@@ -105,7 +118,7 @@ def set_comp(self):
             raise ValueError('No such component in R, E or Q')
 
     def channel_correct(self, switchEN=False, reverseE=False, reverseN=False):
-        """_summary_
+        """Correct channel name for R, E, N components
 
         :param switchEN: _description_, defaults to False
         :type switchEN: bool, optional
@@ -124,9 +137,6 @@ def channel_correct(self, switchEN=False, reverseE=False, reverseN=False):
             self.st.select(channel='*[E2]')[0].stats.channel = chN
             self.st.select(channel='*[N1]')[0].stats.channel = chE
 
-    def __str__(self):
-        return('Event data class {0}'.format(self.datestr))
-
     def write(self, path, evt_datetime):
         for tr in self.st:
             sac = SACTrace.from_obspy_trace(tr)
@@ -140,14 +150,34 @@ def write(self, path, evt_datetime):
             sac.write(fname)
 
     def detrend(self):
+        """Detrend and demean
+        """
         self.st.detrend(type='linear')
         self.st.detrend(type='constant')
         self.fix_channel_name()
 
     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
+        :type freqmax: float, optional
+        :param order: filter order, defaults to 4
+        :type order: int, optional
+        """
         self.st.filter('bandpass', freqmin=freqmin, freqmax=freqmax, corners=order, zerophase=True)
 
     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
+        :type evdp: float
+        :param dis: epicentral distance
+        :type dis: float
+        :param phase: phase name, defaults to 'P'
+        :type phase: str, optional
+        """
         arrivals = model.get_travel_times(evdp, dis, phase_list=[phase])
         if not arrivals:
             raise ValueError('The phase of {} is not exists'.format(phase))
@@ -161,6 +191,12 @@ def get_arrival(self, model, evdp, dis, phase='P'):
             self.phase = phase
 
     def search_inc(self, bazi):
+        """Search incident angle for S wave
+        :param bazi: back azimuth
+        :type bazi: float
+        :return: incident angle
+        :rtype: float
+        """
         inc_range = np.arange(0.1, 90, 0.1)
         s_range = self.trim(20, 20, isreturn=True)
         power = np.zeros(inc_range.shape[0])
@@ -174,6 +210,18 @@ def search_inc(self, bazi):
         self.inc = real_inc
 
     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
+        :type time_b: int, optional
+        :param time_e: time after P arrival, defaults to 20
+        :type time_e: int, optional
+        :param offset: offset for searching, defaults to 90
+        :type offset: int, optional
+        :return: back azimuth and amplitude
+        :rtype: (float, np.ndarray)
+        """
         p_arr = self.arr_correct(write_to_sac=False)
         this_st = self.st.copy()
         this_st.filter('bandpass', freqmin=0.03, freqmax=0.5)
@@ -194,6 +242,8 @@ def search_baz(self, bazi, time_b=10, time_e=20, offset=90):
         return corr_baz, ampt
 
     def fix_channel_name(self):
+        """Fix channel name for R, E, N components
+        """
         if self.st.select(channel='??1') and self.st.select(channel='??Z') and hasattr(self.st.select(channel='*1')[0].stats.sac, 'cmpaz'):
             if self.st.select(channel='*1')[0].stats.sac.cmpaz == 0:
                 self.st.select(channel='*1')[0].stats.channel = self.st.select(channel='*1')[0].stats.channel[:-1] + 'N'
@@ -210,6 +260,18 @@ def fix_channel_name(self):
             pass
 
     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
+        :type inc: float, optional
+        :param method: method for rotation, defaults to 'NE->RT'
+        :type method: str, optional
+        :param search_inc: whether search incident angle, defaults to False
+        :type search_inc: bool, optional
+        :param baz_shift: shift back azimuth, defaults to 0
+        :type baz_shift: int, optional
+        """
         bazi = np.mod(baz + baz_shift, 360)
         if inc is None:
             if self.phase[-1] == 'S' and search_inc:
@@ -231,6 +293,12 @@ def rotate(self, baz, inc=None, method='NE->RT', search_inc=False, baz_shift=0):
             pass
 
     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
+        :rtype: (float, float, float)
+        """
         st_noise = self.trim(length, 0, isreturn=True)
         st_signal = self.trim(0, length, isreturn=True)
         try:
@@ -248,6 +316,10 @@ def snr(self, length=50):
         return snr_E, snr_N, snr_Z
     
     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
+        """
         if event_time is not None and not isinstance(event_time, obspy.core.utcdatetime.UTCDateTime):
             raise TypeError('Event time should be UTCDateTime type in obspy')
         elif event_time is None:
@@ -274,6 +346,18 @@ def _get_time(self, time_before, time_after):
         return t1, t2
 
     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
+        :type time_after: float
+        :param prepick: whether use prepick, defaults to True
+        :type prepick: bool, optional
+        :param stl: short time length for STA/LTA, defaults to 5
+        :type stl: int, optional
+        :param ltl: long time length for STA/LTA, defaults to 10
+        :type ltl: int, optional
+        """
         t1, t2 = self._get_time(time_before, time_after)
         self.st_pick = self.st.copy().trim(t1, t2)
         if len(self.st_pick) == 0:
@@ -360,6 +444,12 @@ def deconvolute(self, shift, time_after, f0=2.0, method='iter', only_r=False,
                         tr.stats.delta = target_dt
     
     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
+        :type tcomp: bool, optional
+        """
         if self.comp == 'lqt':
             win = self.st.select(channel='*L')[0]
             if tcomp:
@@ -377,6 +467,10 @@ def decon_p(self, tshift, tcomp=False, **kwargs):
         self.rf.append(uout)
 
     def decon_s(self, tshift, **kwargs):
+        """Deconvolution for S wave
+        :param tshift: Time shift before P arrival
+        :type tshift: float
+        """
         if self.comp == 'lqt':
             win = self.st.select(channel='*Q')[0]
             uin = self.st.select(channel='*L')[0]
@@ -389,8 +483,33 @@ def decon_s(self, tshift, **kwargs):
         uout.data = np.flip(uout.data)
         self.rf.append(uout)
 
-    def saverf(self, path, evtstr=None, shift=0, evla=-12345., evlo=-12345., evdp=-12345., mag=-12345.,
+    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
+        :type evtstr: str, optional
+        :param shift: time shift before P arrival, defaults to 0
+        :type shift: int, optional
+        :param evla: event latitude, defaults to -12345.
+        :type evla: float, optional
+        :param evlo: event longitude, defaults to -12345.
+        :type evlo: float, optional
+        :param evdp: event depth, defaults to -12345.
+        :type evdp: float, optional
+        :param mag: event magnitude, defaults to -12345.
+        :type mag: float, optional
+        :param gauss: Gaussian factor, defaults to 0
+        :type gauss: float, optional
+        :param baz: back azimuth, defaults to -12345.
+        :type baz: float, optional
+        :param gcarc: epicentral distance, defaults to -12345.
+        :type gcarc: float, optional
+        :param only_r: whether only save R component, defaults to False
+        :type only_r: bool, optional
+        """
         if self.phase[-1] == 'P':
             if self.comp == 'lqt':
                 svcomp = 'Q'
@@ -431,7 +550,7 @@ def saverf(self, path, evtstr=None, shift=0, evla=-12345., evlo=-12345., evdp=-1
             tr.ka = self.phase
             tr.write(filename + '_{0}_{1}.sac'.format(self.phase, tr.kcmpnm[-1]))
 
-    def s_condition(self, trrf, shift):
+    def _s_condition(self, trrf, shift):
         nt0 = int(np.floor((shift)/trrf.stats.delta))
         nt25 = int(np.floor((shift+25)/trrf.stats.delta))
         if rssq(trrf.data[nt0:nt25]) > rssq(trrf.data[nt25:]):
@@ -440,6 +559,20 @@ def s_condition(self, trrf, shift):
             return False
 
     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
+        :type shift: float
+        :param npts: number of points for RF
+        :type npts: int
+        :param criterion: criterion for judging, defaults to 'crust'
+        :type criterion: str, optional
+        :param rmsgate: RMS gate, defaults to None
+        :type rmsgate: float, optional
+        :return: whether RF is valid
+        :rtype: bool
+        """
         if self.phase[-1] == 'P' and self.comp == 'rtz':
             trrfs = self.rf.select(channel='*R')
         elif self.phase[-1] == 'P' and self.comp == 'lqt':
@@ -501,7 +634,7 @@ def judge_rf(self, gauss, shift, npts, criterion='crust', rmsgate=None):
             else:
                 return False
         elif criterion == "lab":
-            return self.s_condition(trrf, shift) and rengpass
+            return self._s_condition(trrf, shift) and rengpass
         elif criterion is None:
             return rmspass and rengpass
         else:
diff --git a/seispy/geo.py b/seispy/geo.py
index 0a36c16a..f9a8de5a 100644
--- a/seispy/geo.py
+++ b/seispy/geo.py
@@ -4,41 +4,89 @@
 from seispy.utils import scalar_instance, array_instance
 
 def sind(deg):
+    """ Sine function with degree as input
+    :param deg: Degree
+    :type deg: float
+    :return: Sine value
+    :rtype: float
+    """
     rad = np.radians(deg)
     return np.sin(rad)
 
 
 def cosd(deg):
+    """ Cosine function with degree as input
+    :param deg: Degree
+    :type deg: float
+    :return: Cosine value
+    :rtype: float
+    """
     rad = np.radians(deg)
     return np.cos(rad)
 
 
 def tand(deg):
+    """ Tangent function with degree as input
+    :param deg: Degree
+    :type deg: float
+    :return: Tangent value
+    :rtype: float
+    """
     rad = np.radians(deg)
     return np.tan(rad)
 
 
 def cotd(deg):
+    """ Cotangent function with degree as input
+    :param deg: Degree
+    :type deg: float
+    :return: Cotangent value
+    :rtype: float
+    """
     rad = np.radians(deg)
     return np.cos(rad) / np.sin(rad)
 
 
 def asind(x):
+    """ Inverse sine function with degree as output
+    :param x: Sine value
+    :type x: float
+    :return: Degree
+    :rtype: float
+    """
     rad = np.arcsin(x)
     return np.degrees(rad)
 
 
 def acosd(x):
+    """ Inverse cosine function with degree as output
+    :param x: Cosine value
+    :type x: float
+    :return: Degree 
+    :rtype: float
+    """
     rad = np.arccos(x)
     return np.degrees(rad)
 
 
 def atand(x):
+    """ Inverse tangent function with degree as output
+    :param x: Tangent value
+    :type x: float
+    :return: Degree
+    :rtype: float
+    """
     rad = np.arctan(x)
     return np.degrees(rad)
 
 
 def km2deg(km):
+    """ Convert km to degree
+    :param km: Distance in km
+    :type km: float
+    :return: Distance in degree
+    :rtype: float
+    """
     radius = 6371
     circum = 2*np.pi*radius
     conv = circum / 360
@@ -47,6 +95,12 @@ def km2deg(km):
 
 
 def deg2km(deg):
+    """ Convert degree to km
+    :param deg: Distance in degree
+    :type deg: float
+    :return: Distance in km
+    :rtype: float
+    """
     radius = 6371
     circum = 2*np.pi*radius
     conv = circum / 360
@@ -55,26 +109,56 @@ def deg2km(deg):
 
 
 def rad2deg(rad):
+    """ Convert radian to degree
+    :param rad: Radian
+    :type rad: float
+    :return: Degree
+    :rtype: float
+    """
     deg = rad*(360/(2*np.pi))
     return deg
 
 
 def skm2sdeg(skm):
+    """ Convert s/km to s/degree
+    :param skm: s/km
+    :type skm: float
+    :return: s/degree
+    :rtype: float
+    """
     sdeg = skm * deg2km(1)
     return sdeg
 
 
 def sdeg2skm(sdeg):
+    """ Convert s/degree to s/km
+    :param sdeg: s/degree
+    :type sdeg: float
+    :return: s/km
+    :rtype: float
+    """
     skm = sdeg / deg2km(1)
     return skm
 
 
 def srad2skm(srad):
+    """ Convert s/rad to s/km
+    :param srad: s/rad
+    :type srad: float
+    :return: s/km
+    :rtype: float
+    """
     sdeg = srad * ((2*np.pi)/360)
     return sdeg / deg2km(1)
 
 
 def skm2srad(skm):
+    """ Convert s/km to s/rad
+    :param skm: s/km
+    :type skm: float
+    :return: s/rad
+    :rtype: float
+    """
     sdeg = skm * deg2km(1)
     return rad2deg(sdeg)
 
@@ -88,8 +172,8 @@ def rot3D(bazi, inc):
     
     :param bazi: back-azimuth of station-event pair
     :param inc: Incidence angle
-    :return: Coefficient matrix M
-    
+    :return: Coefficient matrix m
+    :rtype: np.ndarray
     """
 
     if scalar_instance(inc):
@@ -102,20 +186,44 @@ def rot3D(bazi, inc):
     inc = inc / 180 * np.pi
     bazi = bazi / 180 * np.pi
 
-    M = np.array([[np.cos(inc), -np.sin(inc)*np.sin(bazi), -np.sin(inc)*np.cos(bazi)],
+    m = np.array([[np.cos(inc), -np.sin(inc)*np.sin(bazi), -np.sin(inc)*np.cos(bazi)],
                   [np.sin(inc), np.cos(inc)*np.sin(bazi), np.cos(inc)*np.cos(bazi)],
                   [value31, -np.cos(bazi), np.sin(bazi)]])
-    return M
-
-
-def rotateSeisZENtoLQT(Z, E, N, bazi, inc):
-    M = rot3D(bazi, inc)
-    ZEN = np.array([Z, E, N])
-    LQT = np.dot(M, ZEN)
-    return LQT[0, :], LQT[1, :], LQT[2, :]
+    return m
+
+
+def rotateSeisZENtoLQT(z, e, n, bazi, inc):
+    """ Rotate ZEN to LQT
+    :param z: Vertical component
+    :type z: np.ndarray
+    :param e: East component
+    :type e: np.ndarray
+    :param n: North component
+    :type n: np.ndarray
+    :param bazi: Back-azimuth
+    :type bazi: float
+    :param inc: Incidence angle
+    :type inc: float
+    :return: L, Q and T components
+    :rtype: tuple
+    """
+    m = rot3D(bazi, inc)
+    zen = np.array([z, e, n])
+    lqt = np.dot(m, zen)
+    return lqt[0, :], lqt[1, :], lqt[2, :]
 
 
 def spherical2cartesian(lon, lat, dep):
+    """ Convert spherical coordinates to cartesian coordinates
+    :param lon: Longitude
+    :type lon: float
+    :param lat: Latitude
+    :type lat: float
+    :param dep: Depth
+    :type dep: float
+    :return: Cartesian coordinates
+    :rtype: tuple
+    """
     cola = 90. - lat
     r = 6371 - dep
     x = r * sind(cola) * cosd(lon)
@@ -124,18 +232,42 @@ def spherical2cartesian(lon, lat, dep):
     return x, y, z
 
 
-def rotateSeisENtoTR(E, N, BAZ):
-    angle = np.mod(BAZ+180, 360)
-    R = N*cosd(angle) + E*sind(angle)
-    T = E*cosd(angle) - N*sind(angle)
-    return T, R
+def rotateSeisENtoTR(e, n, baz):
+    """ Rotate EN to TR
+    :param e: East component
+    :type e: np.ndarray
+    :param n: North component
+    :type n: np.ndarray
+    :param baz: Back-azimuth
+    :type baz: float
+    :return: T and R components
+    :rtype: tuple
+    """
+    angle = np.mod(baz+180, 360)
+    r = n*cosd(angle) + e*sind(angle)
+    t = e*cosd(angle) - n*sind(angle)
+    return t, r
 
 
 def rssq(x):
+    """ Root sum square
+    :param x: Input array
+    :type x: np.ndarray
+    :return: Root sum square
+    :rtype: float
+    """
     return np.sqrt(np.sum(x**2)/len(x))
 
 
 def snr(x, y):
+    """ Signal-to-noise ratio
+    :param x: Signal
+    :type x: np.ndarray
+    :param y: Noise
+    :type y: np.ndarray
+    :return: SNR
+    :rtype: float
+    """
     spow = rssq(x)**2
     npow = rssq(y)**2
     if npow == 0:
@@ -225,6 +357,22 @@ def init_lalo(lat0, lon0, npts):
 
 
 def geoproject(lat_p, lon_p, lat1, lon1, lat2, lon2):
+    """ Project a point to a line
+    :param lat_p: Latitude of the point
+    :type lat_p: float
+    :param lon_p: Longitude of the point
+    :type lon_p: float
+    :param lat1: Latitude of the first point of the line
+    :type lat1: float
+    :param lon1: Longitude of the first point of the line
+    :type lon1: float
+    :param lat2: Latitude of the second point of the line
+    :type lat2: float
+    :param lon2: Longitude of the second point of the line
+    :type lon2: float
+    :return: Latitude and longitude of the projected point
+    :rtype: tuple
+    """
     azi = distaz(lat1, lon1, lat2, lon2).baz
     dis_center = distaz(lat1, lon1, lat_p, lon_p).delta
     azi_center = distaz(lat1, lon1, lat_p, lon_p).baz
diff --git a/seispy/harmonics.py b/seispy/harmonics.py
index f3c8a901..8ce7a9ba 100644
--- a/seispy/harmonics.py
+++ b/seispy/harmonics.py
@@ -44,6 +44,8 @@ def cut_trace(self):
         self.time_axis = np.linspace(self.tmin, self.tmax, self.nsamp)
 
     def harmo_trans(self):
+        """Harmonic decomposition for extracting anisotropic and isotropic features from the radial and transverse RFs
+        """
         self.traces = np.vstack((self.datar, self.datat))
         harmonic = np.zeros((self.trace_num*2, 5))
         unmodel = np.zeros((self.trace_num*2, 5))
diff --git a/seispy/para.py b/seispy/para.py
index a245a690..0a69ec00 100644
--- a/seispy/para.py
+++ b/seispy/para.py
@@ -75,6 +75,8 @@ def get_station_from_ws(self, **kwargs):
 
 class RFPara(object):
     def __init__(self):
+        """Parameters for receiver function calculation
+        """
         self.datapath = expanduser('~')
         self.rfpath = expanduser('~')
         self.catalogpath = join(dirname(__file__), 'data', 'EventCMT.dat')
@@ -221,6 +223,10 @@ def decon_method(self, value):
 
     @classmethod
     def read_para(cls, cfg_file):
+        """Read parameters from configure file
+        :param cfg_file: Path to configure file
+        :type cfg_file: str
+        """
         cf = configparser.RawConfigParser(allow_no_value=True)
         pa = cls()
         try:
diff --git a/seispy/plotR.py b/seispy/plotR.py
index 3824d933..e1c7c76e 100644
--- a/seispy/plotR.py
+++ b/seispy/plotR.py
@@ -26,6 +26,8 @@ def __init__(self, stadata, enf=12, xlim=[2, 80], key='bazi'):
         self.init_figure()
 
     def init_figure(self):
+        """Initialize figure
+        """
         self.fig = plt.figure(figsize=(8, 10))
         self.fig.tight_layout()
         gs = GridSpec(15, 3)
@@ -38,6 +40,8 @@ def init_figure(self):
         self.axs.grid(color='gray', linestyle='--', linewidth=0.4, axis='x')
 
     def plot_waves(self):
+        """Plot PRFs with R component
+        """
         bound = np.zeros(self.stadata.rflength)
         for i in range(self.stadata.ev_num):
             datar = self.stadata.data_prime[i] * self.enf + (i + 1)
@@ -57,6 +61,8 @@ def plot_waves(self):
         # return axp
 
     def set_fig(self):
+        """Set figure
+        """
         y_range = np.arange(self.stadata.ev_num) + 1
         x_range = np.arange(0, self.xlim[1]+2, 5)
         space = 2
diff --git a/seispy/rf2depth_makedata.py b/seispy/rf2depth_makedata.py
index 87d7f92f..5e31fb6b 100644
--- a/seispy/rf2depth_makedata.py
+++ b/seispy/rf2depth_makedata.py
@@ -24,27 +24,6 @@ class sta_part, sta_full, _RFInd
 # sta_part is not used
 sta_part = namedtuple('sta_part', ['station', 'stla', 'stlo'])
 sta_full = namedtuple('sta_full', ['station', 'stla', 'stlo', 'stel'])
-_RfInDepth = namedtuple('RfInDepth',
-                       ['station', 'stalat', 'stalon',
-                        'depthrange', 'bazi', 'rayp',
-                        'moveout_correct', 'piercelat', 'piercelon', 'stopindex'])
-class RFInDepth(_RfInDepth):
-    """
-    this class is inheried from a namedtuple, and can be called as a dict
-    mainly to restrict behaviour while output  and
-    provide a better speed while
-
-    warnings: this class is still in develpment
-    """
-    def __getitem__(self, item):
-        """
-        make RFInDepth can be called like dict
-        >>> ll = RFInDepth(station='sta', stalat='stla', stalon='stlo', depthrange='dep', bazi='baz',rayp='rayp', moveout_correct=0, piercelat=0, piercelon=0 )
-        >>> ll['stalat'], ll.stalat
-        ('stla', 'stla')
-        """
-        return getattr(self, item)
-
 
 class Station(object):
     def __init__(self, sta_lst:str):
@@ -62,9 +41,10 @@ def __init__(self, sta_lst:str):
             self.station, self.stla, self.stlo = np.loadtxt(sta_lst, dtype=dtype, unpack=True, ndmin=1)
             self.stel = np.zeros(self.stla.size)
         self.sta_num = self.stla.shape[0]
+
     def __getitem__(self,index):
         """
-        allows  [] do indexing
+        allow for sta[index]
         """
         if hasattr(self,'stel'):
             return sta_full(self.station[index],
@@ -72,9 +52,10 @@ def __getitem__(self,index):
         else:
             return sta_full(self.station[index],
                             self.stla[index], self.stlo[index], 0)
+
     def __iter__(self):
         """
-        allow for...in to iter
+        allow for sta in stalist
         """
         for index in range(self.stla.shape[0]):
             if hasattr(self, 'stel'):
@@ -91,14 +72,14 @@ class RFDepth():
     """Convert receiver function to depth axis
     """
     def __init__(self, cpara:CCPPara, log:Logger=SetupLog().RF2depthlog,
-                 RAYTRACING3D=False, velmod3d=None, modfolder1d=None) -> None:
+                 raytracing3d=False, velmod3d=None, modfolder1d=None) -> None:
         """
         :param cpara: CCPPara object
         :type cpara: CCPPara
         :param log: Log object
         :type log: logging.Logger
-        :param RAYTRACING3D: If True, use 3D ray tracing to calculate the travel time
-        :type RAYTRACING3D: bool
+        :param raytracing3d: If True, use 3D ray tracing to calculate the travel time
+        :type raytracing3d: bool
         :param velmod3d: Path to 3D velocity model in npz file
         :type velmod3d: str
         :param modfolder1d: Folder path to 1D velocity model files with staname.vel as the file name
@@ -108,7 +89,7 @@ def __init__(self, cpara:CCPPara, log:Logger=SetupLog().RF2depthlog,
         self.cpara = cpara
         self.modfolder1d = modfolder1d
         self.log = log
-        self.RAYTRACING3D = RAYTRACING3D
+        self.raytracing3d = raytracing3d
         if velmod3d is not None:
             if isinstance(velmod3d, str):
                 self.mod3d = Mod3DPerturbation(velmod3d, cpara.depth_axis, velmod=cpara.velmod)
@@ -166,7 +147,6 @@ def makedata(self, psphase=1):
                 sphere = False
             self.log.info('the {}th/{} station with {} events'.format(_i + 1, len(self.sta_info), stadatar.ev_num))
 
-
             #### 1d model for each station
             if self.ismod1d:
                 if self.modfolder1d is not None:
@@ -200,15 +180,15 @@ def makedata(self, psphase=1):
 
                 ps_rfdepth, end_index = time2depth(stadatar, self.cpara.depth_axis, newtpds)
 
-            rfdep = self._write_rfdep(self.cpara,stadatar, ps_rfdepth, piercelat, piercelon, end_index)
+            self._write_rfdep(stadatar, ps_rfdepth, piercelat, piercelon, end_index)
         np.save(self.cpara.depthdat, self.rfdepth)
 
-    def _write_rfdep(self,cpara, stadata, amp, pplat, pplon, end_index):
+    def _write_rfdep(self, stadata, amp, pplat, pplon, end_index):
         rfdep = {}
         rfdep['station'] = stadata.staname
         rfdep['stalat'] = stadata.stla
         rfdep['stalon'] = stadata.stlo
-        rfdep['depthrange'] = cpara.depth_axis
+        rfdep['depthrange'] = self.cpara.depth_axis
         rfdep['bazi'] = stadata.bazi
         rfdep['rayp'] = stadata.rayp
         rfdep['moveout_correct'] = amp
@@ -218,7 +198,6 @@ def _write_rfdep(self,cpara, stadata, amp, pplat, pplon, end_index):
         self.rfdepth.append(rfdep)
 
 
-
 def rf2depth():
     """
     CLI for Convert receiver function to depth axis
@@ -229,7 +208,6 @@ def rf2depth():
     3. -d and -r : do moveout correction and raytracing
     4. -m : use {staname}.vel file for RF2depth conversion
 
-
     """
     parser = argparse.ArgumentParser(description="Convert Ps RF to depth axis")
     parser.add_argument('-d',
@@ -271,7 +249,7 @@ def rf2depth():
         modfolder1d = None
     rfd = RFDepth(
         cpara,
-        RAYTRACING3D=raytracing3d,
+        raytracing3d=raytracing3d,
         velmod3d=velmod3d,
         modfolder1d=modfolder1d,
     )