Changeset 753 for trunk/GSASIIplot.py


Ignore:
Timestamp:
Sep 7, 2012 1:32:36 PM (10 years ago)
Author:
vondreele
Message:

change structure rotation mechanisms - now use quaternions

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIplot.py

    r750 r753  
    3131import GSASIIlattice as G2lat
    3232import GSASIIspc as G2spc
     33import GSASIImath as G2mth
    3334import pytexture as ptx
    3435from  OpenGL.GL import *
     
    23152316    Vol = generalData['Cell'][7:8][0]
    23162317    Amat,Bmat = G2lat.cell2AB(cell)         #Amat - crystal to cartesian, Bmat - inverse
     2318    Gmat,gmat = G2lat.cell2Gmat(cell)
    23172319    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
    23182320    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
     
    23532355    shiftDown = False
    23542356    ctrlDown = False
    2355     global sumroll
    2356     sumroll = np.zeros(3)
    23572357   
    23582358    def OnKeyBox(event):
     
    23912391        if key in ['C']:
    23922392            drawingData['viewPoint'] = [[.5,.5,.5],[0,0]]
     2393            VD = np.inner(Bmat,np.array([0,0,1]))
     2394            VD /= np.sqrt(np.sum(VD**2))
     2395            drawingData['viewDir'] = VD
    23932396            drawingData['Rotation'] = [0.0,0.0,0.0,[]]
     2397            drawingData['Quaternion'] = [0.0,0.0,1.0,0.0]
    23942398            SetViewPointText(drawingData['viewPoint'][0])
     2399            SetViewDirText(drawingData['viewDir'])
     2400            Q = drawingData['Quaternion']
     2401            G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
    23952402        elif key in ['N']:
    23962403            drawAtoms = drawingData['Atoms']
     
    24292436            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
    24302437        elif key in ['U','D','L','R'] and mapData['Flip'] == True:
    2431             dirDict = {'U':[0,-1],'D':[0,1],'L':[-1,0],'R':[1,0]}
     2438            dirDict = {'U':[0,1],'D':[0,-1],'L':[-1,0],'R':[1,0]}
    24322439            SetMapRoll(dirDict[key])
    24332440            SetPeakRoll(dirDict[key])
     
    24912498            if event.LeftIsDown():
    24922499                SetRotation(newxy)
    2493                 angX,angY,angZ = drawingData['Rotation'][:3]
    2494                 G2frame.G2plotNB.status.SetStatusText('New rotation: %.2f, %.2f ,%.2f'%(angX,angY,angZ),1)
     2500                Q = drawingData['Quaternion']
     2501                G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
    24952502            elif event.RightIsDown():
    24962503                SetTranslation(newxy)
     
    24992506            elif event.MiddleIsDown():
    25002507                SetRotationZ(newxy)
    2501                 angX,angY,angZ = drawingData['Rotation'][:3]
    2502                 G2frame.G2plotNB.status.SetStatusText('New rotation: %.2f, %.2f, %.2f'%(angX,angY,angZ),1)
     2508                Q = drawingData['Quaternion']
     2509                G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
    25032510            Draw()
    25042511       
     
    25322539                names = [child.GetName() for child in panel]
    25332540                panel[names.index('viewPoint')].SetValue('%.3f, %.3f, %.3f'%(VP[0],VP[1],VP[2]))
     2541               
     2542    def SetViewDirText(VD):
     2543        page = getSelection()
     2544        if page:
     2545            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
     2546                panel = G2frame.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
     2547                names = [child.GetName() for child in panel]
     2548                panel[names.index('viewDir')].SetValue('%.3f, %.3f, %.3f'%(VD[0],VD[1],VD[2]))
    25342549               
    25352550    def SetMapPeaksText(mapPeaks):
     
    25932608       
    25942609    def GetRoll(newxy,rho):
    2595         anglex,angley,anglez,oldxy = drawingData['Rotation']
    2596         Rx = G2lat.rotdMat(anglex,0)
    2597         Ry = G2lat.rotdMat(angley,1)
    2598         Rz = G2lat.rotdMat(anglez,2)
    2599         dxy = np.inner(Bmat,np.inner(Rz,np.inner(Ry,np.inner(Rx,newxy+[0,]))))
    2600         dxy *= np.array([-1,-1,1])
    2601         dxy = np.array(dxy*rho.shape)
     2610        Q = drawingData['Quaternion']
     2611        dxy = G2mth.prodQVQ(G2mth.invQ(Q),np.inner(Bmat,newxy+[0,]))
     2612        dxy = np.array(dxy*rho.shape)       
    26022613        roll = np.where(dxy>0.5,1,np.where(dxy<-.5,-1,0))
    26032614        return roll
    26042615               
    26052616    def SetMapRoll(newxy):
    2606         global sumroll
    26072617        rho = mapData['rho']
    26082618        roll = GetRoll(newxy,rho)
    2609         sumroll += roll
    26102619        mapData['rho'] = np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
    2611 #        print 'sumroll',sumroll,rho.shape      #useful debug?
    26122620        drawingData['Rotation'][3] = list(newxy)
    26132621       
     
    26232631               
    26242632    def SetTranslation(newxy):
     2633#first get translation vector in screen coords.       
     2634        oldxy = drawingData['Rotation'][3]
     2635        if not len(oldxy): oldxy = list(newxy)
     2636        dxy = newxy-oldxy
     2637        drawingData['Rotation'][3] = list(newxy)
     2638        V = np.array([-dxy[0],dxy[1],0.])
     2639#then transform to rotated crystal coordinates & apply to view point       
     2640        Q = drawingData['Quaternion']
     2641        V = np.inner(Bmat,G2mth.prodQVQ(G2mth.invQ(Q),V))
    26252642        Tx,Ty,Tz = drawingData['viewPoint'][0]
    2626         anglex,angley,anglez,oldxy = drawingData['Rotation']
    2627         if not len(oldxy): oldxy = list(newxy)
    2628         Rx = G2lat.rotdMat(anglex,0)
    2629         Ry = G2lat.rotdMat(angley,1)
    2630         Rz = G2lat.rotdMat(anglez,2)
    2631         dxy = list(newxy-oldxy)+[0,]
    2632         dxy = np.inner(Bmat,np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy))))
    2633         Tx -= dxy[0]*0.01
    2634         Ty += dxy[1]*0.01
    2635         Tz -= dxy[2]*0.01
     2643        Tx += V[0]*0.01
     2644        Ty += V[1]*0.01
     2645        Tz += V[2]*0.01
    26362646        drawingData['Rotation'][3] = list(newxy)
    26372647        drawingData['viewPoint'][0] =  Tx,Ty,Tz
    26382648        SetViewPointText([Tx,Ty,Tz])
    26392649       
    2640     def SetRotation(newxy):       
    2641         anglex,angley,anglez,oldxy = drawingData['Rotation']
     2650    def SetRotation(newxy):
     2651#first get rotation vector in screen coords. & angle increment       
     2652        oldxy = drawingData['Rotation'][3]
    26422653        if not len(oldxy): oldxy = list(newxy)
    26432654        dxy = newxy-oldxy
    2644         anglex += dxy[1]*.25
    2645         angley += dxy[0]*.25
    2646         oldxy = newxy
    2647         drawingData['Rotation'] = [anglex,angley,anglez,list(oldxy)]
     2655        drawingData['Rotation'][3] = list(newxy)
     2656        V = np.array([dxy[1],dxy[0],0.])
     2657        A = 0.25*np.sqrt(dxy[0]**2+dxy[1]**2)
     2658# next transform vector back to xtal coordinates via inverse quaternion
     2659# & make new quaternion
     2660        Q = drawingData['Quaternion']
     2661        V = G2mth.prodQVQ(G2mth.invQ(Q),np.inner(Bmat,V))
     2662        DQ = G2mth.AVdeg2Q(A,V)
     2663        Q = G2mth.prodQQ(Q,DQ)
     2664        drawingData['Quaternion'] = Q
     2665# finally get new view vector - last row of rotation matrix
     2666        VD = np.inner(Bmat,G2mth.Q2Mat(Q)[2])
     2667        VD /= np.sqrt(np.sum(VD**2))
     2668        drawingData['viewDir'] = VD
     2669        SetViewDirText(VD)
    26482670       
    26492671    def SetRotationZ(newxy):                       
    2650         anglex,angley,anglez,oldxy = drawingData['Rotation']
     2672#first get rotation vector (= view vector) in screen coords. & angle increment       
     2673        View = glGetIntegerv(GL_VIEWPORT)
     2674        oldxy = drawingData['Rotation'][3]
    26512675        if not len(oldxy): oldxy = list(newxy)
    26522676        dxy = newxy-oldxy
    2653         anglez += (dxy[0]+dxy[1])*.25
    2654         oldxy = newxy
    2655         drawingData['Rotation'] = [anglex,angley,anglez,list(oldxy)]
    2656        
     2677        drawingData['Rotation'][3] = list(newxy)
     2678        V = drawingData['viewDir']
     2679        X0 = drawingData['viewPoint'][0]
     2680        A = (dxy[0]+dxy[1])*.25
     2681# next transform vector back to xtal coordinates via inverse quaternion
     2682# & make new quaternion
     2683        Q = drawingData['Quaternion']
     2684        V = np.inner(Amat,V)
     2685        DQ = G2mth.AVdeg2Q(A,V)
     2686        Q = G2mth.prodQQ(Q,DQ)
     2687        drawingData['Quaternion'] = Q
     2688
    26572689    def RenderBox():
    26582690        glEnable(GL_COLOR_MATERIAL)
     
    28532885        cPos = drawingData['cameraPos']
    28542886        Zclip = drawingData['Zclip']*cPos/200.
    2855         anglex,angley,anglez = drawingData['Rotation'][:3]
     2887        Q = drawingData['Quaternion']
    28562888        Tx,Ty,Tz = drawingData['viewPoint'][0]
    28572889        cx,ct,cs,ci = drawingData['atomPtrs']
     
    28772909        glMatrixMode(GL_MODELVIEW)
    28782910        glLoadIdentity()
    2879         matX = G2lat.rotdMat(anglex,axis=0)
    2880         matY = G2lat.rotdMat(angley,axis=1)
    2881         matZ = G2lat.rotdMat(anglez,axis=2)
    2882         matRot = np.inner(matX,np.inner(matY,matZ))
     2911        matRot = G2mth.Q2Mat(Q)
    28832912        matRot = np.concatenate((np.concatenate((matRot,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
    28842913        glMultMatrixf(matRot.T)
Note: See TracChangeset for help on using the changeset viewer.