Changeset 3320


Ignore:
Timestamp:
Mar 23, 2018 1:41:13 PM (4 years ago)
Author:
vondreele
Message:

implement new delt/sig subplot for PWDR & SASD plots - nvoked with 'w' key
make new MagStructureFactor2 routine - removed magnetism stuff from StructureFactor2

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIplot.py

    r3316 r3320  
    425425            new = False
    426426            plotNum,Page = self.GetTabIndex(label)
    427             if Type == 'mpl' or Type == '3d':         
     427            if Type == 'mpl' or Type == '3d':
     428                Axes = Page.figure.get_axes()
    428429                Plot = Page.figure.gca()          #get previous plot
    429                 limits = Plot.get_xlim(),Plot.get_ylim() # save previous limits
    430 #                print 'Plot limits:',limits
     430                limits = [Plot.get_xlim(),Plot.get_ylim()] # save previous limits
     431                if len(Axes)>1:
     432                    limits[1] = Axes[1].get_ylim()
     433#                    print('Axes[1]',Axes[1].get_ylim())
     434#                print ('Plot limits:',limits,Axes)
    431435                if newImage:
    432436                    Page.figure.clf()
     
    581585            deleted = False
    582586        Toolbar.__init__(self,plotCanvas)
    583         G2path = os.path.split(os.path.abspath(__file__))[0]
     587#        G2path = os.path.split(os.path.abspath(__file__))[0]
    584588        self.updateActions = None # defines a call to be made as part of plot updates
    585589        self.plotCanvas = plotCanvas
     
    599603            fil = ''.join([i[0].lower() for i in direc.split()]+['arrow.ico'])
    600604            self.arrows[direc] = self.AddToolBarTool(sprfx+direc,prfx+direc,fil,self.OnArrow)
    601         G2path = os.path.split(os.path.abspath(__file__))[0]
     605#        G2path = os.path.split(os.path.abspath(__file__))[0]
    602606        if publish:
    603607            self.AddToolBarTool('Publish plot','Create publishable version of plot','publish.ico',publish)
     
    14941498            return
    14951499        newPlot = False
    1496         if event.key == 'w':
     1500        if event.key == 'w' and not G2frame.plotStyle['qPlot'] and not G2frame.plotStyle['dPlot']:  #can't do weight plots when x-axis is different
    14971501            G2frame.Weight = not G2frame.Weight
    14981502            if not G2frame.Weight and 'PWDR' in plottype:
     
    16011605                G2frame.plotStyle['qPlot'] = not G2frame.plotStyle['qPlot']
    16021606                if G2frame.plotStyle['qPlot']:
     1607                    G2frame.Weight = False
    16031608                    G2frame.Contour = False
    16041609                G2frame.plotStyle['dPlot'] = False
     
    16091614            if G2frame.plotStyle['dPlot']:
    16101615                G2frame.Contour = False               
     1616                G2frame.Weight = False
    16111617            G2frame.plotStyle['qPlot'] = False
    16121618            newPlot = True     
     
    21742180            if msg: msg += '\n'
    21752181            msg += " * only when the intensity scale is linear (not log or sqrt)"
     2182        if G2frame.Weight:
     2183            if msg: msg += '\n'
     2184            msg += " * only when weight plot is set to no weight plot"
    21762185        if msg:
    21772186            msg = 'Publication export is only available under limited plot settings\n'+msg
     
    21872196    else:
    21882197        publish = None
    2189     new,plotNum,Page,Plot,limits = G2frame.G2plotNB.FindPlotTab('Powder Patterns','mpl',
    2190                                                                 publish=publish)
     2198    new,plotNum,Page,Plot,limits = G2frame.G2plotNB.FindPlotTab('Powder Patterns','mpl',publish=publish)
    21912199    if not new:
    21922200        G2frame.xylim = limits
     
    22362244                    Page.Choice = (' key press','n: log(I) off',
    22372245                        'c: contour on','q: toggle q plot','t: toggle d-spacing plot',
    2238                             'm: toggle multidata plot','w: toggle divide by sig','+: toggle selection')
     2246                            'm: toggle multidata plot','+: toggle selection')
    22392247                else:
    22402248                    Page.Choice = (' key press','n: log(I) off',
    22412249                        'd: offset down','l: offset left','r: offset right','u: offset up','o: reset offset',
    22422250                        'c: contour on','q: toggle q plot','t: toggle d-spacing plot','f: select data',
    2243                         'm: toggle multidata plot','w: toggle divide by sig','+: toggle selection')
     2251                        'm: toggle multidata plot','+: toggle selection')
    22442252            elif plottype in ['SASD','REFD']:
    22452253                if G2frame.SinglePlot:
     
    22562264                        'b: toggle subtract background','n: log(I) on','s: toggle sqrt plot','c: contour on',
    22572265                        'q: toggle q plot','t: toggle d-spacing plot','m: toggle multidata plot',
    2258                         'w: toggle divide by sig','+: no selection')
     2266                        'w: toggle (Io-Ic)/sig plot','+: no selection')
    22592267                else:
    22602268                    Page.Choice = (' key press','l: offset left','r: offset right','d/D: offset down/10x','u/U: offset up/10x','o: reset offset',
    22612269                        'b: toggle subtract background','n: log(I) on','c: contour on','q: toggle q plot','t: toggle d-spacing plot',
    2262                         'm: toggle multidata plot','f: select data','s: color scheme','w: toggle divide by sig','+: no selection')
     2270                        'm: toggle multidata plot','w: toggle (Io-Ic)/sig plot','f: select data','s: color scheme','+: no selection')
    22632271            elif plottype in ['SASD','REFD']:
    22642272                if G2frame.SinglePlot:
     
    23362344    #Plot.set_title(Title) # show title only w/o magnification
    23372345    if G2frame.plotStyle['qPlot'] or plottype in ['SASD','REFD'] and not G2frame.Contour:
    2338         Plot.set_xlabel(r'$Q, \AA^{-1}$',fontsize=16)
     2346        xLabel = r'$Q, \AA^{-1}$'
    23392347    elif G2frame.plotStyle['dPlot'] and 'PWDR' in plottype and not G2frame.Contour:
    2340         Plot.set_xlabel(r'$d, \AA$',fontsize=16)
     2348        xLabel = r'$d, \AA$'
    23412349    else:
    2342         if 'C' in ParmList[0]['Type'][0]:       
    2343             Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=16)
     2350        if 'C' in ParmList[0]['Type'][0]:
     2351            xLabel = r'$\mathsf{2\theta}$'
    23442352        else:
    23452353            if G2frame.Contour:
    2346                 Plot.set_xlabel(r'Channel no.',fontsize=16)           
     2354                xLabel = r'Channel no.'
    23472355            else:
    2348                 Plot.set_xlabel(r'$TOF, \mathsf{\mu}$s',fontsize=16)           
     2356                xLabel = r'$TOF, \mathsf{\mu}$s'
    23492357    if G2frame.Weight:
     2358        Plot.tick_params('x',length=0,labelbottom=False)
     2359        Plot.tick_params('y',length=0,labelleft=False)
     2360        GS_kw = {'height_ratios':[4, 1],}
     2361        Plot,Plot1 = Page.figure.subplots(2,1,sharex=True,gridspec_kw=GS_kw)
     2362        Plot1.set_ylabel(r'$\mathsf{\Delta(I)/\sigma(I)}$',fontsize=16)
     2363        Plot1.set_xlabel(xLabel,fontsize=16)
     2364        Page.figure.subplots_adjust(left=16/100.,bottom=16/150.,
     2365            right=.98,top=1.-16/200.,hspace=0)
     2366    else:
     2367        Plot.set_xlabel(xLabel,fontsize=16)
     2368    if 'C' in ParmList[0]['Type'][0]:
    23502369        if 'PWDR' in plottype:
    2351             Plot.set_ylabel(r'$\mathsf{I/\sigma(I)}$',fontsize=16)
    2352         elif plottype in ['SASD','REFD']:
    2353             Plot.set_ylabel(r'$\mathsf{\Delta(I)/\sigma(I)}$',fontsize=16)
    2354     else:
    2355         if 'C' in ParmList[0]['Type'][0]:
    2356             if 'PWDR' in plottype:
    2357                 if G2frame.plotStyle['sqrtPlot']:
    2358                     Plot.set_ylabel(r'$\sqrt{Intensity}$',fontsize=16)
    2359                 else:
    2360                     Plot.set_ylabel(r'$Intensity$',fontsize=16)
    2361             elif plottype == 'SASD':
    2362                 if G2frame.plotStyle['sqPlot']:
    2363                     Plot.set_ylabel(r'$S(Q)=I*Q^{4}$',fontsize=16)
    2364                 else:
    2365                     Plot.set_ylabel(r'$Intensity,\ cm^{-1}$',fontsize=16)
    2366             elif plottype == 'REFD':
    2367                 if G2frame.plotStyle['sqPlot']:
    2368                     Plot.set_ylabel(r'$S(Q)=R*Q^{4}$',fontsize=16)
    2369                 else:
    2370                     Plot.set_ylabel(r'$Reflectivity$',fontsize=16)               
    2371         else:       #neutron TOF
    23722370            if G2frame.plotStyle['sqrtPlot']:
    2373                 Plot.set_ylabel(r'$\sqrt{Normalized\ intensity}$',fontsize=16)
     2371                Plot.set_ylabel(r'$\sqrt{Intensity}$',fontsize=16)
    23742372            else:
    2375                 Plot.set_ylabel(r'$Normalized\ intensity$',fontsize=16)
     2373                Plot.set_ylabel(r'$Intensity$',fontsize=16)
     2374        elif plottype == 'SASD':
     2375            if G2frame.plotStyle['sqPlot']:
     2376                Plot.set_ylabel(r'$S(Q)=I*Q^{4}$',fontsize=16)
     2377            else:
     2378                Plot.set_ylabel(r'$Intensity,\ cm^{-1}$',fontsize=16)
     2379        elif plottype == 'REFD':
     2380            if G2frame.plotStyle['sqPlot']:
     2381                Plot.set_ylabel(r'$S(Q)=R*Q^{4}$',fontsize=16)
     2382            else:
     2383                Plot.set_ylabel(r'$Reflectivity$',fontsize=16)               
     2384    else:       #neutron TOF
     2385        if G2frame.plotStyle['sqrtPlot']:
     2386            Plot.set_ylabel(r'$\sqrt{Normalized\ intensity}$',fontsize=16)
     2387        else:
     2388            Plot.set_ylabel(r'$Normalized\ intensity$',fontsize=16)
    23762389    mpl.rcParams['image.cmap'] = G2frame.ContourColor
    23772390    mcolors = mpl.cm.ScalarMappable()       #wants only default as defined in previous line!!
     
    25552568                    else:
    25562569                        Plot.set_ylim(bottom=np.min(np.trim_zeros(YB))/2.,top=np.max(Y)*2.)
     2570                if G2frame.Weight:
     2571                    Ibeg = np.searchsorted(X,limits[1][0])
     2572                    Ifin = np.searchsorted(X,limits[1][1])
     2573                    Plot1.set_yscale("linear")                                                 
     2574                    DZ = (xye[1]-xye[3])*np.sqrt(xye[2])
     2575                    DifLine = Plot1.plot(X[Ibeg:Ifin],DZ[Ibeg:Ifin],colors[3],picker=1.,label='_diff')                    #(Io-Ic)/sig(Io)
     2576                    Plot1.axhline(0.,color='k')
     2577                    Plot1.set_ylim(bottom=np.min(DZ[Ibeg:Ifin])*1.2,top=np.max(DZ[Ibeg:Ifin])*1.2) 
    25572578                if G2frame.logPlot:
    25582579                    if 'PWDR' in plottype:
     
    25652586                        Ibeg = np.searchsorted(X,limits[1][0])
    25662587                        Ifin = np.searchsorted(X,limits[1][1])
    2567                         if G2frame.Weight:
    2568                             Plot.set_yscale("linear")
    2569                             DS = (YB-ZB)*np.sqrt(xye[2])
    2570                             Plot.plot(X[Ibeg:Ifin],DS[Ibeg:Ifin],colors[3],picker=False,label='_diff')
    2571                             Plot.axhline(0.,color='k')
    2572                             Plot.set_ylim(bottom=np.min(DS[Ibeg:Ifin])*1.2,top=np.max(DS[Ibeg:Ifin])*1.2)                                                   
     2588                        Plot.set_yscale("log",nonposy='mask')
     2589                        if G2frame.ErrorBars:
     2590                            if G2frame.plotStyle['sqPlot']:
     2591                                Plot.errorbar(X,YB,yerr=X**4*Sample['Scale'][0]*np.sqrt(1./(Pattern[0]['wtFactor']*xye[2])),
     2592                                    ecolor=colors[0],picker=3.,clip_on=Clip_on)
     2593                            else:
     2594                                Plot.errorbar(X,YB,yerr=Sample['Scale'][0]*np.sqrt(1./(Pattern[0]['wtFactor']*xye[2])),
     2595                                    ecolor=colors[0],picker=3.,clip_on=Clip_on,label='_obs')
    25732596                        else:
    2574                             Plot.set_yscale("log",nonposy='mask')
    2575                             if G2frame.ErrorBars:
    2576                                 if G2frame.plotStyle['sqPlot']:
    2577                                     Plot.errorbar(X,YB,yerr=X**4*Sample['Scale'][0]*np.sqrt(1./(Pattern[0]['wtFactor']*xye[2])),
    2578                                         ecolor=colors[0],picker=3.,clip_on=Clip_on)
    2579                                 else:
    2580                                     Plot.errorbar(X,YB,yerr=Sample['Scale'][0]*np.sqrt(1./(Pattern[0]['wtFactor']*xye[2])),
    2581                                         ecolor=colors[0],picker=3.,clip_on=Clip_on,label='_obs')
    2582                             else:
    2583                                 Plot.plot(X,YB,colors[0]+pP,picker=3.,clip_on=Clip_on,label='_obs')
    2584                             Plot.plot(X,W,colors[2],picker=False,label='_bkg')     #const. background
    2585                             Plot.plot(X,ZB,colors[1],picker=False,label='_calc')
    2586                 elif G2frame.Weight and 'PWDR' in plottype:
    2587                     DY = xye[1]*np.sqrt(xye[2])
    2588                     Ymax = max(DY)
    2589                     DZ = xye[3]*np.sqrt(xye[2])
    2590                     DS = xye[5]*np.sqrt(xye[2])-Ymax*Pattern[0]['delOffset']
    2591                     ObsLine = Plot.plot(X,DY,colors[0]+pP,picker=3.,clip_on=Clip_on,label='_obs')         #Io/sig(Io)
    2592                     Plot.plot(X,DZ,colors[1],picker=False,label='_calc')                    #Ic/sig(Io)
    2593                     DifLine = Plot.plot(X,DS,colors[3],picker=1.,label='_diff')                    #(Io-Ic)/sig(Io)
    2594                     Plot.axhline(0.,color='k')
     2597                            Plot.plot(X,YB,colors[0]+pP,picker=3.,clip_on=Clip_on,label='_obs')
     2598                        Plot.plot(X,W,colors[2],picker=False,label='_bkg')     #const. background
     2599                        Plot.plot(X,ZB,colors[1],picker=False,label='_calc')
    25952600                else:
    25962601                    if G2frame.SubBack:
     
    26082613                            Plot.plot(X,YB,colors[0]+pP,picker=3.,clip_on=Clip_on,label='_obs')
    26092614                            Plot.plot(X,ZB,colors[1],picker=False,label='_calc')
    2610                     if 'PWDR' in plottype:
     2615                    if 'PWDR' in plottype: 
    26112616                        Plot.plot(X,W,colors[2],picker=False,label='_bkg')                 #Ib
    2612                         DifLine = Plot.plot(X,D,colors[3],picker=1.,label='_diff')                 #Io-Ic
     2617                        if not G2frame.Weight: DifLine = Plot.plot(X,D,colors[3],picker=1.,label='_diff')                 #Io-Ic
    26132618                    Plot.axhline(0.,color='k',label='_zero')
    26142619                Page.SetToolTipString('')
     
    29202925        '''
    29212926        hcfigure = mpl.figure.Figure(dpi=plotOpt['dpi'],figsize=(plotOpt['width'],plotOpt['height']))
    2922         hccanvas = hcCanvas(hcfigure)
     2927#        hccanvas = hcCanvas(hcfigure)
    29232928        CopyRietveldPlot(G2frame,Pattern,Plot,Page,hcfigure)
    29242929        longFormatName,typ = plotOpt['format'].split(',')
  • trunk/GSASIIspc.py

    r3317 r3320  
    18821882    phi = f[:Nuniq]
    18831883    return iabsnt,mulp,Uniq,phi
     1884
     1885def MagHKLchk(HKL,SGData):
     1886    SpnFlp = SGData['SpnFlp']
     1887    print(HKL)
     1888    Uniq = GenHKL(HKL,SGData)
    18841889   
    18851890def checkSSLaue(HKL,SGData,SSGData):
     
    29632968#    print('dupDir',dupDir)
    29642969    SplitSytSym = SytSym.split('(')
     2970    if SGData['SGGray']:
     2971        return SytSym+"1'"
    29652972    if SytSym == '1':       #genersl position
    29662973        return SytSym
  • trunk/GSASIIstrMath.py

    r3297 r3320  
    666666    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
    667667    SGT = np.array([ops[1] for ops in SGData['SGOps']])
    668     Ncen = len(SGData['SGCen'])
    669     Nops = len(SGMT)*Ncen*(1+SGData['SGInv'])
    670668    FFtables = calcControls['FFtables']
    671669    BLtables = calcControls['BLtables']
    672     MFtables = calcControls['MFtables']
    673670    Amat,Bmat = G2lat.Gmat2AB(G)
    674671    Flack = 1.0
     
    687684    if not Xdata.size:          #no atoms in phase!
    688685        return
    689     if parmDict[pfx+'isMag']:       #TODO: fix the math - mag moments now along crystal axes
    690         Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
    691         Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
    692         Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
    693         if SGData['SGInv'] and not SGData['SGFixed']:
    694             Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
    695         Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
    696         Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
    697         Mag = np.tile(Mag[:,nxs],len(SGMT)*Ncen).T
    698         if SGData['SGInv'] and not SGData['SGFixed']:
    699             Mag = np.repeat(Mag,2,axis=0)                  #Mag same shape as Gdata
    700686    if 'NC' in calcControls[hfx+'histType']:
    701687        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     
    712698        refDict['FF']['El'] = list(dat.keys())
    713699        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
    714         refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
    715         for iel,El in enumerate(refDict['FF']['El']):
    716             if El in MFtables:
    717                 refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
    718700    else:       #'X'
    719701        dat = G2el.getFFvalues(FFtables,0.)
     
    759741        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    760742        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
    761         if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:       #TODO: math here??
    762             MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
    763             TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
    764             if SGData['SGInv'] and not SGData['SGFixed']:
    765                 mphase = np.hstack((phase,-phase))
    766             else:
    767                 mphase = phase
    768             mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
    769             mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
    770             sinm = np.sin(mphase)                               #ditto - match magstrfc.for
    771             cosm = np.cos(mphase)                               #ditto
    772             HM = np.inner(Bmat.T,H)                             #put into cartesian space
    773             HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
    774             eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
    775             Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
    776             fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
    777             fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
    778             fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
    779             fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
    780             refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
    781             refl.T[7] = np.copy(refl.T[9])               
    782             refl.T[10] = 0.0    #atan2d(fbs[0],fas[0]) - what is phase for mag refl?
     743        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
     744        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
     745            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
     746            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
    783747        else:
    784             Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
    785             if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
    786                 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
    787                 fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
    788             else:
    789                 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    790                 fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
    791             fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
    792             fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
    793             if SGData['SGInv']: #centrosymmetric; B=0
    794                 fbs[0] *= 0.
    795                 fas[1] *= 0.
    796             if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
    797                 refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
     748            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
     749            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
     750        fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
     751        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
     752        if SGData['SGInv']: #centrosymmetric; B=0
     753            fbs[0] *= 0.
     754            fas[1] *= 0.
     755        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
     756            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
     757            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     758        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
     759            if len(TwinLaw) > 1:
     760                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
     761                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
     762                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
     763                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
     764            else:   # checked correct!!
     765                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
     766                refl.T[7] = np.copy(refl.T[9])               
    798767                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    799                 if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
    800                     refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
    801             else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
    802                 if len(TwinLaw) > 1:
    803                     refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
    804                     refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
    805                         np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
    806                     refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
    807                 else:   # checked correct!!
    808                     refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
    809                     refl.T[7] = np.copy(refl.T[9])               
    810                     refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    811768#                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
    812769        iBeg += blkSize
    813770#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
    814771   
     772def MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     773    ''' Compute neutron magnetic structure factors for all h,k,l for phase
     774    puts the result, F^2, in each ref[8] in refList
     775    operates on blocks of 100 reflections for speed
     776    input:
     777   
     778    :param dict refDict: where
     779        'RefList' list where each ref = h,k,l,it,d,...
     780        'FF' dict of form factors - filed in below
     781    :param np.array G:      reciprocal metric tensor
     782    :param str pfx:    phase id string
     783    :param dict SGData: space group info. dictionary output from SpcGroup
     784    :param dict calcControls:
     785    :param dict ParmDict:
     786
     787    '''       
     788#    phfx = pfx.split(':')[0]+hfx
     789    ast = np.sqrt(np.diag(G))
     790    Mast = twopisq*np.multiply.outer(ast,ast)
     791    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     792    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     793    Ncen = len(SGData['SGCen'])
     794    Nops = len(SGMT)*Ncen*(1+SGData['SGInv'])
     795    MFtables = calcControls['MFtables']
     796    Amat,Bmat = G2lat.Gmat2AB(G)
     797    TwinLaw = np.ones(1)
     798#    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
     799#    TwDict = refDict.get('TwDict',{})           
     800#    if 'S' in calcControls[hfx+'histType']:
     801#        NTL = calcControls[phfx+'NTL']
     802#        NM = calcControls[phfx+'TwinNMN']+1
     803#        TwinLaw = calcControls[phfx+'TwinLaw']
     804#        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
     805#        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
     806    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     807        GetAtomFXU(pfx,calcControls,parmDict)
     808    if not Xdata.size:          #no atoms in phase!
     809        return
     810    Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
     811    Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
     812    Gdata = np.inner(Bmat,Gdata.T)              #convert to crystal space
     813    Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
     814    if SGData['SGInv'] and not SGData['SGFixed']:
     815        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
     816    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
     817    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
     818    Gdata = np.inner(Amat,Gdata.T)              #convert back to cart. space MXYZ, Natoms, NOps*Inv*Ncen
     819    Gdata = np.swapaxes(Gdata,1,2)              # put Natoms last
     820    Mag = np.tile(Mag[:,nxs],len(SGMT)*Ncen).T
     821    if SGData['SGInv'] and not SGData['SGFixed']:
     822        Mag = np.repeat(Mag,2,axis=0)                  #Mag same shape as Gdata
     823    Uij = np.array(G2lat.U6toUij(Uijdata))
     824    bij = Mast*Uij.T
     825    blkSize = 100       #no. of reflections in a block - size seems optimal
     826    nRef = refDict['RefList'].shape[0]
     827    SQ = 1./(2.*refDict['RefList'].T[4])**2
     828    refDict['FF']['El'] = list(MFtables.keys())
     829    refDict['FF']['MF'] = np.zeros((nRef,len(MFtables)))
     830    for iel,El in enumerate(refDict['FF']['El']):
     831        refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
     832#reflection processing begins here - big arrays!
     833    iBeg = 0
     834    while iBeg < nRef:
     835        iFin = min(iBeg+blkSize,nRef)
     836        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
     837        H = refl.T[:3]                          #array(blkSize,3)
     838#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
     839#        TwMask = np.any(H,axis=-1)
     840#        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
     841#            for ir in range(blkSize):
     842#                iref = ir+iBeg
     843#                if iref in TwDict:
     844#                    for i in TwDict[iref]:
     845#                        for n in range(NTL):
     846#                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     847#            TwMask = np.any(H,axis=-1)
     848        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
     849        SQfactor = 4.0*SQ*twopisq               #ditto prev.
     850        Uniq = np.inner(H.T,SGMT)
     851        Phi = np.inner(H.T,SGT)
     852        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
     853        biso = -SQfactor*Uisodata[:,nxs]
     854        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
     855        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
     856        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
     857        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     858        MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
     859        TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
     860        if SGData['SGInv'] and not SGData['SGFixed']:
     861            mphase = np.hstack((phase,-phase))
     862        else:
     863            mphase = phase
     864        mphase = np.array([mphase+twopi*np.inner(cen,H.T)[:,nxs,nxs] for cen in SGData['SGCen']])
     865        mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
     866        sinm = np.sin(mphase)                               #ditto - match magstrfc.for
     867        cosm = np.cos(mphase)                               #ditto
     868        HM = np.inner(Bmat.T,H.T)                             #put into cartesian space
     869        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
     870        eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
     871        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
     872        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
     873        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
     874        fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
     875        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
     876        refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
     877        refl.T[7] = np.copy(refl.T[9])               
     878        refl.T[10] = 0.0    #atan2d(fbs[0],fas[0]) - what is phase for mag refl?
     879#        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
     880#            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
     881#            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     882#        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
     883#            if len(TwinLaw) > 1:
     884#                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
     885#                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
     886#                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
     887#                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
     888#            else:   # checked correct!!
     889#                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
     890#                refl.T[7] = np.copy(refl.T[9])               
     891#                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     892##                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
     893        iBeg += blkSize
     894#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
     895
    815896def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    816897    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
     
    31833264            if im:
    31843265                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
     3266            elif parmDict[pfx+'isMag'] and 'N' in calcControls[hfx+'histType']:
     3267                MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    31853268            else:
    31863269                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
Note: See TracChangeset for help on using the changeset viewer.