# Changeset 432

Ignore:
Timestamp:
Dec 5, 2011 4:00:00 PM (12 years ago)
Message:

remove test prints from GSASIIIO
new Gmat2AB routine
work on the DData GUI bug - doesn't crash but bad GUI shown
work on ellipsoidal crystallites - all now work.

Location:
trunk
Files:
6 edited

Unmodified
Removed
• ## trunk/GSASIIIO.py

 r430 File = open(filename,'Ur') cons = Bank.split() print cons Nch = int(cons[2]) if DataType[2] == 'C': y.append(yi) w.append(1.0/vi) print '%8.3f %8.0f %d %8g'%(xi,yi,ni,1./vi) S = File.readline() File.close() N = len(x) print N return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
• ## trunk/GSASIIlattice.py

 r395 G,g = A2Gmat(A) return Gmat2cell(G) def Gmat2AB(G): """Computes orthogonalization matrix from reciprocal metric tensor G returns tuple of two 3x3 numpy arrays (A,B) A for crystal to Cartesian transformations A*x = np.inner(A,x) = X B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B,X) = x """ cellstar = Gmat2cell(G) g = nl.inv(G) cell = Gmat2cell(g) A = np.zeros(shape=(3,3)) # from Giacovazzo (Fundamentals 2nd Ed.) p.75 A[0][0] = cell[0]                # a A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma) A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta) A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma) A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta) A[2][2] = 1/cellstar[2]         # 1/c* B = nl.inv(A) return A,B def cell2AB(cell): returns tuple of two 3x3 numpy arrays (A,B) A for crystal to Cartesian transformations A*x = np.inner(A,x) = X B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B*x) = x B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B,X) = x """ G,g = cell2Gmat(cell)

• ## trunk/GSASIIplot.py

 r431 axes /= nl.norm(axes) Shkl = np.array(coeff[1]) Shape = X.shape[0] XYZ = np.dstack((X,Y,Z)) XYZ = np.nan_to_num(np.apply_along_axis(uniaxCalc,2,XYZ,iso,aniso,axes)) Z = Z[:,:,0] elif coeff[0] == 'ellipsoidal': def ellipseCalc(xyz,E,R): XYZ = xyz*E.T return np.inner(XYZ.T,R) S6 = coeff[4] Sij = G2lat.U6toUij(S6) E,R = nl.eigh(Sij) XYZ = np.dstack((X,Y,Z)) XYZ = np.nan_to_num(np.apply_along_axis(ellipseCalc,2,XYZ,E,R)) X,Y,Z = np.dsplit(XYZ,3) X = X[:,:,0] Y = Y[:,:,0] Z = Z[:,:,0] elif coeff[0] == 'generalized': Shkl = np.array(coeff[4]) if np.any(Shkl): Shape = X.shape[0] XYZ = np.dstack((X,Y,Z)) XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,SGData,A,Shkl)) Z = Z[:,:,0] elif coeff[0] == 'ellipsoidal': print 'ellipsoid plot' if np.any(X) and np.any(Y) and np.any(Z): errFlags = np.seterr(all='ignore') Plot.plot_surface(X,Y,Z,rstride=1,cstride=1,color='g',linewidth=1) np.seterr(all='ignore') xyzlim = np.array([Plot.get_xlim3d(),Plot.get_ylim3d(),Plot.get_zlim3d()]).T XYZlim = [min(xyzlim[0]),max(xyzlim[1])] Plot.set_xlim3d(XYZlim) Plot.set_ylim3d(XYZlim) Plot.set_zlim3d(XYZlim) Plot.set_aspect('equal') if plotType == 'Size': Plot.set_title('Crystallite size for '+phase+'\n'+coeff[0]+' model') Plot.set_ylabel(r'Y, $\mu$m') Plot.set_zlabel(r'Z, $\mu$m') Plot.set_aspect('equal') else: Plot.set_title(r'$\mu$strain for '+phase+'\n'+coeff[0]+' model') Plot.set_ylabel(r'Y, $\mu$strain') Plot.set_zlabel(r'Z, $\mu$strain') Plot.set_aspect('equal') else: h,k,l = generalData['POhkl'] def PlotTexture(self,data,newPlot=False,Start=False): '''Pole figure, inverse pole figure(?), 3D pole distribution and 3D inverse pole distribution(?) plotting; Need way to select pole figure or pole distribution to be displayed - do in key enter menu '''Pole figure, inverse pole figure, 3D pole distribution and 3D inverse pole distribution plotting. dict generalData contains all phase info needed which is in data ''' pf = G2lat.polfcal(ODFln,SamSym[textureData['Model']],np.array([r,]),np.array([p,])) self.G2plotNB.status.SetFields(['','phi =%9.3f, gam =%9.3f, MRD =%9.3f'%(r,p,pf)]) Plot = self.G2plotNB.addMpl('Texture').gca() if self.Projection == '3D display' and 'Axial'  not in SHData['PlotType']: Plot = mp3d.Axes3D(self.G2plotNB.add3D('Texture')) else: Plot = self.G2plotNB.addMpl('Texture').gca() plotNum = self.G2plotNB.plotList.index('Texture') Page = self.G2plotNB.nb.GetPage(plotNum) Page.canvas.mpl_connect('motion_notify_event', OnMotion) if self.Projection == '3D display': pass else: Page.canvas.mpl_connect('motion_notify_event', OnMotion) Page.SetFocus()
• ## trunk/GSASIIpwd.py

 r429 sumDf = np.sum(Df) return Df,dFdp,dFds,dFdg,dFdsh def ellipseSize(H,Sij,GB): HX = np.inner(H.T,GB) lenHX = np.sqrt(np.sum(HX**2)) Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij)) R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal lenR = np.sqrt(np.sum(R**2)) return lenR def ellipseSizeDerv(H,Sij,GB): lenR = ellipseSize(H,Sij,GB) delt = 0.001 dRdS = np.zeros(6) for i in range(6): dSij = Sij[:] dSij[i] += delt dRdS[i] = (ellipseSize(H,dSij,GB)-lenR)/delt return lenR,dRdS def getPeakProfile(parmDict,xdata,varyList,bakType):
• ## trunk/GSASIIstruct.py

 r429 return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse): def GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict): costh = cosd(refl[5]/2.) #crystallite size gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1']*costh) gam *= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2) else:           #ellipsoidal crystallites - wrong not make sense else:           #ellipsoidal crystallites Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] H = np.array(refl[:3]) gam += 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H))) lenR = G2pwd.ellipseSize(H,Sij,GB) gam = 1.8*wave/(np.pi*costh*lenR) #microstrain if calcControls[phfx+'MustrainType'] == 'isotropic': return gam def GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse): def GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict): gamDict = {} costh = cosd(refl[5]/2.) Si = parmDict[phfx+'Size:0'] Sa = parmDict[phfx+'Size:1'] gami = (1.80*wave/np.pi)/(Si*Sa) gami = (1.8*wave/np.pi)/(Si*Sa) sqtrm = np.sqrt((cosP*Sa)**2+(sinP*Si)**2) gam = gami*sqtrm/costh gamDict[phfx+'Size:0'] = gami*Si*sinP**2/(sqtrm*costh)-gam/Si gamDict[phfx+'Size:1'] = gami*Sa*cosP**2/(sqtrm*costh)-gam/Sa else:           #ellipsoidal crystallites - do numerically? - not right not make sense else:           #ellipsoidal crystallites const = 1.8*wave/(np.pi*costh) Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] H = np.array(refl[:3]) gam = 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H))) R,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB) for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]): gamDict[item] = -(const/R**2)*dRdS[i] #microstrain derivatives if calcControls[phfx+'MustrainType'] == 'isotropic': def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): def GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse): def GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict): U = parmDict[hfx+'U'] V = parmDict[hfx+'V'] sig = U*tanPos**2+V*tanPos+W        #save peak sigma sig = max(0.001,sig) gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict) #save peak gamma gam = max(0.001,gam) return sig,gam A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] G,g = G2lat.A2Gmat(A)       #recip & real metric tensors GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies Vst = np.sqrt(nl.det(G))    #V* if 'Pawley' not in Phase['General']['Type']: refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict) sizeEllipse = [] if calcControls[phfx+'SizeType'] == 'ellipsoidal': sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)]) for refl in refList: if 'C' in calcControls[hfx+'histType']: Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift refl[6:8] = GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse)    #peak sig & gam refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13] refl[13] *= Vst*Lorenz A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] G,g = G2lat.A2Gmat(A)       #recip & real metric tensors GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies if 'Pawley' not in Phase['General']['Type']: dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict) sizeEllipse = [] if calcControls[phfx+'SizeType'] == 'ellipsoidal': sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)]) for iref,refl in enumerate(refList): if 'C' in calcControls[hfx+'histType']:        #CW powder elif name in dependentVars: depDerivDict[name] += dDijDict[name]*dervDict['pos'] gamDict = GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict) for name in gamDict: if name in varylist: SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData) #for testing purposes!!! #    file = open('structTestdata.dat','wb') #    cPickle.dump(parmDict,file,1) #    cPickle.dump(varyList,file,1) #    for histogram in Histograms: #        if 'PWDR' in histogram[:4]: #            Histogram = Histograms[histogram] #    cPickle.dump(Histogram,file,1) #    cPickle.dump(Phases,file,1) #    cPickle.dump(calcControls,file,1) #    cPickle.dump(pawleyLookup,file,1) #    file.close() file = open('structTestdata.dat','wb') cPickle.dump(parmDict,file,1) cPickle.dump(varyList,file,1) for histogram in Histograms: if 'PWDR' in histogram[:4]: Histogram = Histograms[histogram] cPickle.dump(Histogram,file,1) cPickle.dump(Phases,file,1) cPickle.dump(calcControls,file,1) cPickle.dump(pawleyLookup,file,1) file.close() def SeqRefine(GPXfile,dlg):
Note: See TracChangeset for help on using the changeset viewer.