From 8881254f960c83715fcdf4c3de65094e9156968f Mon Sep 17 00:00:00 2001
From: Christopher Wipf <wipf@ligo.mit.edu>
Date: Thu, 16 Jul 2020 09:36:41 -0700
Subject: [PATCH] Refactor TF calculation and swap upper and lower joints

The joint order in hForce/vForce must remain consistent between
this code and the thermal noise calculation

Also remove a couple of unused variables

Fixes #68
---
 gwinc/suspension.py | 23 ++++++++---------------
 1 file changed, 8 insertions(+), 15 deletions(-)

diff --git a/gwinc/suspension.py b/gwinc/suspension.py
index 984f868b..3d60edff 100644
--- a/gwinc/suspension.py
+++ b/gwinc/suspension.py
@@ -420,11 +420,9 @@ def suspQuad(f, sus):
     # Complex spring constants
     khr = np.zeros([len(stages), len(w)])
     khi = np.zeros([len(stages), 2, len(w)])
-    kh  = np.zeros([len(stages), len(w)], dtype=np.complex_)
     kvr = np.zeros([len(stages), len(w)])
     kvi = np.zeros([len(stages), 2, len(w)])
-    kv  = np.zeros([len(stages), len(w)], dtype=np.complex_)
-
+    
     for n, stage in enumerate(stages):
 
         ##############################
@@ -570,25 +568,20 @@ def suspQuad(f, sus):
 
     for m in range(2*len(stages)):
         # turn on only the loss of the current joint
-        lossy_region_lower = np.zeros((len(stages), 1))
-        lossy_region_upper = np.zeros((len(stages), 1))
         n = int(m/2)  # stage number
-        if m % 2:
-            lossy_region_upper[n] = 1
-        else:
-            lossy_region_lower[n] = 1
+        isLower = m % 2
+        stage_selection = np.zeros((len(stages), 1))
+        stage_selection[n] = 1
 
         # horizontal
-        # only the imaginary part due to the specified region is used.
-        k = khr + 1j*(khi[:,0,:]*lossy_region_upper +
-                      khi[:,1,:]*lossy_region_lower)
+        # only the imaginary part due to the specified joint is used.
+        k = khr + 1j*khi[:,isLower,:]*stage_selection
         # calculate TFs
         hForce[m,:] = tst_force_to_tst_displ(k, masses, f)
 
         # vertical
-        # only the imaginary part due to the specified stage is used
-        k = kvr + 1j*(kvi[:,0,:]*lossy_region_upper +
-                      kvi[:,1,:]*lossy_region_lower)
+        # only the imaginary part due to the specified joint is used
+        k = kvr + 1j*kvi[:,isLower,:]*stage_selection
         # calculate TFs
         vForce[m,:] = tst_force_to_tst_displ(k, masses, f)
 
-- 
GitLab