Changeset 3888


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

provide contour plot (not finished yet) on structure drawing

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIlattice.py

    r3801 r3888  
    891891    return rdsq
    892892   
    893 def PlaneIntercepts(Amat,Bmat,H,phase,stack):
     893def PlaneIntercepts(Amat,H,phase,stack):
    894894    ''' find unit cell intercepts for a stack of hkl planes
    895895    '''
    … …  
    906906               for j in [0,1,2,3]:
    907907                    hx = [0,0,0]
    908                     intcpt = (phase/360.+step-H[h]*Ux[j,0]-H[k]*Ux[j,1])/H[l]
     908                    intcpt = ((phase)/360.+step-H[h]*Ux[j,0]-H[k]*Ux[j,1])/H[l]
    909909                    if 0. <= intcpt <= 1.:                       
    910910                        hx[h] = Ux[j,0]
    … …  
    924924            A[1:] = [np.dot(DX[1],dx) for dx in DX[1:]]
    925925            HX = HX[np.argsort(A)]
    926 #            GSASIIpath.IPyBreak()
    927926            Stack.append(HX)
    928927    return Stack

  • trunk/GSASIImath.py

    r3877 r3888  
    36503650    return R
    36513651       
     3652def getRhos(XYZ,rho):
     3653    ''' get scattering density at a point by 8-point interpolation
     3654    param xyz:  array coordinates to be probed Nx3
     3655    param: rho: array copy of map (NB: don't use original!)
     3656   
     3657    :returns: density at xyz
     3658    '''
     3659    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
     3660    mapShape = np.array(rho.shape)
     3661    mapStep = 1./mapShape
     3662    X = XYZ%1.    #get into unit cell
     3663    I = np.array(np.rint(X*mapShape),dtype='int')
     3664    R = []
     3665    for i,x in enumerate(X):
     3666        D = x-I[i]*mapStep         #position inside map cell
     3667        D12 = D[0]*D[1]
     3668        D13 = D[0]*D[2]
     3669        D23 = D[1]*D[2]
     3670        D123 = np.prod(D)
     3671        Rho = rollMap(rho,-I[i])    #shifts map so point is in corner
     3672        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]+  \
     3673            Rho[1,1,1]*D123+Rho[0,1,1]*(D23-D123)+Rho[1,0,1]*(D13-D123)+Rho[1,1,0]*(D12-D123)+  \
     3674            Rho[0,0,0]*(D12+D13+D23-D123)-Rho[0,0,1]*(D13+D23-D123)-    \
     3675            Rho[0,1,0]*(D23+D12-D123)-Rho[1,0,0]*(D13+D12-D123))
     3676    return R
     3677       
    36523678def SearchMap(generalData,drawingData,Neg=False):
    36533679    '''Does a search of a density map for peaks meeting the criterion of peak

  • trunk/GSASIIphsGUI.py

    r3850 r3888  
    52185218        if 'showRigidBodies' not in drawingData:
    52195219            drawingData['showRigidBodies'] = True
     5220        if 'showSlice' not in drawingData:
     5221            drawingData['showSlice'] = False
    52205222        if 'Plane' not in drawingData:
    52215223            drawingData['Plane'] = [[0,0,1],False,False,0.0,[255,255,0]]
    … …  
    62496251                G2plt.PlotStructure(G2frame,data)
    62506252               
     6253            def OnShowSlice(event):
     6254                drawingData['showSlice'] = showCS.GetValue()
     6255                G2plt.PlotStructure(G2frame,data)
     6256               
    62516257            def OnViewPoint(event):
    62526258                event.Skip()
    … …  
    63376343           
    63386344            showSizer.Add(line2Sizer)
     6345           
     6346            line3Sizer = wx.BoxSizer(wx.HORIZONTAL)
     6347           
     6348            showCS = wx.CheckBox(drawOptions,-1,label=' Show contour slice?')
     6349            showCS.Bind(wx.EVT_CHECKBOX, OnShowSlice)
     6350            showCS.SetValue(drawingData['showSlice'])           
     6351            line3Sizer.Add(showCS,0,WACV)
     6352           
     6353            showSizer.Add(line3Sizer)
     6354           
    63396355            return showSizer
    63406356           

  • trunk/GSASIIplot.py

    r3868 r3888  
    75577557    wxOrange = wx.Colour(255,128,0)
    75587558    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]])
     7559    eBox = np.array([[.125,.875],[.125,.125],[.9,.125],[.9,.875],])
     7560    eplane = np.array([[-1,-1,0],[-1,1,0],[1,1,0],[1,-1,0]])
    75597561    uEdges = np.array([
    75607562        [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]],
    … …  
    82128214        GL.glPolygonMode(GL.GL_FRONT_AND_BACK,GL.GL_FILL)
    82138215        GL.glFrontFace(GL.GL_CW)
    8214         GL.glBegin(GL.GL_TRIANGLE_FAN)
     8216        GL.glBegin(GL.GL_POLYGON)
    82158217        for vertex in plane:
    82168218            GL.glVertex3fv(vertex)
    82178219        GL.glEnd()
    82188220        GL.glPopMatrix()
     8221        GL.glDisable(GL.GL_BLEND)
     8222        GL.glShadeModel(GL.GL_SMOOTH)
     8223               
     8224    def RenderViewPlane(plane,color,Z,width,height):
     8225        fade = list(color) + [.5,]
     8226        GL.glShadeModel(GL.GL_FLAT)
     8227        ID = GL.glGenTextures(1)
     8228        GL.glBindTexture(GL.GL_TEXTURE_2D, ID)
     8229        GL.glPixelStorei(GL.GL_UNPACK_ALIGNMENT,1)
     8230        GL.glBlendFunc(GL.GL_SRC_ALPHA,GL.GL_ONE_MINUS_SRC_ALPHA)
     8231        GL.glEnable(GL.GL_BLEND)
     8232        GL.glEnable(GL.GL_TEXTURE_2D)
     8233        GL.glPushMatrix()
     8234        GL.glLoadIdentity()
     8235        GL.glShadeModel(GL.GL_SMOOTH)
     8236        GL.glPolygonMode(GL.GL_FRONT_AND_BACK,GL.GL_FILL)
     8237        GL.glFrontFace(GL.GL_CW)
     8238        GL.glTexEnvf(GL.GL_TEXTURE_ENV, GL.GL_TEXTURE_ENV_MODE, GL.GL_REPLACE)
     8239        GL.glBindTexture(GL.GL_TEXTURE_2D, ID)
     8240        GL.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_BASE_LEVEL, 0)
     8241        GL.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_MAX_LEVEL, 0)
     8242        GL.glTexImage2D(GL.GL_TEXTURE_2D,0,GL.GL_RGBA,width,height,0,GL.GL_RGBA,GL.GL_UNSIGNED_BYTE,Z)
     8243        GL.glBegin(GL.GL_POLYGON)
     8244        for vertex,evertex in zip(plane,eBox):
     8245            GL.glTexCoord2fv(evertex)
     8246            GL.glVertex3fv(vertex)
     8247        GL.glEnd()
     8248#        GL.glDrawPixels(width,height,GL.GL_RGBA,GL.GL_UNSIGNED_BYTE,Z)
     8249        GL.glPopMatrix()
     8250        GL.glDisable(GL.GL_TEXTURE_2D)
    82198251        GL.glDisable(GL.GL_BLEND)
    82208252        GL.glShadeModel(GL.GL_SMOOTH)
    … …  
    84678499        GL.glInitNames()
    84688500        GL.glPushName(0)
    8469        
     8501           
    84708502        GL.glMatrixMode(GL.GL_PROJECTION)
    84718503        GL.glLoadIdentity()
    … …  
    84738505        GLU.gluPerspective(20.,aspect,cPos-Zclip,cPos+Zclip)
    84748506        GLU.gluLookAt(0,0,cPos,0,0,0,0,1,0)
    8475         SetLights()           
     8507        SetLights()
    84768508           
    84778509        GL.glMatrixMode(GL.GL_MODELVIEW)
    … …  
    86388670            if drawingData['Plane'][1]:
    86398671                H,phase,stack,phase,color = drawingData['Plane']
    8640                 Planes = G2lat.PlaneIntercepts(Amat,Bmat,H,phase,stack)
     8672                Planes = G2lat.PlaneIntercepts(Amat,H,phase,stack)
    86418673                for plane in Planes:
    86428674                    RenderPlane(plane,color)
     8675            if drawingData['showSlice']:
     8676                rho = generalData['Map']['rho']
     8677                if not len(rho):
     8678                    return
     8679                from matplotlib.backends.backend_agg import FigureCanvasAgg
     8680                import matplotlib.pyplot as plt
     8681                Model = GL.glGetDoublev(GL.GL_MODELVIEW_MATRIX)
     8682                invModel = nl.inv(Model)
     8683                msize = 5.      #-5A - 5A slice
     8684                npts = int(2*msize/0.25)
     8685                VP = np.array(drawingData['viewPoint'][0])
     8686                SX,SY = np.meshgrid(np.linspace(-1.,1.,npts),np.linspace(-1.,1.,npts))
     8687                SXYZ = msize*np.dstack((SX,SY,np.zeros_like(SX)))
     8688                SXYZ = np.reshape(np.inner(SXYZ,invModel[:3,:3].T)+VP[nxs,nxs,:],(-1,3))
     8689                Z = np.reshape(G2mth.getRhos(SXYZ,rho),(npts,npts))
     8690                plt.contour(Z,colors='k',linewidths=1)
     8691                plt.axis("off")
     8692                canvas = plt.get_current_fig_manager().canvas
     8693                agg = canvas.switch_backends(FigureCanvasAgg)
     8694                agg.draw()
     8695                img, (width, height) = agg.print_to_buffer()
     8696                Zimg = np.frombuffer(img, np.uint8).reshape((height, width, 4))
     8697                RenderViewPlane(msize*eplane,Wt,Zimg,width,height)
     8698               
    86438699#        print time.time()-time0
    86448700        try:
Note: See TracChangeset for help on using the changeset viewer.