Changeset 2476


Ignore:
Timestamp:
Sep 22, 2016 11:42:28 AM (5 years ago)
Author:
vondreele
Message:

Add new feature to structure plots - display an hkl plane or stack of planes with a phase offset.
fix lighting issues for planes/polygons

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIlattice.py

    r2473 r2476  
    657657    rdsq = calc_rDsq(H[:3]+(H[3][:,np.newaxis]*vec).T,A)+Z/difC
    658658    return rdsq
     659   
     660def PlaneIntercepts(Amat,Bmat,H,phase,stack):
     661    ''' find unit cell intercepts for a stack of hkl planes
     662    '''
     663    Steps = range(-1,2,2)
     664    if stack:
     665        Steps = range(-10,10,1)
     666    Stack = []
     667    Ux = np.array([[0,0],[1,0],[1,1],[0,1]])
     668    for step in Steps:
     669        HX = []
     670        for i in [0,1,2]:
     671            if H[i]:
     672               h,k,l = [(i+1)%3,(i+2)%3,(i+3)%3]
     673               for j in [0,1,2,3]:
     674                    hx = [0,0,0]
     675                    intcpt = (phase/360.+step-H[h]*Ux[j,0]-H[k]*Ux[j,1])/H[l]
     676                    if 0. <= intcpt <= 1.:                       
     677                        hx[h] = Ux[j,0]
     678                        hx[k] = Ux[j,1]
     679                        hx[l] = intcpt
     680                        HX.append(hx)
     681        if len(HX)> 2:
     682            HX = np.array(HX)
     683            DX = np.inner(HX-HX[0],Amat)
     684            D = np.sqrt(np.sum(DX**2,axis=1))
     685            Dsort = np.argsort(D)
     686            HX = HX[Dsort]
     687            DX = DX[Dsort]
     688            D = D[Dsort]
     689            DX[1:,:] = DX[1:,:]/D[1:,nxs]
     690            A = 2.*np.ones(HX.shape[0])
     691            A[1:] = [np.dot(DX[1],dx) for dx in DX[1:]]
     692            HX = HX[np.argsort(A)]
     693#            GSASIIpath.IPyBreak()
     694            Stack.append(HX)
     695    return Stack
    659696       
    660697def MaxIndex(dmin,A):

  • trunk/GSASIIphsGUI.py

    r2475 r2476  
    38123812            'vdwScale':0.67,'ellipseProb':50,'sizeH':0.50,'unitCellBox':True,
    38133813            'showABC':True,'selectedAtoms':[],'Atoms':[],'oldxy':[],
    3814             'bondList':{},'viewDir':[1,0,0]}
     3814            'bondList':{},'viewDir':[1,0,0],'Plane':[[0,0,1],False,False,0.0,[255,255,0]]}
    38153815        V0 = np.array([0,0,1])
    38163816        V = np.inner(Amat,V0)
    … …  
    38353835        if 'showRigidBodies' not in drawingData:
    38363836            drawingData['showRigidBodies'] = True
     3837        if 'Plane' not in drawingData:
     3838            drawingData['Plane'] = [[0,0,1],False,False,0.0,[255,255,0]]
    38373839        cx,ct,cs,ci = [0,0,0,0]
    38383840        if generalData['Type'] in ['nuclear','faulted',]:
    … …  
    38473849                DrawAtomAdd(drawingData,atom)
    38483850            data['Drawing'] = drawingData
     3851        if len(drawingData['Plane']) < 5:
     3852            drawingData['Plane'].append([255,255,0])
    38493853           
    38503854    def DrawAtomAdd(drawingData,atom):
    … …  
    49824986            radSizer.Add(radFactor,0,WACV)
    49834987            return radSizer
     4988           
     4989        def PlaneSizer():
     4990           
     4991            def OnPlane(event):
     4992                event.Skip()
     4993                vals = plane.GetValue().split()
     4994                try:
     4995                    hkl = [int(vals[i]) for i in range(3)]
     4996                    if not np.any(np.array(hkl)):       #can't be all zeros!
     4997                        raise ValueError
     4998                except (ValueError,IndexError):
     4999                    hkl = drawingData['Plane'][0]
     5000                drawingData['Plane'][0] = hkl
     5001                plane.SetValue('%3d %3d %3d'%(hkl[0],hkl[1],hkl[2]))
     5002                G2plt.PlotStructure(G2frame,data)
     5003               
     5004            def OnShowPlane(event):
     5005                drawingData['Plane'][1] = showPlane.GetValue()
     5006                G2plt.PlotStructure(G2frame,data)
     5007               
     5008            def OnShowStack(event):
     5009                drawingData['Plane'][2] = showStack.GetValue()
     5010                G2plt.PlotStructure(G2frame,data)
     5011               
     5012            def OnPhase(event):
     5013                event.Skip()
     5014                try:
     5015                    val = float(phase.GetValue())
     5016                except ValueError:
     5017                    val = drawingData['Plane'][3]
     5018                drawingData['Plane'][3] = val
     5019                phase.SetValue('%.2f'%(val))
     5020                G2plt.PlotStructure(G2frame,data)
     5021           
     5022            def OnPlaneColor(event):
     5023                drawingData['Plane'][4] = event.GetValue()
     5024                G2plt.PlotStructure(G2frame,data)
     5025
     5026            planeSizer = wx.BoxSizer(wx.VERTICAL)
     5027            planeSizer1 = wx.BoxSizer(wx.HORIZONTAL)
     5028            planeSizer1.Add(wx.StaticText(drawOptions,label=' Plane: '),0,WACV)
     5029            H = drawingData['Plane'][0]
     5030            plane = wx.TextCtrl(drawOptions,value='%3d %3d %3d'%(H[0],H[1],H[2]),
     5031                style=wx.TE_PROCESS_ENTER)
     5032            plane.Bind(wx.EVT_TEXT_ENTER,OnPlane)
     5033            plane.Bind(wx.EVT_KILL_FOCUS,OnPlane)
     5034            planeSizer1.Add(plane,0,WACV)
     5035            showPlane = wx.CheckBox(drawOptions,label=' Show plane?')
     5036            showPlane.SetValue(drawingData['Plane'][1])
     5037            showPlane.Bind(wx.EVT_CHECKBOX, OnShowPlane)
     5038            planeSizer1.Add(showPlane,0,WACV)
     5039            showStack = wx.CheckBox(drawOptions,label=' As a stack?')
     5040            showStack.SetValue(drawingData['Plane'][2])
     5041            showStack.Bind(wx.EVT_CHECKBOX, OnShowStack)
     5042            planeSizer1.Add(showStack,0,WACV)
     5043            planeSizer2 = wx.BoxSizer(wx.HORIZONTAL)
     5044            planeSizer2.Add(wx.StaticText(drawOptions,label=' Phase shift (deg): '),0,WACV)
     5045            phase = wx.TextCtrl(drawOptions,value='%.2f'%(drawingData['Plane'][3]),
     5046                style=wx.TE_PROCESS_ENTER)
     5047            phase.Bind(wx.EVT_TEXT_ENTER,OnPhase)
     5048            phase.Bind(wx.EVT_KILL_FOCUS,OnPhase)
     5049            planeSizer2.Add(phase,0,WACV)
     5050            planeSizer2.Add(wx.StaticText(drawOptions,-1,' Plane color: '),0,WACV)
     5051            planeColor = wcs.ColourSelect(drawOptions, -1,colour=drawingData['Plane'][4],size=wx.Size(25,25))
     5052            planeColor.Bind(wcs.EVT_COLOURSELECT, OnPlaneColor)
     5053            planeSizer2.Add(planeColor,0,WACV)
     5054            planeSizer.Add(planeSizer1)
     5055            planeSizer.Add(planeSizer2)
     5056            return planeSizer
     5057           
    49845058
    49855059        # UpdateDrawOptions exectable code starts here
    … …  
    50055079        mainSizer.Add((5,5),0)
    50065080        mainSizer.Add(RadSizer(),0,)
     5081        mainSizer.Add((5,5),0)
     5082        mainSizer.Add(PlaneSizer(),0,)
    50075083        SetPhaseWindow(G2frame.dataFrame,drawOptions,mainSizer)
    50085084

  • trunk/GSASIIplot.py

    r2475 r2476  
    47494749def PlotStructure(G2frame,data,firstCall=False):
    47504750    '''Crystal structure plotting package. Can show structures as balls, sticks, lines,
    4751     thermal motion ellipsoids and polyhedra
     4751    thermal motion ellipsoids and polyhedra. Magnetic moments shown as black/red
     4752    arrows according to spin state
    47524753    '''
    47534754
    … …  
    52115212        glEnable(GL_LIGHT0)
    52125213        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,0)
    5213         glLightfv(GL_LIGHT0,GL_AMBIENT,[.5,.5,.5,1])
     5214        glLightfv(GL_LIGHT0,GL_AMBIENT,[.8,.8,.8,1])
    52145215        glLightfv(GL_LIGHT0,GL_DIFFUSE,[.8,.8,.8,1])
    5215         glLightfv(GL_LIGHT0,GL_SPECULAR,[1,1,1,1])
     5216#        glLightfv(GL_LIGHT0,GL_SPECULAR,[1,1,1,1])
     5217#        glLightfv(GL_LIGHT0,GL_POSITION,[0,0,1,1])
    52165218       
    52175219    def GetRoll(newxy,rhoshape):
    … …  
    54045406        glColor4ubv([0,0,0,0])
    54055407        glDisable(GL_COLOR_MATERIAL)
     5408       
     5409    def RenderPlane(plane,color):
     5410        fade = list(color) + [.25,]
     5411        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,fade)
     5412        glShadeModel(GL_FLAT)
     5413        glEnable(GL_BLEND)
     5414        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
     5415        glPushMatrix()
     5416        glShadeModel(GL_SMOOTH)
     5417        glPolygonMode(GL_FRONT_AND_BACK,GL_FILL)
     5418        glFrontFace(GL_CW)
     5419        glBegin(GL_TRIANGLE_FAN)
     5420        for vertex in plane:
     5421            glVertex3fv(vertex)
     5422        glEnd()
     5423        glPopMatrix()
     5424        glDisable(GL_BLEND)
     5425        glShadeModel(GL_SMOOTH)
    54065426               
    54075427    def RenderSphere(x,y,z,radius,color):
    … …  
    55185538        glPushMatrix()
    55195539        glTranslate(x,y,z)
    5520         glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
     5540        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,color)
    55215541        glShadeModel(GL_SMOOTH)
    55225542        glMultMatrixf(B4mat.T)
    … …  
    55995619#useful debug?       
    56005620#        if caller:
    5601 #            print caller
     5621#            print caller,generalData['Name']
    56025622# end of useful debug
    56035623        mapData = generalData['Map']
    … …  
    56665686        glMultMatrixf(A4mat.T)
    56675687        glTranslate(-Tx,-Ty,-Tz)
    5668         if drawingData['unitCellBox']:
    5669             RenderBox()
    56705688        if drawingData['showABC']:
    56715689            x,y,z = drawingData['viewPoint'][0]
    … …  
    58205838                Backbone = Backbones[chain]
    58215839                RenderBackbone(Backbone,BackboneColor,bondR)
     5840        if drawingData['unitCellBox']:
     5841            RenderBox()
     5842            if drawingData['Plane'][1]:
     5843                H,phase,stack,phase,color = drawingData['Plane']
     5844                Planes = G2lat.PlaneIntercepts(Amat,Bmat,H,phase,stack)
     5845                for plane in Planes:
     5846                    RenderPlane(plane,color)
    58225847#        print time.time()-time0
    58235848        try:
Note: See TracChangeset for help on using the changeset viewer.