Changeset 3896
 Timestamp:
 Apr 14, 2019 9:23:46 AM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImath.py
r3894 r3896 3657 3657 :returns: density at xyz 3658 3658 ''' 3659 def getBoxes(rho,I): 3660 Rhos = np.zeros((2,2,2)) 3661 Mx,My,Mz = rho.shape 3662 Ix,Iy,Iz = I 3663 Rhos = np.array([[[rho[Ix%Mx,Iy%My,Iz%Mz],rho[Ix%Mx,Iy%My,(Iz+1)%Mz]], 3664 [rho[Ix%Mx,(Iy+1)%My,Iz%Mz],rho[Ix%Mx,(Iy+1)%My,(Iz+1)%Mz]]], 3665 [[rho[(Ix+1)%Mx,Iy%My,Iz%Mz],rho[(Ix+1)%Mx,Iy%My,(Iz+1)%Mz]], 3666 [rho[(Ix+1)%Mx,(Iy+1)%My,Iz%Mz],rho[(Ix+1)%Mx,(Iy+1)%My,(Iz+1)%Mz]]]]) 3667 return Rhos 3668 3659 3669 Blk = 8 #8 seems optimal 3660 3670 nBlk = len(XYZ)//Blk #select Blk so this is an exact divide 3661 rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)[:2,:2,:2]3662 3671 mapShape = np.array(rho.shape) 3663 3672 mapStep = 1./mapShape … … 3665 3674 iBeg = 0 3666 3675 R = np.zeros(len(XYZ)) 3667 #actually a bit faster!3676 #actually a lot faster! 3668 3677 for iblk in range(nBlk): 3669 3678 iFin = iBeg+Blk 3670 3679 Xs = X[iBeg:iFin] 3671 3680 I = np.array(np.rint(Xs*mapShape),dtype='int') 3672 Rhos = np.array([ rollMap(rho,i) for i in I])3681 Rhos = np.array([getBoxes(rho,i) for i in I]) 3673 3682 Ds = XsI*mapStep 3674 3683 RIJs = Rhos[:,0,:2,:2]*(1.Ds[:,0][:,nxs,nxs]) … … 3676 3685 R[iBeg:iFin] = RIs[:,0]*(1.Ds[:,2])+RIs[:,1]*Ds[:,2] 3677 3686 iBeg += Blk 3678 #one at a time3679 # for i,x in enumerate(X):3680 # D = xI[i]*mapStep #position inside map cell3681 # Rho = rollMap(rho,I[i]) #shifts map so point is in corner3682 # RIJ = Rho[0,:2,:2]*(1.D[0])3683 # RI = RIJ[0]*(1.D[1])+RIJ[1]*D[1]3684 # R[i] = RI[0]*(1.D[2])+RI[1]*D[2]3685 3687 return R 3686 3688
Note: See TracChangeset
for help on using the changeset viewer.