debug latlon_from

xumi1993 · Dec 22, 2023 · 3065e16 · 3065e16
1 parent 919f428
commit 3065e16
Show file tree

Hide file tree

Showing 2 changed files with 6 additions and 5 deletions.
diff --git a/seispy/geo.py b/seispy/geo.py
@@ -205,8 +205,8 @@ def init_lalo(lat0, lon0, npts):
             if len(lat0) != len(lon0):
                 raise ValueError('lat0 and lon0 must be in the same length')
             else:
-                azimuth = np.ones(lat0)*azimuth
-                gcarc_dist = np.ones(lat0)*gcarc_dist
+                azimuth = np.ones_like(lat0)*azimuth
+                gcarc_dist = np.ones_like(lat0)*gcarc_dist
         elif scalar_instance(lat0) and scalar_instance(lon0):
             pass
         else:
@@ -219,7 +219,6 @@ def init_lalo(lat0, lon0, npts):
         npts = len(azimuth)
         gcarc_dist = np.ones(lat0, lon0, npts)*gcarc_dist
         lat0, lon0 = init_lalo(lat0, lon0, npts)
-
     g = Geod(ellps=ellps)
     lon, lat, _ = g.fwd(lon0, lat0, azimuth, deg2km(gcarc_dist)*1000)
     return lat, lon

diff --git a/seispy/rf2depth_makedata.py b/seispy/rf2depth_makedata.py
@@ -132,8 +132,10 @@ def makedata(self, psphase=1):
                     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)
-                piercelat, piercelon = latlon_from(self.sta_info.stla[i], self.sta_info.stlo[i],
-                                                            stadatar.bazi, rad2deg(x_s))
+                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:
                 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)