Refs #8116 Refactored FlatAbs

Code/Mantid/scripts/Inelastic/IndirectAbsCor.py

@@ -217,107 +217,107 @@ def AbsRunFeeder(inputWS, geom, beam, ncan, size, density, sigs, siga, avar,
 #  Output parameters :
 #  A1 - Ass ; A2 - Assc ; A3 - Acsc ; A4 - Acc
 
-def Fact(AMU,T,SEC1,SEC2):
-    S = AMU*T*(SEC1-SEC2)
+def Fact(xSection,thickness,sec1,sec2):
+    S = xSection*thickness*(sec1-sec2)
     F = 1.0
     if (S == 0.):
-        F = T
+        F = thickness
     else:
         S = (1-math.exp(-S))/S
-        F = T*S
+        F = thickness*S
     return F
 
 def FlatAbs(ncan, thick, density, sigs, siga, angles, waves):
+
+    #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
+
     PICONV = math.pi/180.
-    ssigs = sigs[0]                             #sam scatt x-sect
-    ssiga = siga[0]                             #sam abs x-sect
-    rhos = density[0]                           #sam density
-    TS = thick[0]                               #sam thicknes
-    T1 = thick[1]                               #can thickness 1
-    T2 = thick[2]                               #can thickness 2
-    csigs = sigs[1]                             #can scatt x-sect
-    csiga = siga[1]                             #can abs x-sect
-    rhoc = density[1]                           #can density
-    TCAN1 = angles[0]                           #angle can to beam
-    TCAN = TCAN1*PICONV
-    THETA1 = angles[1]                          #THETAB value - detector angle
-    THETA = PICONV*THETA1
 
-    AmuS1 = []                                # sample & can cross sections
-    AmuC1 = []
-    nlam = len(waves)
-    for n in range(0,nlam):
-       AS1 = ssigs + ssiga*waves[n]/1.8
-       AmuS1.append(AS1*rhos)
+    #can angle and detector angle
+    tcan1, theta1 = angles
+    canAngle = tcan1*PICONV
+    theta = PICONV*theta1
+
+    # tsec is the angle the scattered beam makes with the normal to the sample surface.
+    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):
-        for n in range(0,nlam):
-            AC1 = csigs + csiga*waves[n]/1.8
-            AmuC1.append(AC1*rhoc)
+        canXSection = (canScatt + canAbs * np.array(waves) /1.8) * canDensity
     else:
-        rhoc=0.
+        canDensity = 0.
+
+    nlam = len(waves)
+
+    ass = np.ones(nlam)
+    assc = np.ones(nlam)
+    acsc = np.ones(nlam)
+    acc = np.ones(nlam)
 
-    SEC1 = 1./math.cos(TCAN)
-    TSEC=THETA1-TCAN1    # TSEC IS THE ANGLE THE SCATTERED BEAM MAKES WITH THE NORMAL TO THE SAMPLE SURFACE.
-    A1 = []
-    A2 = []
-    A3 = []
-    A4 = []
-    if (abs(abs(TSEC)-90.0) < 1.0):            # case where TSEC is close to 90. CALCULATION IS UNRELIABLE
-        ASS = 1.0
-        for n in range(0,nlam):                #start loop over wavelengths
-            A1.append(ASS)
-            A2.append(ASS)
-            A3.append(ASS)
-            A4.append(ASS)
+    # 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))
     else:
-        TSEC = TSEC*PICONV
-        SEC2 = 1./math.cos(TSEC)
-        for n in range(0,nlam):                   #start loop over wavelengths
-            AMUS = AmuS1[n]
-            FS = Fact(AMUS,TS,SEC1,SEC2)
-            ES1=AMUS*TS*SEC1
-            ES2=AMUS*TS*SEC2
-            if (ncan > 1):
-                AMUC = AmuC1[n]
-                F1 = Fact(AMUC,T1,SEC1,SEC2)
-                F2 = Fact(AMUC,T2,SEC1,SEC2)
-                E11 = AMUC*T1*SEC1
-                E12 = AMUC*T1*SEC2
-                E21 = AMUC*T2*SEC1
-                E22 = AMUC*T2*SEC2
-            if (SEC2 < 0.):
-                ASS=FS/TS
-                if(ncan > 1):
-                    ASSC = ASS*math.exp(-(E11-E12))
-                    ACC1 = F1
-                    ACC2 = F2*math.exp(-(E11-E12))
-                    ACSC1 = ACC1
-                    ACSC2 = ACC2*math.exp(-(ES1-ES2))
-                else:
-                    ASSC = 1.0
-                    ACSC = 1.0
-                    ACC = 1.0
-            else:
-                ASS=math.exp(-ES2)*FS/TS
-                if(ncan > 1):
-                    ASSC = math.exp(-(E11+E22))*ASS
-                    ACC1 = math.exp(-(E12+E22))*F1
-                    ACC2 = math.exp(-(E11+E22))*F2
-                    ACSC1 = ACC1*math.exp(-ES2)
-                    ACSC2 = ACC2*math.exp(-ES1)
-                else:
-                    ASSC = 1.0
-                    ACSC = 1.0
-                    ACC = 1.0
-            tsum = T1+T2
-            if(tsum > 0.):
-                ACC = (ACC1+ACC2)/tsum
-                ACSC = (ACSC1+ACSC2)/tsum
-            else:
-                ACC = 1.0
-                ACSC = 1.0
-            A1.append(ASS)
-            A2.append(ASSC)
-            A3.append(ACSC)
-            A4.append(ACC)
-	return A1, A2, A3, A4
+
+        fs = vecFact(sampleXSection, samThickness, sec1, sec2)
+
+        sampleSec1 = sampleXSection * samThickness * sec1
+        sampleSec2 = sampleXSection * samThickness * sec2
+
+        if ( ncan > 1 ):
+            f1 = vecFact(canXSection,canThickness1,sec1,sec2)
+            f2 = vecFact(canXSection,canThickness2,sec1,sec2)
+
+            canThick1Sec1 = canXSection * canThickness1 * sec1
+            canThick1Sec2 = canXSection * canThickness1 * sec2
+            canThick2Sec1 = canXSection * canThickness2 * sec1
+            canThick2Sec2 = canXSection * canThickness2 * sec2
+
+        if (sec2 < 0.):
+            ass = fs / samThickness
+
+            if( ncan > 1 ):
+                val = np.exp(-(canThick1Sec1-canThick1Sec2))
+                assc = ass * val
+
+                acc1 = f1
+                acc2 = f2 * val
+
+                acsc1 = acc1
+                acsc2 = acc2 * np.exp(-(sampleSec1-sampleSec2))
+        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
+
+    return ass, assc, acsc, acc