Changeset 2073 for trunk

Timestamp:

Nov 30, 2015 8:33:57 AM (10 years ago)

Author:

vondreele

Message:

numeric derivatives for ZigZag/Block? waves

File:

• trunk/GSASIImath.py (modified) (5 diffs)

trunk/GSASIImath.py ¶

@@ -962 +962 @@
             XmodZ[iatm][0] += posZigZag(glTau,Tmm,XYZmax).T
         elif 'Block' in waveTypes[iatm]:
-
-            nx = 1
             Tmm = Ax[iatm][0][:2]                        
             XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
             XmodZ[iatm][0] += posBlock(glTau,Tmm,XYZmax).T
         tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
-        XmodA[iatm][nx:] = Ax[iatm,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X 32
-        XmodB[iatm][nx:] = Bx[iatm,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
+        if nx:    
+            XmodA[iatm][:-nx] = Ax[iatm,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X 32
+            XmodB[iatm][:-nx] = Bx[iatm,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
+        else:
+            XmodA[iatm] = Ax[iatm,:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X 32
+            XmodB[iatm] = Bx[iatm,:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
     Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X 3 X 32; sum waves
     Xmod = np.swapaxes(Xmod,1,2)
@@ -1019 +1021 @@
     '''
     glTau,glWt = pwd.pygauleg(0.,1.,32)         #get Gauss-Legendre intervals & weights
+    dT = 1/32.
+    dX = 0.0001
     waveShapes = [FSSdata.T.shape,XSSdata.T.shape,USSdata.T.shape]
     Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
@@ -1036 +1040 @@
         if 'ZigZag' in waveTypes[iatm]:
             nx = 1
-#            Tmm = Ax[iatm][0][:2]                        
-#            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
-#            DT = Tmm[1]-Tmm[0]
-#            slopeUp = 2.*XYZmax/DT
-#            slopeDn = 2.*XYZmax/(1.-DT)
-#            XmodZ[iatm][0] += np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax+slopeUp*((t-Tmm[0])%1.),XYZmax-slopeDn*((t-Tmm[1])%1.)) for t in glTau]).T
+            Tmm = Ax[iatm][0][:2]                        
+            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])            
+            for i in range(2):
+                Tmm[i] += dT
+                Zp = posZigZag(glTau,Tmm,XYZmax).T
+                Tmm[i] -= 2.*dT
+                Zm = posZigZag(glTau,Tmm,XYZmax).T
+                Tmm[i] += dT
+                ZtauXt[iatm][i] = (Zp-Zm)[i]/(2.*dT)
+            for i in range(3):
+                XYZmax[i] += dX
+                Zp = posZigZag(glTau,Tmm,XYZmax).T
+                XYZmax[i] -= 2.*dX
+                Zm = posZigZag(glTau,Tmm,XYZmax).T
+                XYZmax[i] += dX
+                ZtauXx[iatm][i] = (Zp-Zm)[i]/(2.*dX)
         elif 'Block' in waveTypes[iatm]:
             nx = 1
             Tmm = Ax[iatm][0][:2]                        
-            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
-            ZtauXt[iatm] = np.ones(2)[:,nxs]
-            ZtauXx[iatm] += np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-np.ones(3),np.ones(3)) for t in glTau]).T                    
+            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])            
+            for i in range(2):
+                Tmm[i] += dT
+                Zp = posBlock(glTau,Tmm,XYZmax).T
+                Tmm[i] -= 2.*dT
+                Zm = posBlock(glTau,Tmm,XYZmax).T
+                Tmm[i] += dT
+                ZtauXt[iatm][i] = (Zp-Zm)[i]/(2.*dT)
+            for i in range(3):
+                XYZmax[i] += dX
+                Zp = posBlock(glTau,Tmm,XYZmax).T
+                XYZmax[i] -= 2.*dX
+                Zm = posBlock(glTau,Tmm,XYZmax).T
+                XYZmax[i] += dX
+                ZtauXx[iatm][i] = (Zp-Zm)[i]/(2.*dX)
         tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
-        StauX[iatm][nx:] = np.ones_like(Ax)[iatm,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X 32
-        CtauX[iatm][nx:] = np.ones_like(Bx)[iatm,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
+        if nx:    
+            StauX[iatm][:-nx] = np.ones_like(Ax)[iatm,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X 32
+            CtauX[iatm][:-nx] = np.ones_like(Bx)[iatm,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
+        else:
+            StauX[iatm] = np.ones_like(Ax)[iatm,:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X 32
+            CtauX[iatm] = np.ones_like(Bx)[iatm,:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
 #    GSASIIpath.IPyBreak()
     if nWaves[0]:
@@ -1109 +1139 @@
         dGdMuSb = np.sum(part1[:,:,:,:,nxs]*dUdBu*np.sin(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
         dGdMuS = np.concatenate((dGdMuSa,dGdMuSb),axis=-1)   #ops x atoms x waves x 12uij
-#        #GSASIIpath.IPyBreak()
     else:
         HbH = np.ones_like(HdotX)
@@ -1130 +1159 @@
     dGdMzSx = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
     dGdMzS = np.concatenate((dGdMzSt,dGdMzSx),axis=-1)
+#    GSASIIpath.IPyBreak()
 # ops x atoms x waves x 2xyz - imaginary part - good
-#    GSASIIpath.IPyBreak()
     return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS]