include max_mode and fixed error due to NaN of model system

JuliaWaveScattering · Aug 21, 2024 · 4042798 · 4042798
2 parents 0f81ef3 + 05aa027
commit 4042798
Show file tree

Hide file tree

Showing 11 changed files with 804 additions and 184 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "EffectiveWaves"
 uuid = "37e8709b-1ed2-53db-b26a-3571262b3cb4"
 authors = ["Artur L. Gower <[email protected]>"]
-version = "0.3.5"
+version = "0.3.6"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
@@ -30,7 +30,7 @@ HCubature = "1"
 Interpolations = "0.12, ^0.13"
 IterTools = "1"
 MDBM = "0.1"
-MultipleScattering = "^0.1.17"
+MultipleScattering = "^0.1.21"
 OffsetArrays = "^0.23, 1"
 Optim = "0.22, 1"
 ParticleCorrelations = "^0.0.1"

diff --git a/docs/src/theory/reflection_transmission_low_freq.nb b/docs/src/theory/reflection_transmission_low_freq.nb
@@ -10,10 +10,17 @@
 NotebookFileLineBreakTest
 NotebookFileLineBreakTest
 NotebookDataPosition[       158,          7]
+<<<<<<< HEAD
 NotebookDataLength[     46591,       1235]
 NotebookOptionsPosition[     43605,       1177]
 NotebookOutlinePosition[     43993,       1193]
 CellTagsIndexPosition[     43950,       1190]
+=======
+NotebookDataLength[     37111,       1018]
+NotebookOptionsPosition[     34220,        961]
+NotebookOutlinePosition[     34615,        977]
+CellTagsIndexPosition[     34572,        974]
+>>>>>>> 05aa02771abe37e27e9b199224b40ee7f45810f4
 WindowFrame->Normal*)
 
 (* Beginning of Notebook Content *)
@@ -86,7 +93,7 @@ Cell[BoxData[{
  RowBox[{
   RowBox[{"X0", " ", "=", " ", 
    RowBox[{"{", 
-    RowBox[{"0", ",", "y"}], "}"}]}], ";"}], "\[IndentingNewLine]", 
+    RowBox[{"x0", ",", "y"}], "}"}]}], ";"}], "\[IndentingNewLine]", 
  RowBox[{"X", " ", "=", " ", 
   RowBox[{"{", 
    RowBox[{"x", ",", "y"}], "}"}]}], "\[IndentingNewLine]", 
@@ -104,12 +111,9 @@ Cell[BoxData[{
         FractionBox["1", "\[Rho]1"], " ", 
         RowBox[{"D", "[", 
          RowBox[{
-          RowBox[{"uin", "[", "X", "]"}], ",", "x"}], "]"}]}], " ", "+", 
-       RowBox[{
-        FractionBox["1", "\[Rho]1"], 
-        RowBox[{"D", "[", " ", 
-         RowBox[{
-          RowBox[{"uR", "[", "X", "]"}], ",", "x"}], "]"}]}], " ", "-", 
+          RowBox[{
+           RowBox[{"uin", "[", "X", "]"}], "+", " ", 
+           RowBox[{"uR", "[", "X", "]"}]}], ",", "x"}], "]"}]}], " ", "-", 
        RowBox[{
         FractionBox["1", "\[Rho]2"], 
         RowBox[{"D", "[", 
@@ -139,43 +143,61 @@ Cell[BoxData[{
         RowBox[{"Thread", "[", 
          RowBox[{"eqs", "==", " ", "0"}], "]"}], "/.", "subSnells"}], ",", 
        RowBox[{"{", 
-        RowBox[{"R", ",", "T"}], "}"}]}], "]"}]}], "/.", " ", 
-    RowBox[{"k1", " ", "->", " ", 
+        RowBox[{"R", ",", "T"}], "}"}]}], "]"}]}], "/.", "  ", 
+    RowBox[{"\[Rho]2", "->", " ", 
      RowBox[{"Q", " ", "k2", " ", 
-      RowBox[{"\[Rho]1", "/", "\[Rho]2"}]}]}]}], " ", "//", 
+      RowBox[{"\[Rho]1", "/", "k1"}]}]}]}], " ", "//", 
    "Simplify"}]}]}], "Input",
  CellChangeTimes->{{3.906541843416709*^9, 3.906542195715724*^9}, {
   3.906542229477569*^9, 3.906542268783965*^9}, {3.906542314159116*^9, 
   3.906542336427779*^9}, {3.906542451667676*^9, 3.9065424769550533`*^9}, {
   3.906542728456843*^9, 3.906542728602861*^9}, {3.906542769766479*^9, 
-  3.906542771001995*^9}, {3.9163110611893473`*^9, 3.916311101609243*^9}},
- CellLabel->"In[54]:=",ExpressionUUID->"505ed892-b0ee-4b93-bf20-789175b2a272"],
+  3.906542771001995*^9}, {3.9163110611893473`*^9, 3.916311101609243*^9}, {
+  3.916390798484353*^9, 3.9163907988358917`*^9}, {3.916391202260542*^9, 
+  3.916391250295336*^9}},
+ CellLabel->"In[68]:=",ExpressionUUID->"505ed892-b0ee-4b93-bf20-789175b2a272"],
 
 Cell[BoxData[
  RowBox[{"{", 
   RowBox[{"x", ",", "y"}], "}"}]], "Output",
  CellChangeTimes->{{3.9065420170841618`*^9, 3.906542045835774*^9}, {
    3.906542135646899*^9, 3.906542196986534*^9}, 3.906542269001123*^9, 
    3.906542336756316*^9, 3.906542477180331*^9, 3.906542729035215*^9, 
-   3.906542772845778*^9, {3.916311047778927*^9, 3.916311102578001*^9}},
- CellLabel->"Out[62]=",ExpressionUUID->"fc2f5679-1c2b-49fa-882b-1731f835a4da"],
+   3.906542772845778*^9, {3.916311047778927*^9, 3.916311102578001*^9}, 
+   3.916390630140895*^9, 3.916390802680944*^9, {3.91639120552874*^9, 
+   3.9163912508145723`*^9}},
+ CellLabel->"Out[76]=",ExpressionUUID->"335fb237-18c0-43d1-ad5d-de3928c8561a"],
 
 Cell[BoxData[
  RowBox[{"{", 
   RowBox[{
    RowBox[{"R", "\[Rule]", 
     FractionBox[
      RowBox[{
-      RowBox[{"Q", " ", 
-       RowBox[{"Cos", "[", "\[Theta]1", "]"}]}], "-", 
-      RowBox[{"Cos", "[", "\[Theta]2", "]"}]}], 
+      SuperscriptBox["\[ExponentialE]", 
+       RowBox[{"2", " ", "\[ImaginaryI]", " ", "k1", " ", "x0", " ", 
+        RowBox[{"Cos", "[", "\[Theta]1", "]"}]}]], " ", 
+      RowBox[{"(", 
+       RowBox[{
+        RowBox[{"Q", " ", 
+         RowBox[{"Cos", "[", "\[Theta]1", "]"}]}], "-", 
+        RowBox[{"Cos", "[", "\[Theta]2", "]"}]}], ")"}]}], 
      RowBox[{
       RowBox[{"Q", " ", 
        RowBox[{"Cos", "[", "\[Theta]1", "]"}]}], "+", 
       RowBox[{"Cos", "[", "\[Theta]2", "]"}]}]]}], ",", 
    RowBox[{"T", "\[Rule]", 
     FractionBox[
-     RowBox[{"2", " ", "Q", " ", 
+     RowBox[{"2", " ", 
+      SuperscriptBox["\[ExponentialE]", 
+       RowBox[{"\[ImaginaryI]", " ", "x0", " ", 
+        RowBox[{"(", 
+         RowBox[{
+          RowBox[{"k1", " ", 
+           RowBox[{"Cos", "[", "\[Theta]1", "]"}]}], "-", 
+          RowBox[{"k2", " ", 
+           RowBox[{"Cos", "[", "\[Theta]2", "]"}]}]}], ")"}]}]], " ", "Q", 
+      " ", 
       RowBox[{"Cos", "[", "\[Theta]1", "]"}]}], 
      RowBox[{
       RowBox[{"Q", " ", 
@@ -184,8 +206,10 @@ Cell[BoxData[
  CellChangeTimes->{{3.9065420170841618`*^9, 3.906542045835774*^9}, {
    3.906542135646899*^9, 3.906542196986534*^9}, 3.906542269001123*^9, 
    3.906542336756316*^9, 3.906542477180331*^9, 3.906542729035215*^9, 
-   3.906542772845778*^9, {3.916311047778927*^9, 3.916311102600688*^9}},
- CellLabel->"Out[65]=",ExpressionUUID->"fbf24fc6-1eef-40cf-8680-2c99a27abd01"]
+   3.906542772845778*^9, {3.916311047778927*^9, 3.916311102578001*^9}, 
+   3.916390630140895*^9, 3.916390802680944*^9, {3.91639120552874*^9, 
+   3.9163912508973093`*^9}},
+ CellLabel->"Out[79]=",ExpressionUUID->"e5079cb5-0728-4534-9e5c-542711c2dcc0"]
 }, Open  ]],
 
 Cell[CellGroupData[{
@@ -544,13 +568,19 @@ Cell[BoxData[{
    3.916311027162562*^9}, {3.916311115609462*^9, 3.916311284904895*^9}, {
    3.9163113799386*^9, 3.916311394482524*^9}, {3.916311426640787*^9, 
    3.9163114270258093`*^9}, 3.916311474276861*^9, {3.916382402189849*^9, 
+<<<<<<< HEAD
    3.9163824234690104`*^9}, 3.9169793377431507`*^9},
  CellLabel->"In[23]:=",ExpressionUUID->"55751fef-0974-4177-ba5a-bc71fb3567d7"],
+=======
+   3.9163824234690104`*^9}},
+ CellLabel->"In[1]:=",ExpressionUUID->"55751fef-0974-4177-ba5a-bc71fb3567d7"],
+>>>>>>> 05aa02771abe37e27e9b199224b40ee7f45810f4
 
 Cell[BoxData[
  RowBox[{"{", 
   RowBox[{"x", ",", "y"}], "}"}]], "Output",
  CellChangeTimes->{3.9163114748345833`*^9, 3.916311852898614*^9, 
+<<<<<<< HEAD
   3.9163824296590137`*^9, 3.916979292269546*^9, 3.916979338036211*^9},
  CellLabel->"Out[36]=",ExpressionUUID->"d23a16fe-8428-427c-93eb-78158897f762"],
 
@@ -803,6 +833,10 @@ Cell[BoxData[
  CellChangeTimes->{3.9163114748345833`*^9, 3.916311852898614*^9, 
   3.9163824296590137`*^9, 3.916979292269546*^9, 3.9169793380771847`*^9},
  CellLabel->"Out[39]=",ExpressionUUID->"92b5dd1d-1291-41b8-b15d-6eada467d5ca"]
+=======
+  3.9163824296590137`*^9, 3.916390196741238*^9},
+ CellLabel->"Out[14]=",ExpressionUUID->"3563e61b-7090-41c9-a6f5-d334eb325176"]
+>>>>>>> 05aa02771abe37e27e9b199224b40ee7f45810f4
 }, Open  ]],
 
 Cell[CellGroupData[{
@@ -841,7 +875,11 @@ Cell[BoxData[
  CellChangeTimes->{{3.9163115735325537`*^9, 3.916311684079542*^9}, {
   3.9163819903119593`*^9, 3.9163820323579607`*^9}, {3.916382138712728*^9, 
   3.916382178342121*^9}, {3.916382455244184*^9, 3.916382455708098*^9}},
+<<<<<<< HEAD
  CellLabel->"In[40]:=",ExpressionUUID->"fa57b4ab-8700-426a-812d-7d22e665f9e6"],
+=======
+ CellLabel->"In[18]:=",ExpressionUUID->"fa57b4ab-8700-426a-812d-7d22e665f9e6"],
+>>>>>>> 05aa02771abe37e27e9b199224b40ee7f45810f4
 
 Cell[BoxData[
  FractionBox[
@@ -891,8 +929,13 @@ Cell[BoxData[
  CellChangeTimes->{{3.916311625302472*^9, 3.916311687922065*^9}, 
    3.91631186514664*^9, {3.9163820128719797`*^9, 3.916382032936818*^9}, {
    3.916382172674266*^9, 3.9163821786813383`*^9}, {3.9163824325172873`*^9, 
+<<<<<<< HEAD
    3.916382455928849*^9}, 3.9169792945538483`*^9, 3.916979359088944*^9},
  CellLabel->"Out[40]=",ExpressionUUID->"a40f6309-5623-4495-a241-bc245d1675eb"],
+=======
+   3.916382455928849*^9}, 3.916390201114299*^9},
+ CellLabel->"Out[18]=",ExpressionUUID->"18a0d1d7-f308-4a70-bb31-950eaec89494"],
+>>>>>>> 05aa02771abe37e27e9b199224b40ee7f45810f4
 
 Cell[BoxData[
  FractionBox[
@@ -918,8 +961,13 @@ Cell[BoxData[
  CellChangeTimes->{{3.916311625302472*^9, 3.916311687922065*^9}, 
    3.91631186514664*^9, {3.9163820128719797`*^9, 3.916382032936818*^9}, {
    3.916382172674266*^9, 3.9163821786813383`*^9}, {3.9163824325172873`*^9, 
+<<<<<<< HEAD
    3.916382455928849*^9}, 3.9169792945538483`*^9, 3.916979359090062*^9},
  CellLabel->"Out[41]=",ExpressionUUID->"dcbc02b8-c450-447c-88f1-42d77fdfeb7d"]
+=======
+   3.916382455928849*^9}, 3.916390201290867*^9},
+ CellLabel->"Out[19]=",ExpressionUUID->"5f51da47-faef-4571-a812-984913826dbe"]
+>>>>>>> 05aa02771abe37e27e9b199224b40ee7f45810f4
 }, Open  ]],
 
 Cell[CellGroupData[{
@@ -1175,9 +1223,9 @@ Cell[BoxData[
 }, Open  ]]
 }, Open  ]]
 },
-WindowSize->{1440, 762},
-WindowMargins->{{0, Automatic}, {0, Automatic}},
-FrontEndVersion->"13.3 for Linux x86 (64-bit) (June 3, 2023)",
+WindowSize->{1152., 605.4},
+WindowMargins->{{0, Automatic}, {Automatic, 0}},
+FrontEndVersion->"13.2 for Linux x86 (64-bit) (December 7, 2022)",
 StyleDefinitions->"Default.nb",
 ExpressionUUID->"ec8f41f7-53e2-40ba-8f40-226433cfd18b"
 ]
@@ -1193,26 +1241,27 @@ CellTagsIndex->{}
 (*NotebookFileOutline
 Notebook[{
 Cell[CellGroupData[{
-Cell[580, 22, 189, 3, 54, "Subsection",ExpressionUUID->"fc2d53ca-f1df-4305-9756-5a2ef579fb84"],
+Cell[580, 22, 189, 3, 64, "Subsection",ExpressionUUID->"fc2d53ca-f1df-4305-9756-5a2ef579fb84"],
 Cell[CellGroupData[{
-Cell[794, 29, 4141, 122, 362, "Input",ExpressionUUID->"505ed892-b0ee-4b93-bf20-789175b2a272"],
-Cell[4938, 153, 434, 7, 33, "Output",ExpressionUUID->"fc2f5679-1c2b-49fa-882b-1731f835a4da"],
-Cell[5375, 162, 1032, 25, 58, "Output",ExpressionUUID->"fbf24fc6-1eef-40cf-8680-2c99a27abd01"]
+Cell[794, 29, 4136, 121, 457, "Input",ExpressionUUID->"505ed892-b0ee-4b93-bf20-789175b2a272"],
+Cell[4933, 152, 533, 9, 40, "Output",ExpressionUUID->"335fb237-18c0-43d1-ad5d-de3928c8561a"],
+Cell[5469, 163, 1677, 41, 73, "Output",ExpressionUUID->"e5079cb5-0728-4534-9e5c-542711c2dcc0"]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[6444, 192, 1162, 30, 92, "Input",ExpressionUUID->"5fa6d3ee-515e-4d4f-b620-3b3d999753b5"],
-Cell[7609, 224, 279, 4, 33, "Output",ExpressionUUID->"3cb4d945-6e85-4b2c-9d5a-f43d2f9542d4"],
-Cell[7891, 230, 277, 4, 33, "Output",ExpressionUUID->"3bb59479-0d48-4621-9aaa-1626fcf721a3"]
+Cell[7183, 209, 1162, 30, 117, "Input",ExpressionUUID->"5fa6d3ee-515e-4d4f-b620-3b3d999753b5"],
+Cell[8348, 241, 279, 4, 40, "Output",ExpressionUUID->"3cb4d945-6e85-4b2c-9d5a-f43d2f9542d4"],
+Cell[8630, 247, 277, 4, 40, "Output",ExpressionUUID->"3bb59479-0d48-4621-9aaa-1626fcf721a3"]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[8205, 239, 297, 7, 71, "Subsubsection",ExpressionUUID->"d1cb2ce3-807f-4d39-b7c9-8c9054dbea11"],
+Cell[8944, 256, 297, 7, 85, "Subsubsection",ExpressionUUID->"d1cb2ce3-807f-4d39-b7c9-8c9054dbea11"],
 Cell[CellGroupData[{
-Cell[8527, 250, 1504, 43, 78, "Input",ExpressionUUID->"54e62f58-fdf6-4b56-9122-ce13eef202da"],
-Cell[10034, 295, 683, 18, 64, "Output",ExpressionUUID->"76976d20-2cfb-47b0-b25a-51817d64923f"]
+Cell[9266, 267, 1504, 43, 118, "Input",ExpressionUUID->"54e62f58-fdf6-4b56-9122-ce13eef202da"],
+Cell[10773, 312, 683, 18, 79, "Output",ExpressionUUID->"76976d20-2cfb-47b0-b25a-51817d64923f"]
 }, Open  ]],
-Cell[10732, 316, 384, 9, 58, "Input",ExpressionUUID->"0028ab2b-a57b-4f16-ae70-6bbe3d5fc8e7"],
-Cell[11119, 327, 1202, 30, 70, "Input",ExpressionUUID->"f518e691-2764-48d5-93d6-9b4a69e0aebd"],
+Cell[11471, 333, 384, 9, 72, "Input",ExpressionUUID->"0028ab2b-a57b-4f16-ae70-6bbe3d5fc8e7"],
+Cell[11858, 344, 1202, 30, 117, "Input",ExpressionUUID->"f518e691-2764-48d5-93d6-9b4a69e0aebd"],
 Cell[CellGroupData[{
+<<<<<<< HEAD
 Cell[12346, 361, 6335, 186, 621, "Input",ExpressionUUID->"55751fef-0974-4177-ba5a-bc71fb3567d7"],
 Cell[18684, 549, 289, 5, 33, "Output",ExpressionUUID->"d23a16fe-8428-427c-93eb-78158897f762"],
 Cell[18976, 556, 10858, 248, 248, "Output",ExpressionUUID->"92b5dd1d-1291-41b8-b15d-6eada467d5ca"]
@@ -1234,6 +1283,28 @@ Cell[39041, 1054, 1465, 35, 60, "Output",ExpressionUUID->"faf97621-e52d-46d6-8ec
 Cell[CellGroupData[{
 Cell[40543, 1094, 1837, 48, 121, "Input",ExpressionUUID->"1e0f5f1e-e6b5-4cb1-8547-d22264967da6"],
 Cell[42383, 1144, 1182, 28, 61, "Output",ExpressionUUID->"e811459d-cd1f-4826-97ac-1c2cb268777d"]
+=======
+Cell[13085, 378, 6343, 187, 784, "Input",ExpressionUUID->"55751fef-0974-4177-ba5a-bc71fb3567d7"],
+Cell[19431, 567, 267, 5, 40, "Output",ExpressionUUID->"3563e61b-7090-41c9-a6f5-d334eb325176"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[19735, 577, 1289, 34, 95, "Input",ExpressionUUID->"fa57b4ab-8700-426a-812d-7d22e665f9e6"],
+Cell[21027, 613, 2136, 49, 78, "Output",ExpressionUUID->"18a0d1d7-f308-4a70-bb31-950eaec89494"],
+Cell[23166, 664, 1018, 25, 94, "Output",ExpressionUUID->"5f51da47-faef-4571-a812-984913826dbe"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[24221, 694, 1275, 34, 71, "Input",ExpressionUUID->"ddfbfa82-f116-4c7e-a778-31889c380ac4"],
+Cell[25499, 730, 765, 23, 73, "Output",ExpressionUUID->"7fac7d07-c1fd-48bf-8214-dd25ff83690b"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[26301, 758, 875, 19, 66, "Input",ExpressionUUID->"2f57086d-3898-4086-82cd-a2df5438364d"],
+Cell[27179, 779, 2891, 68, 84, "Output",ExpressionUUID->"bc59f864-54e0-4a98-88f3-f19ae677cbc5"],
+Cell[30073, 849, 1048, 24, 71, "Output",ExpressionUUID->"6f3a18d5-a39c-450e-9089-7df0e80659f1"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[31158, 878, 1837, 48, 153, "Input",ExpressionUUID->"1e0f5f1e-e6b5-4cb1-8547-d22264967da6"],
+Cell[32998, 928, 1182, 28, 61, "Output",ExpressionUUID->"e811459d-cd1f-4826-97ac-1c2cb268777d"]
+>>>>>>> 05aa02771abe37e27e9b199224b40ee7f45810f4
 }, Open  ]]
 }, Open  ]]
 }, Open  ]]

diff --git a/src/EffectiveWaves.jl b/src/EffectiveWaves.jl
@@ -22,7 +22,7 @@ export  match_error, x_mesh_match
 ## effective waves
 export wavenumbers, wavenumber, asymptotic_monopole_wavenumbers, eigenvectors, convert_eigenvector_basis, eigenvector_length
 export WaveModes, WaveMode, wavemode_wienerhopf
-export solve_boundary_condition, scattering_field, material_scattering_coefficients, material_scattered_waves
+export solve_boundary_condition, scattering_field, material_scattering_coefficients, material_scattered_waves, material_effective_tmatrix
 
 export dispersion_equation, dispersion_complex, eigensystem # supplies a matrix used for the disperision equation and effective eignvectors
 

diff --git a/src/acoustics/discrete_wave/regular/material_scattering_coefficients.jl b/src/acoustics/discrete_wave/regular/material_scattering_coefficients.jl
@@ -74,3 +74,70 @@ function material_scattering_coefficients(scat_field::ScatteringCoefficientsFiel
     return [v]
 
 end
+
+function material_effective_tmatrix(ω::T, material::Material{Sphere{T,2}};
+    basis_order::Int = 2,
+    basis_order_contraction::Int = basis_order,
+    monopole_approximation::Bool = false
+    ) where T
+
+    # Extracting important parameters
+    medium = material.microstructure.medium
+    k = ω / medium.c
+    R = outer_radius(material.shape)
+    species = material.microstructure.species
+    Rtildas = R .- outer_radius.(species)
+
+    # Computing effective wavenumber
+    if basis_order_contraction < 2
+        k_eff = wavenumbers(ω, medium, species; basis_order = 2)[1]
+    else
+        k_eff = wavenumbers(ω, medium, species; basis_order = basis_order_contraction)[1]
+    end
+
+    # Solving artificial eigensystem
+    F = eigenvectors(ω, k_eff, material.microstructure, PlanarSymmetry{2}(); basis_order = basis_order_contraction)
+
+    _, n_λ, ps = size(F)
+
+    if ps > 1 @error "More than one eigenvector was found for a planar system! Using first eigenvector dfor calculations." end
+
+    L = basis_order_contraction+basis_order
+
+    n_densities = [number_density(sp) for sp in species]
+
+    J = [besselj(n,k*Rtildas[λ]) 
+        for n = -L-1:L, λ in 1:n_λ]
+
+    Jstar = [besselj(n,k_eff*Rtildas[λ])*n_densities[λ] 
+        for n = -L-1:L, λ in 1:n_λ]
+
+    H = [hankelh1(n,k*Rtildas[λ])
+        for n = -L-1:L, λ in 1:n_λ]
+
+    pre_num = k*J[1:1+2*L,:].*Jstar[2:2*(1+L),:] - k_eff*J[2:2*(1+L),:].*Jstar[1:1+2*L,:]
+    pre_denom = k*H[1:1+2*L,:].*Jstar[2:2*(1+L),:] - k_eff*H[2:2*(1+L),:].*Jstar[1:1+2*L,:]
+
+    Matrix_Num = complex(zeros(2*basis_order_contraction+1,2*basis_order+1))
+    Matrix_Denom = complex(zeros(2*basis_order_contraction+1,2*basis_order+1))
+    for n = 1:1+2*basis_order        
+        Matrix_Num[:,n] = vec(sum(pre_num[n+2*basis_order_contraction:-1:n,:].*F[:,:,1],dims=2))
+        Matrix_Denom[:,n] = vec(sum(pre_denom[n+2*basis_order_contraction:-1:n,:].*F[:,:,1],dims=2))
+    end
+
+    Tmats = - vec(sum(Matrix_Num,dims=1)./sum(Matrix_Denom,dims=1))
+
+    if monopole_approximation
+
+        TmatMs = complex(zeros(2basis_order+1))
+        for i in eachindex(TmatMs)
+            TmatMs[i] = -sum(pre_num[i+L-basis_order,:] .* F[basis_order_contraction+1,:,1])/sum(pre_denom[i+L-basis_order,:] .* F[basis_order_contraction+1,:,1])
+        end
+
+        return [Tmats, TmatMs]
+
+    else
+        return Tmats
+    end
+
+end