adapted tds2el harvesting to trigonals cells

ab2tds/TDS_Reading.py

@@ -109,8 +109,13 @@ class CalcDatas:
         cellVolume = CellVolume( self.cellVects )
         return cellVolume
         
-    def GetBrillVects(self):
-        Brillvectors= BrillVects(self.cellVects)
+    def GetBrillVects(self, CTransf=None):
+        if CTransf is None:
+            Brillvectors= BrillVects(self.cellVects)
+        else:
+            cellVects = numpy.dot(CTransf, self.cellVects)
+            Brillvectors= BrillVects(cellVects)
+            
         return Brillvectors
 
 

ab2tds/make_TDS_tds2elHarvest_Clone.py

-
 #/*##########################################################################
 # Copyright (C) 2011-2014 European Synchrotron Radiation Facility
 #
@@ -66,9 +65,10 @@
         * DeltaQ_min = 0.01
         * fitSpots = [[-2,-1,-2],[-1,0,4],[0,1,-8] ,[2,0,-2], [-2,-1,-2],  ]
 
-    * Interactions
+      * Cell transformation
+
+        * Ctransf       = [ [1,0,0],[0,1,0],[0,0,1]  ]
 
-        * right mouse click on graph plots intensity profile and displays position in reciprocal space
 
 The input variables are documented with docstrings below
   
@@ -149,6 +149,24 @@ DeltaQ_min = 0.01
 fitSpots = [[-2,-1,-2],[-1,0,4],[0,1,-8] ,[2,0,-2], [-2,-1,-2],  ]
 """
 """
+
+
+
+Ctransf       = None
+"""
+    if Ctransf is None, no transformation is done. This is equivalent to do the transformation [ [1,0,0],[0,1,0],[0,0,1]  ].
+    The transformation is done on the data, passing to a new unit cell UCp such that ::
+
+         UCp = Ctransf * UC
+
+    As an example ::
+
+      Ctransf = [ [0,1,-1],[1,-1,0],[1,1,1]  ]
+
+    transforms a trigonal UC into an hexagonal one
+
+"""
+
 if(sys.argv[0][-12:]!="sphinx-build"):
 
     s=open(sys.argv[2], "r").read()
@@ -219,8 +237,17 @@ if(sys.argv[0][-12:]!="sphinx-build"):
 
     BV=calculatedDatas.GetBrillVects()
     BVinv=numpy.linalg.inv(BV)
+
+    if Ctransf is not None:
+        BV_simple=calculatedDatas.GetBrillVects(Ctransf)
+        BVinv_simple=numpy.linalg.inv(BV_simple)
+    else:
+        BV_simple = BV
+        BVinv_simple = BVinv
+
     
-    maxBrill = numpy.max(numpy.sqrt(    (BV*BV).sum(axis=-1) ))
+    maxBrill = numpy.max(numpy.sqrt(    (BV_simple*BV_simple).sum(axis=-1) ))
+
     qradius_min = maxBrill * DeltaQ_min
     qradius_max = maxBrill * DeltaQ_max   ## in ab2tds i vettori Q hanno il 2pi
     
@@ -231,11 +258,16 @@ if(sys.argv[0][-12:]!="sphinx-build"):
     print "   DeltaQ_min  ,   DeltaQ_max  " , DeltaQ_min  ,   DeltaQ_max
     print  "  maxBrill  "  , maxBrill
     print  " norme.max(), norme.min()  ",  norme.max(), norme.min()
+    
+    reds = numpy.abs(numpy.tensordot(  Qs, BVinv_simple, axes = [[-1],[-1]]    ))
+    reds =  numpy.max(reds, axis=-1)
 
-    Qs = Qs[  numpy.less(qradius_min, norme)*numpy.less(norme, qradius_max )  ]
-    frequencies  = calculatedDatas.frequencies[numpy.less(qradius_min, norme)*numpy.less(norme, qradius_max )]
-    eigenvectors = calculatedDatas.eigenvectors[numpy.less(qradius_min, norme)*numpy.less(norme, qradius_max )]
+    Qs = Qs[  numpy.less(qradius_min, norme)*numpy.less(norme, qradius_max ) *numpy.less(reds, 0.5)  ]
 
+
+    frequencies  = calculatedDatas.frequencies[numpy.less(qradius_min, norme)*numpy.less(norme, qradius_max )*numpy.less(reds, 0.5)  ]
+    eigenvectors = calculatedDatas.eigenvectors[numpy.less(qradius_min, norme)*numpy.less(norme, qradius_max )*numpy.less(reds, 0.5) ]
+    
     dic_ints = {}
     hkl_max=0
     countspot=-1
@@ -250,10 +282,13 @@ if(sys.argv[0][-12:]!="sphinx-build"):
         # qsLine_notShifted = numpy.reshape(   Qs ,  [ (2*Nqs+1)* (2*Nqs+1), 3 ] )
         # intensity_line    =  numpy.reshape( intensity,  [ (2*Nqs+1)* (2*Nqs+1) ])
 
+
         qsLine_notShifted  =   Qs 
         intensity_line     =   intensity
-        qsLine_shifted  =   Qs  + numpy.dot(spot, BV)
 
+
+        qsLine_shifted  =   Qs  + numpy.dot(spot, BV_simple)
+        
         if NEUTRONCALC:
             scatterers = numpy.array([  CohB[aname] for aname in  calculatedDatas.atomNames ] , "F" )   
         else: