fix toeplitz issues by just using existing code and pasting it into the full array

Author: xzackli

pspy/mcm_fortran/mcm_compute.mod

@@ -414,9 +414,9 @@ subroutine calc_coupling_block_spin0(wcl, l_exact, l_band, l_toeplitz, coupling)
     nlmax = size(coupling,1) - 1
 
     !$omp parallel do private(l2, l1) schedule(dynamic)
-    do l1 = 0, min(nlmax, l_exact)
+    do l1 = 2, min(nlmax, l_exact)
         do l2 = l1, nlmax
-            call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1+1, l2+1))
+            call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1-1, l2-1))
         end do
     end do
 
@@ -430,14 +430,14 @@ subroutine calc_coupling_block_spin0(wcl, l_exact, l_band, l_toeplitz, coupling)
             end if
 
             do l2 = l1, lmax_band
-                call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1+1, l2+1))
+                call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1-1, l2-1))
             end do
         end do
 
         if (l_toeplitz .lt. nlmax) then
             !$omp parallel do
             do l1 = l_toeplitz + 1, nlmax
-                call calc_coupling_elem_spin0(wcl, l1, l1, coupling(l1+1, l1+1))
+                call calc_coupling_elem_spin0(wcl, l1, l1, coupling(l1-1, l1-1))
             end do
         end if
     end if
@@ -457,11 +457,9 @@ subroutine calc_coupling_block_spin02(wcl, l_exact, l_band, l_toeplitz, coupling
     !$omp parallel do private(l2, l1) schedule(dynamic)
     do l1 = 2, min(nlmax, l_exact)
         do l2 = l1, nlmax
-            call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1+1, l2+1))
+            call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1-1, l2-1))
         end do
     end do
-    coupling(1,1) = 1
-    coupling(2,2) = 1
 
     if (l_exact .lt. nlmax) then
         !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic)
@@ -473,14 +471,14 @@ subroutine calc_coupling_block_spin02(wcl, l_exact, l_band, l_toeplitz, coupling
             end if
 
             do l2 = l1, lmax_band
-                call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1+1, l2+1))
+                call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1-1, l2-1))
             end do
         end do
 
         if (l_toeplitz .lt. nlmax) then
             !$omp parallel do
             do l1 = l_toeplitz + 1, nlmax
-                call calc_coupling_elem_spin02(wcl, l1, l1, coupling(l1+1, l1+1))
+                call calc_coupling_elem_spin02(wcl, l1, l1, coupling(l1-1, l1-1))
             end do
         end if
     end if
@@ -500,11 +498,9 @@ subroutine calc_coupling_block_spin2_pp(wcl, l_exact, l_band, l_toeplitz, coupli
     !$omp parallel do private(l2, l1) schedule(dynamic)
     do l1 = 2, min(nlmax, l_exact)
         do l2 = l1, nlmax
-            call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1+1, l2+1))
+            call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1-1, l2-1))
         end do
     end do
-    coupling(1,1) = 1
-    coupling(2,2) = 1
 
     if (l_exact .lt. nlmax) then
         !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic)
@@ -516,14 +512,14 @@ subroutine calc_coupling_block_spin2_pp(wcl, l_exact, l_band, l_toeplitz, coupli
             end if
 
             do l2 = l1, lmax_band
-                call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1+1, l2+1))
+                call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1-1, l2-1))
             end do
         end do
 
         if (l_toeplitz .lt. nlmax) then
             !$omp parallel do
             do l1 = l_toeplitz + 1, nlmax
-                call calc_coupling_elem_spin2_pp(wcl, l1, l1, coupling(l1+1, l1+1))
+                call calc_coupling_elem_spin2_pp(wcl, l1, l1, coupling(l1-1, l1-1))
             end do
         end if
     end if
@@ -543,11 +539,9 @@ subroutine calc_coupling_block_spin2_mm(wcl, l_exact, l_band, l_toeplitz, coupli
     !$omp parallel do private(l2, l1) schedule(dynamic)
     do l1 = 2, min(nlmax, l_exact)
         do l2 = l1, nlmax
-            call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1+1, l2+1))
+            call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1-1, l2-1))
         end do
     end do
-    coupling(1,1) = 1
-    coupling(2,2) = 1
 
     if (l_exact .lt. nlmax) then
         !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic)
@@ -559,14 +553,14 @@ subroutine calc_coupling_block_spin2_mm(wcl, l_exact, l_band, l_toeplitz, coupli
             end if
 
             do l2 = l1, lmax_band
-                call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1+1, l2+1))
+                call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1-1, l2-1))
             end do
         end do
 
         if (l_toeplitz .lt. nlmax) then
             !$omp parallel do
             do l1 = l_toeplitz + 1, nlmax
-                call calc_coupling_elem_spin2_mm(wcl, l1, l1, coupling(l1+1, l1+1))
+                call calc_coupling_elem_spin2_mm(wcl, l1, l1, coupling(l1-1, l1-1))
             end do
         end if
     end if

pspy/mcm_fortran/mcm_fortran.f90

@@ -522,8 +522,7 @@ def coupling_block(
     l = np.arange(len(wcl))
     wcl *= (2 * l + 1)
 
-
-    mcm = np.zeros((lmax+1, lmax+1))
+    mcm = np.zeros((lmax, lmax))
 
     if l_toep is None: l_toep = lmax
     if l_band is None: l_band = lmax
@@ -542,8 +541,13 @@ def coupling_block(
         mcm = format_toepliz_fortran2(mcm, l_toep, l_exact, lmax)
 
     mcm_fortran.fill_upper(mcm.T)
+    
+    mcm_full = np.zeros((lmax+1, lmax+1))
+    mcm_full[2:, 2:] = mcm[:-1,:-1]
+    mcm_full[0,0] = 1.0
+    mcm_full[1,1] = 1.0
 
     if legacy_ell_range:
-        return mcm[2:-1, 2:-1]
-    else:
         return mcm
+    else:
+        return mcm_full