Changeset 3818


Ignore:
Timestamp:
Feb 12, 2019 11:58:33 AM (5 years ago)
Author:
vondreele
Message:

implement sqrt(I) plots for multiplots & contour plots - key 's' to toggle. Old 's' for select color scheme changed to 'S' (capital S).

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImath.py

    r3808 r3818  
    12961296    Mast: array orthogonalization matrix for Uij
    12971297    '''
    1298     ngl = 32
     1298    ngl = 36                    #selected for integer steps for 1/6,1/4,1/3...
    12991299    glTau,glWt = pwd.pygauleg(0.,1.,ngl)         #get Gauss-Legendre intervals & weights
    13001300    Ax = np.array(XSSdata[:3]).T   #atoms x waves x sin pos mods
    … …  
    13491349        MmodA = Am[:,:,:,nxs]*np.sin(twopi*tauM[nxs,:,nxs,:]) #atoms X waves X 3 X ngl
    13501350        MmodB = Bm[:,:,:,nxs]*np.cos(twopi*tauM[nxs,:,nxs,:]) #ditto
    1351         Mmod = np.array([MmodB,MmodA])
    1352         Mmod = np.sum(Mmod,axis=2)                #atoms X 3 X ngl; sum waves
    1353         Mmod = np.swapaxes(Mmod,1,-1)       #ReIm,mxyz,Ntau,Natm
     1351        Mmod = np.sum(MmodA+MmodB,axis=1)
     1352        Mmod = np.swapaxes(Mmod,1,2)    #Mxyz,Ntau,Natm
    13541353    else:
    13551354        Mmod = 1.0

  • trunk/GSASIIplot.py

    r3817 r3818  
    18501850                newPlot = True
    18511851        elif event.key == 's' and 'PWDR' in plottype:
    1852             if G2frame.SinglePlot:  #toggle sqrt plot
    1853                 Page.plotStyle['sqrtPlot'] = not Page.plotStyle['sqrtPlot']
    1854                 if Page.plotStyle['sqrtPlot']:
    1855                     Page.plotStyle['logPlot'] = False
    1856                 Ymax = max(Pattern[1][1])
    1857                 if Page.plotStyle['sqrtPlot']:
    1858                     Page.plotStyle['delOffset'] = .002*np.sqrt(Ymax)
    1859                     Page.plotStyle['refOffset'] = -0.1*np.sqrt(Ymax)
    1860                     Page.plotStyle['refDelt'] = .1*np.sqrt(Ymax)
    1861                 else:
    1862                     Page.plotStyle['delOffset'] = .02*Ymax
    1863                     Page.plotStyle['refOffset'] = -0.1*Ymax
    1864                     Page.plotStyle['refDelt'] = .1*Ymax
    1865             else:   #select color scheme for multiplots & contour plots
    1866                 choice = [m for m in mpl.cm.datad.keys()]   # if not m.endswith("_r")
    1867                 choice.sort()
    1868                 dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
    1869                 if dlg.ShowModal() == wx.ID_OK:
    1870                     sel = dlg.GetSelection()
    1871                     G2frame.ContourColor = choice[sel]
    1872                 else:
    1873                     G2frame.ContourColor = GSASIIpath.GetConfigValue('Contour_color','Paired')
    1874                 dlg.Destroy()
     1852            Page.plotStyle['sqrtPlot'] = not Page.plotStyle['sqrtPlot']
     1853            if Page.plotStyle['sqrtPlot']:
     1854                Page.plotStyle['logPlot'] = False
     1855            Ymax = max(Pattern[1][1])
     1856            if Page.plotStyle['sqrtPlot']:
     1857                Page.plotStyle['delOffset'] = .02*np.sqrt(Ymax)
     1858                Page.plotStyle['refOffset'] = -0.1*np.sqrt(Ymax)
     1859                Page.plotStyle['refDelt'] = .1*np.sqrt(Ymax)
     1860            else:
     1861                Page.plotStyle['delOffset'] = .02*Ymax
     1862                Page.plotStyle['refOffset'] = -0.1*Ymax
     1863                Page.plotStyle['refDelt'] = .1*Ymax
     1864            newPlot = True
     1865        elif event.key == 'S' and 'PWDR' in plottype:
     1866            choice = [m for m in mpl.cm.datad.keys()]   # if not m.endswith("_r")
     1867            choice.sort()
     1868            dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
     1869            if dlg.ShowModal() == wx.ID_OK:
     1870                sel = dlg.GetSelection()
     1871                G2frame.ContourColor = choice[sel]
     1872            else:
     1873                G2frame.ContourColor = GSASIIpath.GetConfigValue('Contour_color','Paired')
     1874            dlg.Destroy()
    18751875            newPlot = True
    18761876        elif event.key == 'u' and (G2frame.Contour or not G2frame.SinglePlot):
    … …  
    19571957                newPlot = True     
    19581958        elif event.key == 'm':
    1959             Page.plotStyle['sqrtPlot'] = False
     1959#            Page.plotStyle['sqrtPlot'] = False
    19601960            if not G2frame.Contour:               
    19611961                G2frame.SinglePlot = not G2frame.SinglePlot               
    … …  
    19951995                    break
    19961996        else:
    1997             print('no binding for key',event.key)
     1997            #print('no binding for key',event.key)
    19981998            #GSASIIpath.IPyBreak()
    19991999            return
    … …  
    25102510            'PWDR' in G2frame.GPXtree.GetItemText(PickId)) and xpos:
    25112511            Id = G2gd.GetGPXtreeItemId(G2frame,PatternId,'Reflection Lists')
    2512 #            GSASIIpath.IPyBreak()
    25132512            if Id:     
    2514                 #Phases = G2frame.GPXtree.GetItemPyData(Id)
    25152513                pick = str(G2frame.itemPicked).split('(',1)[1][:-1]
    25162514                if 'line' not in pick:       #avoid data points, etc.
    … …  
    25742572    if G2frame.Contour:
    25752573        Page.Choice = (' key press','d: lower contour max','u: raise contour max','o: reset contour max','g: toggle grid',
    2576             'i: interpolation method','s: color scheme','c: contour off','t: temperature for y-axis')
     2574            'i: interpolation method','S: color scheme','c: contour off','t: temperature for y-axis','s: toggle sqrt plot')
    25772575    else:
    25782576        if Page.plotStyle['logPlot']:
    … …  
    26052603                    Page.Choice = (' key press','l: offset left','r: offset right','d/D: offset down/10x','u/U: offset up/10x','o: reset offset',
    26062604                        'b: toggle subtract background','n: log(I) on','c: contour on','q: toggle q plot','t: toggle d-spacing plot','g: toggle grid',
    2607                         'm: toggle multidata plot','w: toggle (Io-Ic)/sig plot','f: select data','s: color scheme','+: no selection')
     2605                        'm: toggle multidata plot','w: toggle (Io-Ic)/sig plot','s: toggle sqrt plot','f: select data','S: color scheme','+: no selection')
    26082606            elif plottype in ['SASD','REFD']:
    26092607                if G2frame.SinglePlot:
    … …  
    26982696    if Page.plotStyle['logPlot']:
    26992697        Title = 'log('+Title+')'
     2698    elif Page.plotStyle['sqrtPlot']:
     2699        Title = r'$\sqrt{I}$ for '+Title
     2700        Ymax = np.sqrt(Ymax)
    27002701    if Page.plotStyle['qPlot'] or plottype in ['SASD','REFD'] and not G2frame.Contour:
    27012702        xLabel = r'$Q, \AA^{-1}$'
    … …  
    27882789                Page.plotStyle['logPlot'] or Page.plotStyle['sqrtPlot'] or G2frame.Contour):
    27892790            magLineList = data[0].get('Magnification',[])
    2790             if ('C' in ParmList[0]['Type'][0] and Page.plotStyle['dPlot']
    2791                 ) or ('T' in ParmList[0]['Type'][0] and Page.plotStyle['qPlot']
    2792                 ): # reversed regions relative to data order
     2791            if ('C' in ParmList[0]['Type'][0] and Page.plotStyle['dPlot']) or \
     2792                ('T' in ParmList[0]['Type'][0] and Page.plotStyle['qPlot']): # reversed regions relative to data order
    27932793                tcorner = 1
    27942794                tpos = 1.0
    … …  
    28342834            if Page.plotStyle['sqrtPlot']:
    28352835                olderr = np.seterr(invalid='ignore') #get around sqrt(-ve) error
    2836                 Y = np.where(xye[1]+bxye>=0.,np.sqrt(xye[1]+bxye),-np.sqrt(-xye[1]-bxye))
     2836                Y = np.where(xye[1]+bxye>=0.,np.sqrt(xye[1]+bxye),-np.sqrt(-xye[1]-bxye))+bxye+NoffY*Ymax/100.0
    28372837                np.seterr(invalid=olderr['invalid'])
    28382838            elif 'PWDR' in plottype and G2frame.SinglePlot and not (

  • trunk/GSASIIstrMath.py

    r3812 r3818  
    15111511    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
    15121512
    1513     if parmDict[pfx+'isMag']:       #TODO: fix the math       
    1514 #1st try at this
    1515         GSdata = Gdata[nxs,nxs,:,:]+Mmod    #ReIm,ntau,Mxyz,Natm
    1516         GSdata = np.inner(np.swapaxes(GSdata,2,3),np.swapaxes(SGMT,1,2)).T  #apply sym. ops.--> Mxyz,Nops,Natm,Ntau,ReIm
    1517         if SGData['SGInv'] and not SGData['SGFixed']:
    1518             GSdata = np.hstack((GSdata,-GSdata))       #inversion if any
     1513    if parmDict[pfx+'isMag']:       #TODO: fix the math
     1514        GSdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2))  #apply sym. ops.--> Natm,Nops,Nxyz
     1515        MSmod = np.array([np.roll(Mmod,int(ngl*ssgt[3]),2) for ssgt in SSGT])   #Nops,Natm,Ntau,Mxyz
     1516        if SGData['SGInv'] and not SGData['SGFixed']:   #inversion if any
     1517            GSdata = np.hstack((GSdata,-GSdata))     
     1518            MSmod = np.vstack((MSmod,-MSmod))
    15191519        GSdata = np.hstack([GSdata for icen in range(Ncen)])        #dup over cell centering
    1520         GSdata = SGData['MagMom'][nxs,:,nxs,nxs,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
    1521         Kdata = np.inner(GSdata.T,uAmat).T     #Cartesian unit vectors
    1522         SMag = np.sqrt(np.sum(np.sum(Kdata**2,axis=0),axis=-1))
    1523         Kmean = np.mean(SMag,axis=0)    #normalization --> unit vectors
    1524         Kdata /= Kmean[nxs,nxs,:,:,nxs]      #mxyz,nops,natm,ntau,ReIm
     1520        MSmod = np.vstack([MSmod for icen in range(Ncen)])        #dup over cell centering - right??
     1521       
     1522        GSdata = GSdata[:,:,nxs,:]+np.swapaxes(MSmod,0,1)         #Natm,Nops,Ntau,Mxyz
     1523        GSdata = SGData['MagMom'][nxs,:,nxs,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
     1524       
     1525        Kdata = np.inner(GSdata,uAmat).T     #Cartesian unit vectors
     1526        SMag = np.sqrt(np.sum(Kdata**2,axis=0))
     1527        Kdata /= SMag      #mxyz,ntau,nops,natm
    15251528
    15261529    FF = np.zeros(len(Tdata))
    … …  
    15911594            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
    15921595                     
    1593             D = twopi*H.T[:,3:]*glTau[nxs,:]
    1594             mphase = phase[:,:,nxs,:]+D[:,nxs,:,nxs]
    1595             mphase = np.array([mphase+twopi*np.inner(cen,HP.T)[:,nxs,nxs,nxs] for cen in SGData['SGCen']])
     1596            mphase = np.array([phase+twopi*np.inner(cen,HP.T)[:,nxs,nxs] for cen in SGData['SGCen']])
    15961597            mphase = np.concatenate(mphase,axis=1)    #remove extra axis; Nref,Nop,Natm
    1597             sinm = np.swapaxes(np.sin(mphase),2,3)    #--> Nref,Nop,Natm,Ntau
    1598             cosm = np.swapaxes(np.cos(mphase),2,3)                               #ditto
     1598            sinm = np.sin(mphase)    #--> Nref,Ntau,Nop,Natm
     1599            cosm = np.cos(mphase)                               #ditto
    15991600           
    16001601            HM = np.inner(Bmat,HP.T)                             #put into cartesian space
    16011602            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
    1602             eDotK = np.sum(HM[:,:,nxs,nxs,nxs,nxs]*Kdata[:,nxs,:,:,:,:],axis=0)
    1603             Q = HM[:,:,nxs,nxs,nxs,nxs]*eDotK[nxs,:,:,:,:,:]-Kdata[:,nxs,:,:,:,:] #Mxyz,Nref,Nop,Natm,Ntau,ReIm
     1603            eDotK = np.sum(HM[:,:,nxs,nxs,nxs]*Kdata[:,nxs,:,:,:],axis=0)
     1604            Q = HM[:,:,nxs,nxs,nxs]*eDotK[nxs,:,:,:,:]-Kdata[:,nxs,:,:,:] #Mxyz,Nref,Ntau,Nop,Natm
    16041605
    1605             fam = (Q*TMcorr[nxs,:,nxs,:,nxs,nxs]*cosm[nxs,:,:,:,:,nxs]*SMag[nxs,nxs,:,:,:,nxs])   #Mxyz,Nref,Nop,Natm,Ntau,ReIm
    1606             fbm = (Q*TMcorr[nxs,:,nxs,:,nxs,nxs]*sinm[nxs,:,:,:,:,nxs]*SMag[nxs,nxs,:,:,:,nxs])
     1606            fam = (Q*TMcorr[nxs,:,nxs,nxs,:]*cosm[nxs,:,nxs,:,:]*SMag[nxs,nxs,:,:,:])   #Mxyz,Nref,Ntau,Nop,Natm
     1607            fbm = (Q*TMcorr[nxs,:,nxs,nxs,:]*sinm[nxs,:,nxs,:,:]*SMag[nxs,nxs,:,:,:])
    16071608           
    1608             fams = np.sum(np.sum(np.sum(fam,axis=-1),axis=2),axis=2)      #xyz,Nref,ntau; sum ops & atoms
    1609             fbms = np.sum(np.sum(np.sum(fbm,axis=-1),axis=2),axis=2)      #ditto
     1609            fas = np.sum(np.sum(fam,axis=-1),axis=-1)      #xyz,Nref,ntau; sum ops & atoms
     1610            fbs = np.sum(np.sum(fbm,axis=-1),axis=-1)      #ditto
    16101611           
    1611             fas = np.sum(fams*glWt[nxs,nxs,:],axis=-1)
    1612             fbs = np.sum(fbms*glWt[nxs,nxs,:],axis=-1)
    1613 
     1612            refl.T[10] = np.sum(np.sum(fas**2,axis=0)*glWt[nxs,nxs,:],axis=-1)+     \
     1613                np.sum(np.sum(fbs**2,axis=0)*glWt[nxs,nxs,:],axis=-1)    #square of sums
     1614            refl.T[11] = atan2d(fbs[:,0],fas[:,0])  #ignore f' & f"
     1615           
    16141616        else:
    16151617            GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
    … …  
    16251627            fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
    16261628           
    1627         refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
    1628         refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     1629            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
     1630            refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    16291631        if 'P' not in calcControls[hfx+'histType']:
    16301632            refl.T[8] = np.copy(refl.T[10])               
Note: See TracChangeset for help on using the changeset viewer.