Refs #8116 Refactor code into functions

Code/Mantid/scripts/Inelastic/IndirectAbsCor.py

@@ -214,12 +214,10 @@ def AbsRunFeeder(inputWS, geom, beam, ncan, size, density, sigs, siga, avar,
 #  sigs - list of scattering  cross-sections
 #  siga - list of absorption cross-sections
 #  density - list of density
-#  ncan - =0 no can, >1 with can
+#  ncan - = 0 no can, >1 with can
 #  thick - list of thicknesses: sample thickness, can thickness1, can thickness2
 #  angles - list of angles
 #  waves - list of wavelengths
-#  Output parameters :
-#  A1 - Ass ; A2 - Assc ; A3 - Acsc ; A4 - Acc
 
 def Fact(xSection,thickness,sec1,sec2):
     S = xSection*thickness*(sec1-sec2)
@@ -231,16 +229,60 @@ def Fact(xSection,thickness,sec1,sec2):
         F = thickness*S
     return F
 
-def FlatAbs(ncan, thick, density, sigs, siga, angles, waves):
+def calcThicknessAtSec(xSection, thickness, sec):
+    sec1, sec2 = sec
+
+    thickSec1 = xSection * thickness * sec1
+    thickSec2 = xSection * thickness * sec2
+
+    return thickSec1, thickSec2
+
+def calcFlatAbsCan(ass, canXSection, thick, sec):
+    assc = np.ones(nlam)
+    acsc = np.ones(nlam)
+    acc = np.ones(nlam)
+
+    #thickness of the can
+    canThickness1, canThickness2 = thick
+
+    sec1, sec2 = sec
+
+    #vector version of fact
+    vecFact = np.vectorize(Fact)
+    f1 = vecFact(canXSection,canThickness1,sec1,sec2)
+    f2 = vecFact(canXSection,canThickness2,sec1,sec2)
+
+    canThick1Sec1, canThick1Sec2 = calcThicknessAtSec(canXSection, canThickness1, sec)
+    canThick2Sec1, canThick2Sec2 = calcThicknessAtSec(canXSection, canThickness2, sec)
+
+    if (sec2 < 0.):
+        val = np.exp(-(canThick1Sec1-canThick1Sec2))
+        assc = ass * val
+
+        acc1 = f1
+        acc2 = f2 * val
+
+        acsc1 = acc1
+        acsc2 = acc2 * np.exp(-(sampleSec1-sampleSec2))
+    else:
+        val = np.exp(-(canThick1Sec1+canThick2Sec2))
+        assc = ass * val
+
+        acc1 = f1 * np.exp(-(canThick1Sec2+canThick2Sec2))
+        acc2 = f2 * val
+
+        acsc1 = acc1 * np.exp(-sampleSec2)
+        acsc2 = acc2 * np.exp(-sampleSec1)
+
+    canThickness = canThickness1+canThickness2
+
+    if(canThickness > 0.):
+        acc = (acc1+acc2)/canThickness
+        acsc = (acsc1+acsc2)/canThickness
 
-    #sample & can scattering x-section
-    sampleScatt, canScatt = sigs[:2]
-    #sample & can absorption x-section                           
-    sampleAbs, canAbs = siga[:2]
-    #sample & can density                           
-    sampleDensity, canDensity = density[:2]
-    #thickness of the sample and can
-    samThickness, canThickness1, canThickness2 = thick
+    return assc, acsc, acc
+
+def FlatAbs(ncan, thick, density, sigs, siga, angles, waves):
 
     PICONV = math.pi/180.
 
@@ -253,20 +295,6 @@ def FlatAbs(ncan, thick, density, sigs, siga, angles, waves):
     tsec = theta1-tcan1
     tsec = tsec*PICONV
 
-    sec1 = 1./math.cos(canAngle)
-    sec2 = 1./math.cos(tsec)
-
-    #vector version of fact
-    vecFact = np.vectorize(Fact)
-
-    #sample & can cross sections
-    sampleXSection = (sampleScatt + sampleAbs * np.array(waves) /1.8) * sampleDensity
-
-    if (ncan > 1):
-        canXSection = (canScatt + canAbs * np.array(waves) /1.8) * canDensity
-    else:
-        canDensity = 0.
-
     nlam = len(waves)
 
     ass = np.ones(nlam)
@@ -277,51 +305,41 @@ def FlatAbs(ncan, thick, density, sigs, siga, angles, waves):
     # case where tsec is close to 90 degrees. CALCULATION IS UNRELIABLE
     if (abs(abs(tsec)-90.0) < 1.0):
         #default to 1 for everything
-        return np.ones((4,nlam))
+        return ass, assc, acsc, acc
     else:
-
-        fs = vecFact(sampleXSection, samThickness, sec1, sec2)
+        #sample & can scattering x-section
+        sampleScatt, canScatt = sigs[:2]
+        #sample & can absorption x-section                           
+        sampleAbs, canAbs = siga[:2]
+        #sample & can density                           
+        sampleDensity, canDensity = density[:2]
+        #thickness of the sample and can
+        samThickness, canThickness1, canThickness2 = thick
 
-        sampleSec1 = sampleXSection * samThickness * sec1
-        sampleSec2 = sampleXSection * samThickness * sec2
+        sec1 = 1./math.cos(canAngle)
+        sec2 = 1./math.cos(tsec)
 
-        if ( ncan > 1 ):
-            f1 = vecFact(canXSection,canThickness1,sec1,sec2)
-            f2 = vecFact(canXSection,canThickness2,sec1,sec2)
+        #list of wavelengths
+        waves = np.array(waves)
 
-            canThick1Sec1 = canXSection * canThickness1 * sec1
-            canThick1Sec2 = canXSection * canThickness1 * sec2
-            canThick2Sec1 = canXSection * canThickness2 * sec1
-            canThick2Sec2 = canXSection * canThickness2 * sec2
+        #sample cross section
+        sampleXSection = (sampleScatt + sampleAbs * waves /1.8) * sampleDensity
 
-        if (sec2 < 0.):
-            ass = fs / samThickness
-
-            if( ncan > 1 ):
-                val = np.exp(-(canThick1Sec1-canThick1Sec2))
-                assc = ass * val
+        #vector version of fact
+        vecFact = np.vectorize(Fact)
+        fs = vecFact(sampleXSection, samThickness, sec1, sec2)
 
-                acc1 = f1
-                acc2 = f2 * val
+        sampleSec1, sampleSec2 = calcThicknessAtSec(sampleXSection, samThickness, [sec1, sec2])
 
-                acsc1 = acc1
-                acsc2 = acc2 * np.exp(-(sampleSec1-sampleSec2))
+        if (sec2 < 0.):
+            ass = fs / samThickness
         else:
             ass= np.exp(-sampleSec2) * fs / samThickness
 
-            if(ncan > 1):
-                val = np.exp(-(canThick1Sec1+canThick2Sec2))
-                assc = ass * val
-
-                acc1 = f1 * np.exp(-(canThick1Sec2+canThick2Sec2))
-                acc2 = f2 * val
-
-                acsc1 = acc1 * np.exp(-sampleSec2)
-                acsc2 = acc2 * np.exp(-sampleSec1)
-
-        canThickness = canThickness1+canThickness2
-        if(canThickness > 0.):
-            acc = (acc1+acc2)/canThickness
-            acsc = (acsc1+acsc2)/canThickness
+        useCan = (ncan > 1)
+        if useCan:
+            #calculate can cross section
+            canXSection = (canScatt + canAbs * waves /1.8) * canDensity
+            assc, acsc, acc = calcFlatAbsCan(ass, canXSection, [canThickness1, canThickness2], [sec1, sec2])
 
     return ass, assc, acsc, acc