From 74eec48ca942ad1fb90f1726e78118345b69313f Mon Sep 17 00:00:00 2001
From: pbienst <pbienst@78e232d7-6d57-4136-945a-988d5fb676e4>
Date: Tue, 4 Mar 2008 10:41:22 +0000
Subject: [PATCH] Better handling of near singular matrices in GARCLED.py

---
 ChangeLog        |  2 ++
 camfr/GARCLED.py | 34 +++++++++++++++++++++++++++++++---
 2 files changed, 33 insertions(+), 3 deletions(-)

diff --git a/ChangeLog b/ChangeLog
index 0ab8860a..68dc355f 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,5 @@
+ - Better handling of near-singular matrices in GARCLED.py.
+
 CAMFR 1.3:
 
  - Added code to model spontaneous emission and light extraction in light
diff --git a/camfr/GARCLED.py b/camfr/GARCLED.py
index df0da38f..8eb5cb05 100644
--- a/camfr/GARCLED.py
+++ b/camfr/GARCLED.py
@@ -270,6 +270,34 @@ def horizontal_y(cav, kx0, ky0):
 
 
 
+#############################################################################
+#
+# Safe matrix inversion, uses svd if matrix is singular.
+#
+#############################################################################
+
+def inv_safe(A):
+
+  try:
+
+    return inv(A)
+
+  except:
+
+    print "Singularities dectected in matrix inversion."
+
+    # A = U * sigma * Vh
+
+    U, sigma, Vh = linalg.svd(A)
+
+    # inv_A = V * 1/sigma * Uh, but set 1/0 in sigma to 0.
+
+    inv_sigma = where(abs(sigma) <= 1e-11, 0, 1/sigma) 
+    
+    return dot(Vh.conj().T, dot(diag(inv_sigma), U.conj().T))
+
+
+
 #############################################################################
 #
 # Calculate extracted field through a cavity, but only for waves
@@ -306,8 +334,8 @@ def P(cav, sources, kx0, ky0, single_pass_substrate=True, add_all=True,
 
   N = cav.top.inc().N()
 
-  inv_1 = inv(identity(N) - dot(R_bot, R_top)) 
-  inv_2 = inv(identity(N) - dot(R_top, R_bot))
+  inv_1 = inv_safe(identity(N) - dot(R_bot, R_top)) 
+  inv_2 = inv_safe(identity(N) - dot(R_top, R_bot))
 
   # Loop over the different sources.
 
@@ -398,7 +426,7 @@ def P(cav, sources, kx0, ky0, single_pass_substrate=True, add_all=True,
       # Calculate transmitted power density.
         
       if len(propagating) > 0: # LAPACK crashes on empty matrix.
-        inv_  = inv(identity(len(propagating)) - dot(R_cav, R_sub))
+        inv_  = inv_safe(identity(len(propagating)) - dot(R_cav, R_sub))
         P_out = dot(T_sub, dot(inv_, P_sub_))
       else:
         P_out = [0.0]