fixed hough transform for gaussian grids

martalmeida · Nov 13, 2020 · 7c00ef3 · 7c00ef3
1 parent 45eec3e
commit 7c00ef3
Show file tree

Hide file tree

Showing 3 changed files with 18 additions and 8 deletions.
diff --git a/nmf3d/transforms.py b/nmf3d/transforms.py
@@ -241,14 +241,15 @@ def hough(HOUGHs_UVZ,W_nk,Lat,ws0):
     for n in range(nM): # wavenumber index
       for l in range(nL): # meridional index
         for t in range(nTimes): # time index
-          Aux = W_nk[:,k,:,n,t]*np.conjugate(HOUGHs_UVZ[:,n,l,k,:])*cosTheta   # Aux(3,Lat)
-          y1=Aux.sum(0) # Integrand -> y1(Lat)
-
           if latType=='linear':
+            Aux = W_nk[:,k,:,n,t]*np.conjugate(HOUGHs_UVZ[:,n,l,k,:])*cosTheta # Aux(3,Lat)
+            y1=Aux.sum(0) # Integrand -> y1(Lat)
             # Computes latitude integral of the Integrand using trapezoidal method
             aux1=(y1[:-1]+y1[1:])*Dl/2. # len Lat.size-1; complex
-
           elif latType=='gaussian':
+            Aux = W_nk[:,k,:,n,t]*np.conjugate(HOUGHs_UVZ[:,n,l,k,:]) # Aux(3,Lat)
+            y1=Aux.sum(0) # Integrand -> y1(Lat)
+            # Computes latitude integral of the Integrand using gaussian quadrature
             aux1=gw*y1
 
           w_nlk[k,n,l,t] = aux1.sum()

diff --git a/nmf3d_mat/hough_transform.m b/nmf3d_mat/hough_transform.m
@@ -33,14 +33,16 @@
   for n=1:nM % wavenumber index
     for l=1:nL % meridional index
       for t=1:nTimes % time index
-        Aux=squeeze(var_nk(:,k,:,n,t)).*squeeze(conj(HOUGHs_UVZ(:,n,l,k,:))).*cosTheta; % Aux(3,Lat)
-        y1=sum(Aux); % Integrand -> y1(Lat)
-
         if isequal(latType,'linear')
+          Aux=squeeze(var_nk(:,k,:,n,t)).*squeeze(conj(HOUGHs_UVZ(:,n,l,k,:))).*cosTheta; % Aux(3,Lat)
+          y1=sum(Aux); % Integrand -> y1(Lat)
           % Computes latitude integral of the Integrand using trapezoidal method
           aux1=(y1(1:end-1)+y1(2:end))*Dl/2.; % len Lat.size-1; complex
         elseif isequal(latType,'gaussian')
-          aux1=gw*y1;
+          Aux=squeeze(var_nk(:,k,:,n,t)).*squeeze(conj(HOUGHs_UVZ(:,n,l,k,:))); % Aux(3,Lat)
+          y1=sum(Aux); % Integrand -> y1(Lat)
+          % Computes latitude integral of the Integrand using gaussian quadrature
+          aux1=gw.*y1.';
         end
 
         var_nlk(k,n,l,t) = sum(aux1);

diff --git a/nmf3d_mat/lgwt.m b/nmf3d_mat/lgwt.m
@@ -11,6 +11,9 @@
 % the definite integral using sum(f.*w);
 %
 % Written by Greg von Winckel - 02/25/2004
+%
+% changed, nov-2020, to return increasing x
+
 N=N-1;
 N1=N+1; N2=N+2;
 
@@ -54,5 +57,9 @@
 % Linear map from[-1,1] to [a,b]
 x=(a*(1-y)+b*(1+y))/2;      
 
+% changed nov 2020:
+% note that w is always symmetric, no need to change
+x=x(end:-1:1);
+
 % Compute the weights
 w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;