From ec1a456f7ba3d4d455e2e876fa246b54780237b5 Mon Sep 17 00:00:00 2001 From: Luca De Siena Date: Thu, 19 Sep 2024 10:50:42 +0200 Subject: [PATCH] Interpolation on segments Changes to MuRAT_segments to + remove detection trade offs + correct length computation for segments --- MuRAT3.mlx | Bin 4972 -> 4968 bytes bin/Murat_segments.m | 93 +++++++++++++++++++++++-------------------- 2 files changed, 50 insertions(+), 43 deletions(-) diff --git a/MuRAT3.mlx b/MuRAT3.mlx index 678843b0d606ed2663eacbbb56d58ca7e5e4ac27..b61c01452d4d72054e6bb61c48ead8e23a34791e 100644 GIT binary patch delta 363 zcmV-x0hIpiCg>)xhYST$tqe;VlZ*^Lf5RXYhVT0ng!>lgtWMH)7n)sZVq%TPJ3`Mk zVxa`kZf`HETXnj`OY(o`%i;XxVyip!0nV%5RG4KHBQUaWm04H#sad9TjC>HL6}TT2$csrmkvaL7DWgD+U=M;&u(s;C9ZRt z)0E{Yy>3{+3&vUTT9RM$Cr%wJ7+R_Sq!-bcva5#8Iisi2qaapE{&ZOz%Rud2g J8W8{h005Tyotppv delta 392 zcmaE%_C{?(I}huJJ%+w7swFTJ}Bm_CXGl+dD24EbZ(xmU=(y9Yc5a z^(|Kq-i}uD2@%qt{ogv?@BIAr-*&mmKd8Aplj(CyB8yvTY4SGKx>eUw8Yl8!y0Yml zkH`-Oxy7NDC#Q<8D*o`mid)a*L7<7G{uJF0dMw%+U$fYweyIqlwaq@4Df#Gf?#ik4 zuYB_wm-U=o*`|;!sbb63Dck#F$(pMxw!B{_%d&ORg(fy(o2Q)1@2|-_bnwvC&f*G3 z_JUjSt?el)-}gzb-1mNRomD|$y}U2bdp1}mMjoBkbpA+5p`tB! z($p{h)1N*O-u!^qg#{W`nnH2RPk0$7|K*jN94jO+nTwx~p8*UQArwDIPHOUdA#-LX b0SMbvfL|=Yo0ScufCC6Em>3x31VKCiV}7Nq diff --git a/bin/Murat_segments.m b/bin/Murat_segments.m index 4c221b6..09d31f1 100644 --- a/bin/Murat_segments.m +++ b/bin/Murat_segments.m @@ -6,8 +6,8 @@ % CREATES variables to compute ray-dependent inversion matrices and more % % Input parameters: -% modv: velocity model -% rma: ray in MuRAt format +% modv: velocity model +% rma: ray in MuRAT format % % Output parameters: % lunpar: the lenght of the segments in each block @@ -25,7 +25,8 @@ % The code takes two points of the ray and checks if a % boundary is crossed. If not, the length of the segment is added to % the previous part. It the bounday is crossed, the segment is broken, -% and a new segment starts. +% and a new segment starts. The ray gets interpolated into smaller +% segments for hit check. lunpar = zeros(100,1); blocch = zeros(100,1); @@ -63,13 +64,18 @@ % Check that all x-y elements are more than zero fallr = find(rma(:,2)>0 & rma(:,3)>0); -% Define the vectors -xpolo = rma(fallr,2); -ypolo = rma(fallr,3); +% Define the interpolated vectors for hit count +xpolo_noint = rma(fallr,2); +ypolo_noint = rma(fallr,3); +zpolo_noint = rma(fallr,4); -% Need to change as rma is with positive depths -zpolo = rma(fallr,4); -tot = length(xpolo); +tot = 1000; +xpolo =... + linspace(xpolo_noint(1),xpolo_noint(end),tot); +ypolo =... + linspace(ypolo_noint(1),ypolo_noint(end),tot); +zpolo =... + linspace(zpolo_noint(1),zpolo_noint(end),tot); % Read hypocenter coordinates sorg = [xpolo(1) ypolo(1) zpolo(1)]; @@ -117,16 +123,18 @@ cx = (ewS <= BLx & ewR > BLx) |... (ewS >=BLx & ewR < BLx); + [a_x, b_x] = find(cx); cy = (nsS <= BLy & nsR > BLy) |... (nsS >=BLy & nsR < BLy); + [a_y, b_y] = find(cy); cv = (vS <= BLv & vR > BLv) |... (vS >= BLv & vR < BLv); + [a_z, b_z] = find(cv); % x crossing - if find(cx) > 0 - n = find(cx>0, 1, 'first'); + if a_x > 0 index = index + 1; - inte(index,1:5) = [BLx(n),nsR,vR,index-1,0]; + inte(index,1:5) = [BLx(b_x),nsR,vR,index-1,0]; lung = sum(diffLoc(k1:k)); k1 = k; inte(index,5) = lung; @@ -143,13 +151,12 @@ break end end - end - - % y crossing - if find(cy) > 0 - n = find(cy>0, 1, 'first'); + end + + % y crossing + if a_y > 0 index = index + 1; - inte(index,1:5) = [ewR,BLy(n),vR,index-1,0]; + inte(index,1:5) = [ewR,BLy(b_y),vR,index-1,0]; lung = sum(diffLoc(k1:k)); k1 = k; inte(index,5) = lung; @@ -166,14 +173,14 @@ break end end - end - - % z crossing - if find(cv) > 0 - n = find(cv>0, 1, 'first'); + end + + % z crossing + if a_z > 0 index = index + 1; - inte(index,1:5) = [ewR,nsR,BLv(n),index-1,0]; + inte(index,1:5) = [ewR,nsR,BLv(b_z),index-1,0]; lung = sum(diffLoc(k1:k)); + k1 = k; inte(index,5) = lung; if vS <= vR if inte(index,6)+passoz < max(bv) @@ -188,43 +195,43 @@ break end end - + end continue end - + % In case the receiver is in the grid. if k == tot - + % receiver location - bR = xv < ricev(1) & ricev(1)0); + bRR = bv(bR>0); if isempty(bRR) == 0 - index = index + 1; - lung = sum(diffLoc(k1:end)); - inte(index,1:5) = [ricev(1),ricev(2),ricev(3),index-1,lung]; + index = index + 1; + lung = sum(diffLoc(k1:end)); + inte(index,1:5) = [ricev(1),ricev(2),ricev(3),index-1,lung]; end end % Variable inte monitors the points where the rays cross the grid -linte = length(inte(:,1)); +linte = length(inte(:,1)); -lunpar(1:linte,1) = inte(:,5); -blocch(1:linte,1) = inte(:,6); -lunto(1,1) = sum(lunpar(:,1)); +lunpar(1:linte,1) = inte(:,5); +blocch(1:linte,1) = inte(:,6); +lunto(1,1) = sum(lunpar(:,1)); -s = zeros(100,1); -no = blocch>length(modv(:,1)); -blocch(no) = 0; -lunpar(no) = 0; -bff = find(blocch); +s = zeros(100,1); +no = blocch>length(modv(:,1)); +blocch(no) = 0; +lunpar(no) = 0; +bff = find(blocch); if bff > 0 - s(bff) = 1./modv(blocch(bff),4); + s(bff) = 1./modv(blocch(bff),4); rayCrossing(blocch(bff))= rayCrossing(blocch(bff))+1; end - + end