Merge branch 'feature/8116_fltabs_code_tidy' into develop

Conflicts:
	Code/Mantid/scripts/Inelastic/IndirectAbsCor.py

Refs #8116

Code/Mantid/scripts/Inelastic/IndirectAbsCor.py

@@ -204,121 +204,142 @@ def AbsRunFeeder(inputWS, geom, beam, ncan, size, density, sigs, siga, avar,
         graph.activeLayer().setAxisTitle(mp.Layer.Bottom, 'Angle')
 
 
-# FlatAbs - CALCULATE FLAT PLATE ABSORPTION FACTORS
+# FlatAbs - calculate flat plate absorption factors
+# 
+# For more information See:
+#   - MODES User Guide: http://www.isis.stfc.ac.uk/instruments/iris/data-analysis/modes-v3-user-guide-6962.pdf  
+#   - C J Carlile, Rutherford Laboratory report, RL-74-103 (1974)  
 #
 #  Input parameters :
 #  sigs - list of scattering  cross-sections
 #  siga - list of absorption cross-sections
 #  density - list of density
-#  ncan - =0 no can, >1 with can
-#  thick - list of thicknesses ts,t1,t2
+#  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(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 calcThicknessAtSec(xSection, thickness, sec):
+    sec1, sec2 = sec
+
+    thickSec1 = xSection * thickness * sec1
+    thickSec2 = xSection * thickness * sec2
+
+    return thickSec1, thickSec2
+
+def calcFlatAbsCan(ass, canXSection, canThickness1, canThickness2, sampleSec1, sampleSec2, sec):
+    nlam = len(canXSection)
+
+    assc = np.ones(nlam)
+    acsc = np.ones(nlam)
+    acc = np.ones(nlam)
+
+    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
+
+    return assc, acsc, acc
+
 def FlatAbs(ncan, thick, density, sigs, siga, angles, waves):
+
     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 = []
+
+    #can angle and detector angle
+    tcan1, theta1 = angles
+    canAngle = tcan1*PICONV
+    theta = theta1*PICONV
+
+    # tsec is the angle the scattered beam makes with the normal to the sample surface.
+    tsec = theta1-tcan1
+
     nlam = len(waves)
-    for n in range(0,nlam):
-       AS1 = ssigs + ssiga*waves[n]/1.8
-       AmuS1.append(AS1*rhos)
-    if (ncan > 1):
-        for n in range(0,nlam):
-            AC1 = csigs + csiga*waves[n]/1.8
-            AmuC1.append(AC1*rhoc)
-    else:
-        rhoc=0.
-
-    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)
+
+    ass = np.ones(nlam)
+    assc = np.ones(nlam)
+    acsc = np.ones(nlam)
+    acc = np.ones(nlam)
+
+    # 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 ass, assc, acsc, acc
     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
+        #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
+
+        tsec = tsec*PICONV
+
+        sec1 = 1./math.cos(canAngle)
+        sec2 = 1./math.cos(tsec)
+
+        #list of wavelengths
+        waves = np.array(waves)
+
+        #sample cross section
+        sampleXSection = (sampleScatt + sampleAbs * waves /1.8) * sampleDensity
+
+        #vector version of fact
+        vecFact = np.vectorize(Fact)
+        fs = vecFact(sampleXSection, samThickness, sec1, sec2)
+
+        sampleSec1, sampleSec2 = calcThicknessAtSec(sampleXSection, samThickness, [sec1, sec2])
+
+        if (sec2 < 0.):
+            ass = fs / samThickness
+        else:
+            ass= np.exp(-sampleSec2) * fs / samThickness
+
+        useCan = (ncan > 1)
+        if useCan:
+            #calculate can cross section
+            canXSection = (canScatt + canAbs * waves /1.8) * canDensity
+            assc, acsc, acc = calcFlatAbsCan(ass, canXSection, canThickness1, canThickness2, sampleSec1, sampleSec2, [sec1, sec2])
+
+    return ass, assc, acsc, acc