# Changeset 3888

Ignore:
Timestamp:
Apr 11, 2019 3:59:48 PM (5 years ago)
Message:

provide contour plot (not finished yet) on structure drawing

Location:
trunk
Files:
4 edited

Unmodified
Removed
• ## trunk/GSASIIlattice.py

 r3801 return rdsq def PlaneIntercepts(Amat,Bmat,H,phase,stack): def PlaneIntercepts(Amat,H,phase,stack): ''' find unit cell intercepts for a stack of hkl planes ''' for j in [0,1,2,3]: hx = [0,0,0] intcpt = (phase/360.+step-H[h]*Ux[j,0]-H[k]*Ux[j,1])/H[l] intcpt = ((phase)/360.+step-H[h]*Ux[j,0]-H[k]*Ux[j,1])/H[l] if 0. <= intcpt <= 1.: hx[h] = Ux[j,0] A[1:] = [np.dot(DX[1],dx) for dx in DX[1:]] HX = HX[np.argsort(A)] #            GSASIIpath.IPyBreak() Stack.append(HX) return Stack
• ## trunk/GSASIImath.py

 r3877 return R def getRhos(XYZ,rho): ''' get scattering density at a point by 8-point interpolation param xyz:  array coordinates to be probed Nx3 param: rho: array copy of map (NB: don't use original!) :returns: density at xyz ''' rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2) mapShape = np.array(rho.shape) mapStep = 1./mapShape X = XYZ%1.    #get into unit cell I = np.array(np.rint(X*mapShape),dtype='int') R = [] for i,x in enumerate(X): D = x-I[i]*mapStep         #position inside map cell D12 = D[0]*D[1] D13 = D[0]*D[2] D23 = D[1]*D[2] D123 = np.prod(D) Rho = rollMap(rho,-I[i])    #shifts map so point is in corner R.append(Rho[0,0,0]*(1.-np.sum(D))+Rho[1,0,0]*D[0]+Rho[0,1,0]*D[1]+Rho[0,0,1]*D[2]+  \ Rho[1,1,1]*D123+Rho[0,1,1]*(D23-D123)+Rho[1,0,1]*(D13-D123)+Rho[1,1,0]*(D12-D123)+  \ Rho[0,0,0]*(D12+D13+D23-D123)-Rho[0,0,1]*(D13+D23-D123)-    \ Rho[0,1,0]*(D23+D12-D123)-Rho[1,0,0]*(D13+D12-D123)) return R def SearchMap(generalData,drawingData,Neg=False): '''Does a search of a density map for peaks meeting the criterion of peak
• ## trunk/GSASIIphsGUI.py

 r3850 if 'showRigidBodies' not in drawingData: drawingData['showRigidBodies'] = True if 'showSlice' not in drawingData: drawingData['showSlice'] = False if 'Plane' not in drawingData: drawingData['Plane'] = [[0,0,1],False,False,0.0,[255,255,0]] G2plt.PlotStructure(G2frame,data) def OnShowSlice(event): drawingData['showSlice'] = showCS.GetValue() G2plt.PlotStructure(G2frame,data) def OnViewPoint(event): event.Skip() showSizer.Add(line2Sizer) line3Sizer = wx.BoxSizer(wx.HORIZONTAL) showCS = wx.CheckBox(drawOptions,-1,label=' Show contour slice?') showCS.Bind(wx.EVT_CHECKBOX, OnShowSlice) showCS.SetValue(drawingData['showSlice']) line3Sizer.Add(showCS,0,WACV) showSizer.Add(line3Sizer) return showSizer
• ## trunk/GSASIIplot.py

 r3868 wxOrange = wx.Colour(255,128,0) uBox = np.array([[0,0,0],[1,0,0],[1,1,0],[0,1,0],[0,0,1],[1,0,1],[1,1,1],[0,1,1]]) eBox = np.array([[.125,.875],[.125,.125],[.9,.125],[.9,.875],]) eplane = np.array([[-1,-1,0],[-1,1,0],[1,1,0],[1,-1,0]]) uEdges = np.array([ [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]], GL.glPolygonMode(GL.GL_FRONT_AND_BACK,GL.GL_FILL) GL.glFrontFace(GL.GL_CW) GL.glBegin(GL.GL_TRIANGLE_FAN) GL.glBegin(GL.GL_POLYGON) for vertex in plane: GL.glVertex3fv(vertex) GL.glEnd() GL.glPopMatrix() GL.glDisable(GL.GL_BLEND) GL.glShadeModel(GL.GL_SMOOTH) def RenderViewPlane(plane,color,Z,width,height): fade = list(color) + [.5,] GL.glShadeModel(GL.GL_FLAT) ID = GL.glGenTextures(1) GL.glBindTexture(GL.GL_TEXTURE_2D, ID) GL.glPixelStorei(GL.GL_UNPACK_ALIGNMENT,1) GL.glBlendFunc(GL.GL_SRC_ALPHA,GL.GL_ONE_MINUS_SRC_ALPHA) GL.glEnable(GL.GL_BLEND) GL.glEnable(GL.GL_TEXTURE_2D) GL.glPushMatrix() GL.glLoadIdentity() GL.glShadeModel(GL.GL_SMOOTH) GL.glPolygonMode(GL.GL_FRONT_AND_BACK,GL.GL_FILL) GL.glFrontFace(GL.GL_CW) GL.glTexEnvf(GL.GL_TEXTURE_ENV, GL.GL_TEXTURE_ENV_MODE, GL.GL_REPLACE) GL.glBindTexture(GL.GL_TEXTURE_2D, ID) GL.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_BASE_LEVEL, 0) GL.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_MAX_LEVEL, 0) GL.glTexImage2D(GL.GL_TEXTURE_2D,0,GL.GL_RGBA,width,height,0,GL.GL_RGBA,GL.GL_UNSIGNED_BYTE,Z) GL.glBegin(GL.GL_POLYGON) for vertex,evertex in zip(plane,eBox): GL.glTexCoord2fv(evertex) GL.glVertex3fv(vertex) GL.glEnd() #        GL.glDrawPixels(width,height,GL.GL_RGBA,GL.GL_UNSIGNED_BYTE,Z) GL.glPopMatrix() GL.glDisable(GL.GL_TEXTURE_2D) GL.glDisable(GL.GL_BLEND) GL.glShadeModel(GL.GL_SMOOTH) GL.glInitNames() GL.glPushName(0) GL.glMatrixMode(GL.GL_PROJECTION) GL.glLoadIdentity() GLU.gluPerspective(20.,aspect,cPos-Zclip,cPos+Zclip) GLU.gluLookAt(0,0,cPos,0,0,0,0,1,0) SetLights() SetLights() GL.glMatrixMode(GL.GL_MODELVIEW) if drawingData['Plane'][1]: H,phase,stack,phase,color = drawingData['Plane'] Planes = G2lat.PlaneIntercepts(Amat,Bmat,H,phase,stack) Planes = G2lat.PlaneIntercepts(Amat,H,phase,stack) for plane in Planes: RenderPlane(plane,color) if drawingData['showSlice']: rho = generalData['Map']['rho'] if not len(rho): return from matplotlib.backends.backend_agg import FigureCanvasAgg import matplotlib.pyplot as plt Model = GL.glGetDoublev(GL.GL_MODELVIEW_MATRIX) invModel = nl.inv(Model) msize = 5.      #-5A - 5A slice npts = int(2*msize/0.25) VP = np.array(drawingData['viewPoint'][0]) SX,SY = np.meshgrid(np.linspace(-1.,1.,npts),np.linspace(-1.,1.,npts)) SXYZ = msize*np.dstack((SX,SY,np.zeros_like(SX))) SXYZ = np.reshape(np.inner(SXYZ,invModel[:3,:3].T)+VP[nxs,nxs,:],(-1,3)) Z = np.reshape(G2mth.getRhos(SXYZ,rho),(npts,npts)) plt.contour(Z,colors='k',linewidths=1) plt.axis("off") canvas = plt.get_current_fig_manager().canvas agg = canvas.switch_backends(FigureCanvasAgg) agg.draw() img, (width, height) = agg.print_to_buffer() Zimg = np.frombuffer(img, np.uint8).reshape((height, width, 4)) RenderViewPlane(msize*eplane,Wt,Zimg,width,height) #        print time.time()-time0 try:
Note: See TracChangeset for help on using the changeset viewer.