scrivendo contributo secondo ordine : header della routine

ab2tds/make_TDS_tds2elHarvest_Clone.py

@@ -74,6 +74,7 @@ The input variables are documented with docstrings below
   
 """ 
 NOREPLICA=False
+DOSECONDORDER = True
 
 import startmessage
 
@@ -289,10 +290,8 @@ if(sys.argv[0][-12:]!="sphinx-build"):
 
     if NOREPLICA:
         Qs =     calculatedDatas.GetNonClippedQs()
-
         frequencies  = calculatedDatas.frequencies
         eigenvectors = calculatedDatas.eigenvectors
-        
     else:
         res = TDS_Simmetry.GetAllReplicated (  calculatedDatas, simmetries_dict ,  filename=filename , Nfour_interp=None,
                         md5postfix=md5postfix,
@@ -327,8 +326,7 @@ if(sys.argv[0][-12:]!="sphinx-build"):
     reds = numpy.abs(numpy.tensordot(  Qs, BVinv_simple, axes = [[-1],[-1]]    ))
     reds =  numpy.max(reds, axis=-1)
 
-
-    lut = LUT(Qs ) 
+    lut = LUT(calculatedDatas.QsReduce,frequencies ,eigenvectors) 
     
     Qs_roi           = Qs[  numpy.less(qradius_min, norme)*numpy.less(norme, qradius_max ) *numpy.less(reds, 0.5)  ]
     frequencies_roi  = frequencies[numpy.less(qradius_min, norme)*numpy.less(norme, qradius_max )*numpy.less(reds, 0.5)  ]
@@ -434,52 +432,39 @@ if(sys.argv[0][-12:]!="sphinx-build"):
                         dwfacts_mass_kin_scatt_perq
                         ):
           count+=1
-
           # if count%nprocs!=myrank:
           #   continue
-
-          # dm = TDS_Simmetry.GetDM_fromF(q,  Replica["Four_Rvects"],   Replica["Four_DMs"]  )
-          # evals, evects = numpy.linalg.eigh(dm)
-
           evals = frequencies_roi[count]
           evects= eigenvectors_roi[count]
-
           # originally the procedure beklow which has been copied,
           # was getting vects obtained from numpy.linalg.eig
           # which gives vectors as columns
           evects=numpy.transpose(evects)
-
           evals=(numpy.maximum(evals, numpy.array([bottom_meV/(0.0001239852 *1000)], numpy.float32) ))
-
           qvect_s = numpy.tensordot(  qshifted ,      evects.reshape([-1,3, Nbranches]) , [[0], [1]   ]  )
-
           Fjs = numpy.sum(dwfacts_mass_kin_scatt* (qvect_s.T), axis=-1)  # sommato sugli ioni
           Fjs_debug = numpy.sum(debugfacts[count] * (qvect_s.T), axis=-1)  # sommato sugli ioni
-
           exp_plus =  numpy.exp( factor_forcoth * evals)
           exp_minus = numpy.exp(-factor_forcoth * evals)
-          # print evals[:3], factor_forcoth ,  evals[:3]* factor_forcoth 
-          # print "MIN MAX ", exp_plus[0], exp_plus[2]
-
           intensity_array = ( Fjs*(Fjs.conjugate())).real
           intensity_array_debug = ( Fjs_debug*(Fjs_debug.conjugate())).real
-
           if NEUTRONCALC:
               kinstokes = numpy.sqrt(numpy.maximum( (NeutronE-evals*(0.0001239852 *1000))/NeutronE, 0.0 ))
               kinantistk= numpy.sqrt((NeutronE+evals*(0.0001239852 *1000))/NeutronE)
           else:
               kinstokes = 1.0
               kinantistk= 1.0
-
           deno = (exp_plus-exp_minus)*evals
-
           totalstats = ( exp_plus*kinstokes+exp_minus*kinantistk)/deno
-          # totalstats[3:] = 0 
-
-          # intensity_line[count] = numpy.sum( intensity_array_debug*  totalstats      )
           intensity_line[count] = numpy.sum( intensity_array* totalstats      )
-          # intensity_line[count] = evals.min()
-          # intensity_line[count] = fake_int[count]
+          if DOSECONDORDER:
+
+              red = numpy.dot( BVinv , qnotshifted )
+              add = secondorder(qnotshifted,red,   qshifted, Qs, calculatedDatas.QsReduced, 
+                                lut , massfactors, dwfacts_mass_kin_scatt,factor_forcoth  )
+
+          
+
 
           if countspot == 0 and count == -21971 : 
               file.write("count=%d\n"%count ) 
@@ -498,7 +483,6 @@ if(sys.argv[0][-12:]!="sphinx-build"):
               file.write("tmp**2\n") 
               file.write(str(intensity_array_debug[:3])+"\n") 
               file.close()
-        
         hkl_max = max( hkl_max, (numpy.abs( numpy.array( spot )   )).max())
         dic_ints[tuple(spot)] = intensity_line