Changeset 1958


Ignore:
Timestamp:
Aug 21, 2015 2:57:01 PM (8 years ago)
Author:
vondreele
Message:

add r-values by SS index
add a fade for frac modulation of atoms in drawing
work on modulation structure factor

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIgrid.py

    r1943 r1958  
    29622962        for value in data[0]:
    29632963            if 'Nref' in value:
    2964                 mainSizer.Add((5,5),)
    29652964                pfx = value.split('Nref')[0]
    29662965                name = data[0].get(pfx.split(':')[0]+'::Name','?')
    2967                 mainSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,' For phase '+name+':'))
    2968                 mainSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,
    2969                     u' Unweighted phase residuals RF\u00b2: %.3f%%, RF: %.3f%% on %d reflections  '% \
    2970                     (data[0][pfx+'Rf^2'],data[0][pfx+'Rf'],data[0][value])))
     2966                if 'SS' in value:
     2967                    mainSizer.Add((5,5),)
     2968                    mainSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,' For incommensurate phase '+name+':'))
     2969                    for m,(Rf2,Rf,Nobs) in enumerate(zip(data[0][pfx+'Rf^2'],data[0][pfx+'Rf'],data[0][value])):
     2970                        mainSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,
     2971                            u' m = +/- %d: RF\u00b2: %.3f%%, RF: %.3f%% on %d reflections  '% \
     2972                            (m,Rf2,Rf,Nobs)))
     2973                else:
     2974                    mainSizer.Add((5,5),)
     2975                    mainSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,' For phase '+name+':'))
     2976                    mainSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,
     2977                        u' Unweighted phase residuals RF\u00b2: %.3f%%, RF: %.3f%% on %d reflections  '% \
     2978                        (data[0][pfx+'Rf^2'],data[0][pfx+'Rf'],data[0][value])))
     2979                   
    29712980    mainSizer.Add((5,5),)
    29722981    mainSizer.Layout()   

  • trunk/GSASIImath.py

    r1957 r1958  
    929929    Smult,TauT = SStauM             # both atoms x SGops
    930930    m = SSUniq.T[3]
    931     nh = np.zeros(1)
     931    A = np.array(XSSdata[:3])   #sin mods x waves x atoms
     932    B = np.array(XSSdata[3:])   #cos mods...
    932933    if XSSdata.ndim > 2:
    933934        nh = np.arange(XSSdata.shape[1])        #0..max harmonic order-1
    934     M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis])     # waves x SGops
    935     A = np.array(XSSdata[:3])   #sin mods x waves x atoms
    936     B = np.array(XSSdata[3:])   #cos mods...
    937     HdotA = (np.inner(A.T,SSUniq.T[:3].T)+SSPhi+TauT[:,np.newaxis,:])    # atoms X waves x SGops
    938     HdotB = (np.inner(B.T,SSUniq.T[:3].T)+SSPhi)    #+TauT[:,np.newaxis,:])    # ditto
    939     GpA = sp.jn(M,twopi*HdotA)                      # atoms x waves x SGops
    940     GpB = sp.jn(M,twopi*HdotB)*(1.j)**M             # ditto
    941     Gp = np.sum(GpA+GpB,axis=1)                     # atoms x SGops
    942     return np.real(Gp).T,np.imag(Gp).T              # SGops x atoms
     935        M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis])     # waves x SGops
     936        HdotA = (np.inner(A.T,SSUniq.T[:3].T))    #atoms X waves x SGops
     937        HdotB = (np.inner(B.T,SSUniq.T[:3].T))    # ditto
     938        Hdot = np.sqrt(HdotA**2+HdotB**2)
     939        Hphi = np.sum(np.arctan2(HdotB,HdotA)*M,axis=1)     
     940        GpA = np.sum(sp.jn(-M,twopi*Hdot),axis=1)                      # sum on waves
     941    else:
     942        HdotA = np.inner(A.T,SSUniq.T[:3].T)    #atoms X SGops
     943        HdotB = np.inner(B.T,SSUniq.T[:3].T)    # ditto
     944        Hdot = np.sqrt(HdotA**2+HdotB**2)
     945        Hphi = np.arctan2(HdotB,HdotA)       
     946        GpA = sp.jn(-m,twopi*Hdot)                      # atoms x SGops
     947       
     948#    GSASIIpath.IPyBreak()
     949    return GpA.T,Hphi.T             # SGops x atoms
    943950   
    944951def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM):
    945952    import scipy.special as sp
    946     m = SSUniq.T[3]
     953    Smult,TauT = SStauM             # both atoms x SGops
    947954    nh = np.zeros(1)
    948955    if XSSdata.ndim > 2:
    … …  
    950957    M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis])
    951958    A = np.array([[a,b] for a,b in zip(XSSdata[:3],XSSdata[3:])])
    952     HdotA = twopi*(np.inner(SSUniq.T[:3].T,A.T)+SSPhi)
     959    HdotA = twopi*(np.inner(A.T,SSUniq.T[:3].T)+SSPhi+TauT[:,np.newaxis,:])
    953960    Gpm = sp.jn(M-1,HdotA)
    954961    Gpp = sp.jn(M+1,HdotA)
    … …  
    9951002    atoms = data['Atoms']
    9961003    drawAtoms = drawingData['Atoms']
     1004    Fade = np.zeros(len(drawAtoms))
    9971005    for atom in atoms:   
    9981006        atxyz = G2spc.MoveToUnitCell(np.array(atom[cx:cx+3]))[0]
    9991007        atuij = np.array(atom[cia+2:cia+8])
    10001008        waveType = atom[-1]['SS1']['waveType']
     1009        Sfrac = atom[-1]['SS1']['Sfrac']
    10011010        Spos = atom[-1]['SS1']['Spos']
    10021011        Sadp = atom[-1]['SS1']['Sadp']
    … …  
    10071016            sop,ssop,icent = G2spc.OpsfromStringOps(opr,SGData,SSGData)
    10081017            sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(tau,sop,ssop,atxyz)
    1009             smul = sdet*ssdet
    1010 #            tauT *= icent
     1018            smul = sdet         # n-fold vs m operator on wave
     1019            tauT *= icent       #invert wave on -1
    10111020            wave = np.zeros(3)
    10121021            uwave = np.zeros(6)
     1022            #how do I handle Sfrac? - fade the atoms?
     1023            if len(Sfrac):
     1024                scof = []
     1025                ccof = []
     1026                for i,sfrac in enumerate(Sfrac):
     1027                    if not i and 'Crenel' in waveType:
     1028                        Fade[ind] += fracCrenel(tauT,sfrac[0][0],sfrac[0][1])
     1029                    else:
     1030                        scof.append(sfrac[0][0])
     1031                        ccof.append(sfrac[0][1])
     1032                    if len(scof):
     1033                        Fade[ind] += np.sum(fracFourier(tauT,scof,ccof))                           
    10131034            if len(Spos):
    10141035                scof = []
    … …  
    10401061                X = G2spc.ApplyStringOps(opr,SGData,atxyz+wave)
    10411062                drawatom[dcx:dcx+3] = X
    1042     return drawAtoms
     1063    return drawAtoms,Fade
    10431064   
    10441065   

  • trunk/GSASIIphsGUI.py

    r1957 r1958  
    24422442                    waveSizer.Add(wx.StaticText(waveData,label=' %s  parameters: %s'%(waveName,str(names).rstrip(']').lstrip('[').replace("'",''))),0,WACV)
    24432443                    for ival,val in enumerate(wave[0]):
    2444                         if any(CSI[Stype][0][ival]):
     2444                        if np.any(CSI[Stype][0][ival]):
    24452445                            waveVal = wx.TextCtrl(waveData,value='%.4f'%(val),style=wx.TE_PROCESS_ENTER)
    24462446                            waveVal.Bind(wx.EVT_TEXT_ENTER,OnWaveVal)

  • trunk/GSASIIplot.py

    r1957 r1958  
    30313031                scof.append(spos[0][:3])
    30323032                ccof.append(spos[0][3:])
    3033         wave += G2mth.posFourier(tau,np.array(scof),np.array(ccof),-1)  #hmm, why -1 work for Na2CO3?
     3033        wave += G2mth.posFourier(tau,np.array(scof),np.array(ccof),1)  #hmm, why -1 work for Na2CO3?
    30343034    if mapData['Flip']:
    30353035        Title = 'Charge flip'
    … …  
    44144414            if keyBox:
    44154415                OnKeyPressed(event)
    4416         Draw('key')
     4416            return
     4417        Draw('key up')
    44174418       
    44184419    def OnKeyPressed(event):    #On key down for repeating operation - used to change tau...
    … …  
    44334434            G2frame.tau %= 1.   #force 0-1 range
    44344435            G2frame.G2plotNB.status.SetStatusText('Modulation tau = %.2f'%(G2frame.tau),1)
    4435             data['Drawing']['Atoms'] = G2mth.ApplyModulation(data,G2frame.tau)     #modifies drawing atom array!         
     4436            data['Drawing']['Atoms'],Fade = G2mth.ApplyModulation(data,G2frame.tau)     #modifies drawing atom array!         
    44364437            SetDrawAtomsText(data['Drawing']['Atoms'])
    44374438            G2phG.FindBondsDraw(data)           #rebuild bonds & polygons
    4438             Draw('key')
     4439            if not np.any(Fade):
     4440                Fade += 1
     4441            Draw('key down',Fade)
    44394442           
    44404443    def GetTruePosition(xy,Add=False):
    … …  
    50015004        glShadeModel(GL_SMOOTH)
    50025005                           
    5003     def Draw(caller=''):
     5006    def Draw(caller='',Fade=None):
    50045007#useful debug?       
    50055008#        if caller:
    … …  
    50755078#        glEnable(GL_BLEND)
    50765079#        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
     5080        if Fade == None:
     5081            atmFade = np.ones(len(drawingData['Atoms']))
     5082        else:
     5083            atmFade = Fade
    50775084        for iat,atom in enumerate(drawingData['Atoms']):
    50785085            x,y,z = atom[cx:cx+3]
    … …  
    50845091                atNum = -1
    50855092            CL = atom[cs+2]
    5086             atColor = np.array(CL)/255.
     5093            if not atmFade[iat]:
     5094                continue
     5095            atColor = atmFade[iat]*np.array(CL)/255.
    50875096            if drawingData['showRigidBodies'] and atom[ci] in rbAtmDict:
    50885097                bndColor = Or

  • trunk/GSASIIspc.py

    r1957 r1958  
    15061506        CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
    15071507        if 'Crenel' in waveType:
    1508             dF = fracCrenel(tau,delt2[:1],delt2[1:]).squeeze()
     1508            dF = np.zeros_like(tau)
    15091509        else:
    15101510            dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
    … …  
    15171517            fsc = np.ones(2,dtype='i')
    15181518            if 'Crenel' in waveType:
    1519                 dFT = fracCrenel(tauT,delt2[:1],delt2[1:]).squeeze()
     1519                dFT = np.zeros_like(tau)
    15201520                fsc = [1,1]
    15211521            else:   #Fourier
    … …  
    15711571            sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
    15721572            xsc = np.ones(6,dtype='i')
    1573             if waveType == 'Fourier':
     1573            if 'Fourier' in waveType:
    15741574                dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,np.newaxis,np.newaxis]
    15751575            elif waveType == 'Sawtooth':

  • trunk/GSASIIstrMath.py

    r1957 r1958  
    640640        for isym,(sop,ssop) in enumerate(zip(SGOps,SSOps)):
    641641            sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(0,sop,ssop,xyz)
    642             Smult[ix][isym] = sdet*ssdet
     642            Smult[ix][isym] = sdet
    643643            TauT[ix][isym] = tauT
    644644    return Smult,TauT
    … …  
    971971   
    972972    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     973    SGInv = SGData['SGInv']
    973974    SGT = np.array([ops[1] for ops in SGData['SGOps']])
    974975    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
    … …  
    978979    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
    979980    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    980     SStauM = GetSSTauM(SGData['SGOps'],SSGData['SSGOps'],pfx,calcControls,Xdata)
     981    SStauM = list(GetSSTauM(SGData['SGOps'],SSGData['SSGOps'],pfx,calcControls,Xdata))
     982    if SGData['SGInv']:
     983        SStauM[0] = np.hstack((SStauM[0],SStauM[0]))
     984        SStauM[1] = np.hstack((SStauM[1],SStauM[1]))
     985    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
    981986    FF = np.zeros(len(Tdata))
    982987    if 'NC' in calcControls[hfx+'histType']:
    … …  
    10111016        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    10121017        FF = refDict['FF']['FF'][iref][Tindx]
    1013         Uniq = np.inner(H[:3],SGMT)
    1014         SSUniq = np.inner(H,SSGMT)
    1015         Phi = np.inner(H[:3],SGT)
    1016         SSPhi = np.inner(H,SSGT)
    1017         GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM)       
     1018        HM = np.zeros(4)
     1019        HM[:3] += H[3]*modQ
     1020        HM += H
     1021        Uniq = np.inner(HM[:3],SGMT)
     1022        SSUniq = np.inner(HM,SSGMT)
     1023        Phi = np.inner(HM[:3],SGT)
     1024        SSPhi = np.inner(HM,SSGT)
     1025        if SGInv:   #if centro - expand HKL sets
     1026            Uniq = np.vstack((Uniq,-Uniq))
     1027            SSUniq = np.vstack((SSUniq,-SSUniq))
     1028            Phi = np.hstack((Phi,-Phi))
     1029            SSPhi = np.hstack((SSPhi,-SSPhi))
     1030        GfpuA,Hphi = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM)       
    10181031        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
    10191032        sinp = np.sin(phase)
    … …  
    10251038        Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
    10261039        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])     #2 x sym x atoms
    1027         fb = np.zeros_like(fa)                                      #ditto
    1028         if not SGData['SGInv']:
    1029             fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
    1030         fa = fa*GfpuA[np.newaxis,:,:]-fb*GfpuB[np.newaxis,:,:]
    1031         fb = fb*GfpuA[np.newaxis,:,:]+fa*GfpuB[np.newaxis,:,:]
    1032         fas = np.real(np.sum(np.sum(fa,axis=1),axis=1))        #real
     1040        fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
     1041        fa *= GfpuA
     1042        fb *= GfpuA       
     1043        fas = np.real(np.sum(np.sum(fa,axis=1),axis=1))
    10331044        fbs = np.real(np.sum(np.sum(fb,axis=1),axis=1))
    10341045        fasq = fas**2
    … …  
    10371048        refl[7+im] = np.sum(fasq)+np.sum(fbsq)
    10381049        refl[10+im] = atan2d(fbs[0],fas[0])
     1050#        GSASIIpath.IPyBreak()
    10391051   
    10401052def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
    … …  
    27722784            sumdF = 0
    27732785            sumdF2 = 0
     2786            if im:
     2787                sumSSFo = np.zeros(10)
     2788                sumSSFo2 = np.zeros(10)
     2789                sumSSdF = np.zeros(10)
     2790                sumSSdF2 = np.zeros(10)
     2791                sumSSwYo = np.zeros(10)
     2792                sumSSwdf2 = np.zeros(10)
     2793                SSnobs = np.zeros(10)
    27742794            nobs = 0
    27752795            nrej = 0
    27762796            next = 0
     2797            maxH = 0
    27772798            if calcControls['F**2']:
    27782799                for i,ref in enumerate(refDict['RefList']):
    … …  
    27922813                            df[i] = -w*(ref[5+im]-ref[7+im])
    27932814                            sumwYo += (w*ref[5+im])**2      #w*Fo^2
     2815                            if im:  #accumulate super lattice sums
     2816                                ind = int(abs(ref[3]))
     2817                                sumSSFo[ind] += Fo
     2818                                sumSSFo2[ind] += ref[5+im]
     2819                                sumSSdF[ind] += abs(Fo-np.sqrt(ref[7+im]))
     2820                                sumSSdF2[ind] += abs(ref[5+im]-ref[7+im])
     2821                                sumSSwYo[ind] += (w*ref[5+im])**2      #w*Fo^2
     2822                                sumSSwdf2[ind] +=  df[i]**2
     2823                                SSnobs[ind] += 1
     2824                                maxH = max(maxH,ind)                           
    27942825                        else:
    27952826                            if ref[3+im]:
    … …  
    28162847                            df[i] = -w*(Fo-Fc)
    28172848                            sumwYo += (w*Fo)**2
     2849                            if im:
     2850                                ind = int(abs(ref[3]))
     2851                                sumSSFo[ind] += Fo
     2852                                sumSSFo2[ind] += ref[5+im]
     2853                                sumSSdF[ind] += abs(Fo-Fc)
     2854                                sumSSdF2[ind] += abs(ref[5+im]-ref[7+im])
     2855                                sumSSwYo[ind] += (w*Fo)**2                                                           
     2856                                sumSSwdf2[ind] +=  df[i]**2
     2857                                SSnobs[ind] += 1                           
     2858                                maxH = max(maxH,ind)                           
    28182859                        else:
    28192860                            if ref[3+im]:
    … …  
    28312872            Histogram['Residuals'][phfx+'Nrej'] = nrej
    28322873            Histogram['Residuals'][phfx+'Next'] = next
     2874            if im:
     2875                Histogram['Residuals'][phfx+'SSRf'] = 100.*sumSSdF[:maxH+1]/sumSSFo[:maxH+1]
     2876                Histogram['Residuals'][phfx+'SSRf^2'] = 100.*sumSSdF2[:maxH+1]/sumSSFo2[:maxH+1]
     2877                Histogram['Residuals'][phfx+'SSNref'] = SSnobs[:maxH+1]
     2878                Histogram['Residuals']['SSwR'] = np.sqrt(sumSSwdf2[:maxH+1]/sumSSwYo[:maxH+1])*100.               
    28332879            Nobs += nobs
    28342880            Nrej += nrej
Note: See TracChangeset for help on using the changeset viewer.