Changeset 4519


Ignore:
Timestamp:
Jul 13, 2020 1:27:06 PM (3 years ago)
Author:
vondreele
Message:

G2fpaGUI - minor cleanup
G2lattice - test on 'T' instead of 'C' so that 'B' PWDR types properly handled
G2math - add getPinkalpha, getPinkbeta % deriv routines for pink beam peak shapes

modify setPeakparams to use them

G2plot - test on 'T' instead of 'C' so that 'B' PWDR types properly handled
G2pwd - add pink beam peak shape function - same as TOF peak shape; scale sig & gam to be in centidegrees

  • add wtFactor to FitPeaks? fix bad PWDR weights for peak fitting so reduced chi2 is more appropriate

G2pwdGUI - test on 'T' instead of 'C' so that 'B' PWDR types properly handled

  • add wtFactor to DoPeakFit? to fix bad PWDR weights for peak fitting so reduced chi2 is more appropriate

G2strMath - latest in incommensurate mag str factors - no real improvement

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIfpaGUI.py

    r4403 r4519  
    587587        except:
    588588            pass
    589         GoOn = pgbar.Update(5,newmsg='Creating peak list')
     589        pgbar.Update(5,newmsg='Creating peak list')
    590590        pgbar.Raise()
    591591        for pos in peaklist:
    … …  
    593593            area = sum(intArr[max(0,i-maxPtsHM):min(len(intArr),i+maxPtsHM)])
    594594            peakData['peaks'].append(G2mth.setPeakparms(Parms,Parms2,pos,area))
    595         GoOn = pgbar.Update(10,newmsg='Refining peak positions')
     595        pgbar.Update(10,newmsg='Refining peak positions')
    596596        histData = G2frame.GPXtree.GetItemPyData(histId)
    597597        # refine peak positions only
    … …  
    600600                                            bkg,limits[1],
    601601                                            Parms,Parms2,histData[1],bxye,[],
    602                                            False,controldat,None)[0]
    603         GoOn = pgbar.Update(20,newmsg='Refining peak positions && areas')
     602                                           False,controldat)[0]
     603        pgbar.Update(20,newmsg='Refining peak positions && areas')
    604604        # refine peak areas as well
    605605        for pk in peakData['peaks']:
    … …  
    609609                                            Parms,Parms2,histData[1],bxye,[],
    610610                                           False,controldat)[0]
    611         GoOn = pgbar.Update(40,newmsg='Refining profile function')
     611        pgbar.Update(40,newmsg='Refining profile function')
    612612        # refine profile function
    613613        for p in ('U', 'V', 'W', 'X', 'Y'):
    … …  
    617617                                            Parms,Parms2,histData[1],bxye,[],
    618618                                           False,controldat)[0]
    619         GoOn = pgbar.Update(70,newmsg='Refining profile function && asymmetry')
     619        pgbar.Update(70,newmsg='Refining profile function && asymmetry')
    620620        # add in asymmetry
    621621        Parms['SH/L'][2] = True
    … …  
    624624                                            Parms,Parms2,histData[1],bxye,[],
    625625                                           False,controldat)[0]
    626         GoOn = pgbar.Update(100,newmsg='Done')
     626        pgbar.Update(100,newmsg='Done')
    627627        # reset "initial" profile
    628628        for p in Parms:
    … …  
    688688            txt = open(filename[0],'r').read()
    689689            NISTparms.clear()
    690             array = np.array
    691690            d = eval(txt)
    692691            NISTparms.update(d)

  • trunk/GSASIIlattice.py

    r4440 r4519  
    842842    ''' convert powder pattern position (2-theta or TOF, musec) to d-spacing
    843843    '''
    844     if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
     844    if 'T' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
     845        return TOF2dsp(Inst,pos)
     846    else:   #'C' or 'B'
    845847        wave = G2mth.getWave(Inst)
    846848        return wave/(2.0*sind((pos-Inst.get('Zero',[0,0])[1])/2.0))
    847     else:   #'T'OF - ignore difB
    848         return TOF2dsp(Inst,pos)
    849849       
    850850def TOF2dsp(Inst,Pos):
    … …  
    868868    ''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
    869869    '''
    870     if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
     870    if 'T' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
     871        pos = Inst['difC'][1]*dsp+Inst['Zero'][1]+Inst['difA'][1]*dsp**2+Inst.get('difB',[0,0,False])[1]/dsp
     872    else:   #'C' or 'B'
    871873        wave = G2mth.getWave(Inst)
    872874        val = min(0.995,wave/(2.*dsp))  #set max at 168deg
    873875        pos = 2.0*asind(val)+Inst.get('Zero',[0,0])[1]             
    874     else:   #'T'OF
    875         pos = Inst['difC'][1]*dsp+Inst['Zero'][1]+Inst['difA'][1]*dsp**2+Inst.get('difB',[0,0,False])[1]/dsp
    876876    return pos
    877877   
    … …  
    879879    ''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
    880880    '''
    881     if 'C' in dataType:
     881    if 'T' in dataType:
     882        pos = parmdict['difC']*dsp+parmdict['difA']*dsp**2+parmdict['difB']/dsp+parmdict['Zero']
     883    else:   #'C' or 'B'
    882884        pos = 2.0*asind(parmdict['Lam']/(2.*dsp))+parmdict['Zero']
    883     else:   #'T'OF
    884         pos = parmdict['difC']*dsp+parmdict['difA']*dsp**2+parmdict['difB']/dsp+parmdict['Zero']
    885885    return pos
    886886                   

  • trunk/GSASIImath.py

    r4514 r4519  
    14631463    psin = np.sin(twopi*phase)
    14641464    pcos = np.cos(twopi*phase)
    1465     MmodA = np.sum(Am[nxs,nxs,:,:,:]*pcos[:,:,:,nxs,nxs],axis=3)    #cos term
    1466     MmodB = np.sum(Bm[nxs,nxs,:,:,:]*psin[:,:,:,nxs,nxs],axis=3)    #sin term
     1465    MmodA = np.sum(Am[nxs,nxs,:,:,:]*pcos[:,:,:,nxs,nxs],axis=3)/2.    #cos term
     1466    MmodB = np.sum(Bm[nxs,nxs,:,:,:]*psin[:,:,:,nxs,nxs],axis=3)/2.    #sin term
    14671467    MmodA = np.sum(SGMT[nxs,:,nxs,:,:]*MmodA[:,:,:,nxs,:],axis=-1)*SGData['SpnFlp'][nxs,:,nxs,nxs]
    14681468    MmodB = np.sum(SGMT[nxs,:,nxs,:,:]*MmodB[:,:,:,nxs,:],axis=-1)*SGData['SpnFlp'][nxs,:,nxs,nxs]
    14691469    return MmodA,MmodB    #Ntau,Nops,Natm,Mxyz; cos & sin parts; sum matches drawn atom moments
    14701470       
    1471 def MagMod2(m,glTau,XYZ,modQ,MSSdata,SGData,SSGData):
     1471def MagMod2(m,XYZ,modQ,MSSdata,SGData,SSGData):
    14721472    '''
    14731473    this needs to make magnetic moment modulations & magnitudes as
    … …  
    14771477    Bm = np.array(MSSdata[:3]).T[:,0,:]  #...sin pos mods
    14781478    SGMT = np.array([ops[0] for ops in SGData['SGOps']])        #not .T!!
     1479    SSGMT = np.array([ops[0] for ops in SSGData['SSGOps']])        #not .T!!
    14791480    Sinv = np.array([nl.inv(ops[0]) for ops in SSGData['SSGOps']])
    14801481    SGT = np.array([ops[1] for ops in SSGData['SSGOps']])
    14811482    if SGData['SGInv']:
    14821483        SGMT = np.vstack((SGMT,-SGMT))
     1484        SSGMT = np.vstack((SSGMT,-SSGMT))
    14831485        Sinv = np.vstack((Sinv,-Sinv))
    14841486        SGT = np.vstack((SGT,-SGT))
    14851487    SGMT = np.vstack([SGMT for cen in SGData['SGCen']])
     1488    SSGMT = np.vstack([SSGMT for cen in SGData['SGCen']])
    14861489    Sinv = np.vstack([Sinv for cen in SGData['SGCen']])
    14871490    SGT = np.vstack([SGT+cen for cen in SSGData['SSGCen']])%1.
    14881491    if SGData['SGGray']:
    14891492        SGMT = np.vstack((SGMT,SGMT))
     1493        SSGMT = np.vstack((SSGMT,SSGMT))
    14901494        Sinv = np.vstack((Sinv,Sinv))
    14911495        SGT = np.vstack((SGT,SGT+.5))%1.
    1492     mst = Sinv[:,3,:3]
    14931496    epsinv = Sinv[:,3,3]
    1494     phase = np.inner(XYZ,modQ).T+(np.inner(mst,modQ)-epsinv)[:,nxs]+glTau
    1495    
    1496     psin = np.sin(twopi*m*phase).T
    1497     pcos = np.cos(twopi*m*phase).T
    1498     MmodA = Am[nxs,nxs,:,:]*pcos[:,:,nxs,nxs]    #cos term
    1499     MmodB = Bm[nxs,nxs,:,:]*psin[:,:,nxs,nxs]    #sin term
    1500     MmodA = np.sum(SGMT[nxs,:,nxs,:,:]*MmodA[:,:,:,nxs,:],axis=-1)*SGData['SpnFlp'][nxs,:,nxs,nxs]
    1501     MmodB = np.sum(SGMT[nxs,:,nxs,:,:]*MmodB[:,:,:,nxs,:],axis=-1)*SGData['SpnFlp'][nxs,:,nxs,nxs]
    1502     return MmodA,MmodB    #Nref,Ntau,Nops,Natm,Mxyz; cos & sin parts; sum matches drawn atom moments
     1497    phi = np.inner(XYZ,modQ).T
     1498    TA = phi+(epsinv*(np.inner(modQ,SGT[:,:3])-SGT[:,3]))[:,nxs]    #Nops,Natm
     1499    phase = phi+(np.inner(modQ,SGT[:,:3])-SGT[:,3])[:,nxs]
     1500   
     1501    pcos = np.cos(-twopi*m[:,nxs,nxs]*phase[nxs,:,:])      #Nref,Nops,Natm
     1502    psin = np.sin(-twopi*m[:,nxs,nxs]*phase[nxs,:,:])
     1503    MmodA = TA[nxs,:,:,nxs]*(Am[nxs,nxs,:,:]*pcos[:,:,:,nxs]-Bm[nxs,nxs,:,:]*psin[:,:,:,nxs])/2.    #Nref,Nops,Natm,Mxyz
     1504    MmodB = TA[nxs,:,:,nxs]*(Am[nxs,nxs,:,:]*psin[:,:,:,nxs]+Bm[nxs,nxs,:,:]*pcos[:,:,:,nxs])/2.    #Nref,Nops,Natm,Mxyz
     1505    MmodA = np.sum(SGMT[nxs,:,nxs,:,:]*MmodA[:,:,:,nxs,:],axis=-1)*SGData['MagMom'][nxs,:,nxs,nxs]
     1506    MmodB = np.sum(SGMT[nxs,:,nxs,:,:]*MmodB[:,:,:,nxs,:],axis=-1)*SGData['MagMom'][nxs,:,nxs,nxs]
     1507    return MmodA,MmodB    #Nref,Nops,Natm,Mxyz; cos & sin parts
    15031508       
    15041509def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
    … …  
    41634168    return 1./dsp
    41644169   
     4170def getPinkalpha(ins,tth):
     4171    '''get TOF peak profile alpha
     4172   
     4173    :param dict ins: instrument parameters with at least 'alpha'
     4174      as values only
     4175    :param float tth: 2-theta of peak
     4176   
     4177    :returns: flaot getPinkalpha: peak alpha
     4178   
     4179    '''
     4180    return ins['alpha-0']+ ins['alpha-1']*tand(tth/2.)
     4181   
     4182def getPinkalphaDeriv(tth):
     4183    '''get derivatives of TOF peak profile beta wrt alpha
     4184   
     4185    :param float dsp: d-spacing of peak
     4186   
     4187    :returns: float getTOFalphaDeriv: d(alp)/d(alpha-0), d(alp)/d(alpha-1)
     4188   
     4189    '''
     4190    return 1.0,tand(tth/2.)
     4191   
     4192def getPinkbeta(ins,tth):
     4193    '''get TOF peak profile beta
     4194   
     4195    :param dict ins: instrument parameters with at least 'beat-0' & 'beta-1'
     4196      as values only
     4197    :param float tth: 2-theta of peak
     4198   
     4199    :returns: float getaPinkbeta: peak beta
     4200   
     4201    '''
     4202    return ins['beta-0']+ins['beta-1']*tand(tth/2.)
     4203   
     4204def getPinkbetaDeriv(tth):
     4205    '''get derivatives of TOF peak profile beta wrt beta-0 & beta-1
     4206   
     4207    :param float dsp: d-spacing of peak
     4208   
     4209    :returns: list getTOFbetaDeriv: d(beta)/d(beta-0) & d(beta)/d(beta-1)
     4210   
     4211    '''
     4212    return 1.0,tand(tth/2.)
     4213   
    41654214def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
    41664215    '''set starting peak parameters for single peak fits from plot selection or auto selection
    … …  
    41764225        for CW: [pos,0,mag,1,sig,0,gam,0]
    41774226        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
     4227        for Pink: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
    41784228        NB: mag refinement set by default, all others off
    41794229   
    … …  
    41834233        ind = 1
    41844234    ins = {}
    4185     if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an else
    4186         for x in ['U','V','W','X','Y','Z']:
    4187             ins[x] = Parms.get(x,[0.0,0.0])[ind]
    4188         if ifQ:                              #qplot - convert back to 2-theta
    4189             pos = 2.0*asind(pos*getWave(Parms)/(4*math.pi))
    4190         sig = getCWsig(ins,pos)
    4191         gam = getCWgam(ins,pos)           
    4192         XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
    4193     else:
     4235    if 'T' in Parms['Type'][0]:
    41944236        if ifQ:
    41954237            dsp = 2.*np.pi/pos
    … …  
    42114253        gam = getTOFgamma(ins,dsp)
    42124254        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
     4255    elif 'C' in Parms['Type'][0]:                            #CW data - TOF later in an else
     4256        for x in ['U','V','W','X','Y','Z']:
     4257            ins[x] = Parms.get(x,[0.0,0.0])[ind]
     4258        if ifQ:                              #qplot - convert back to 2-theta
     4259            pos = 2.0*asind(pos*getWave(Parms)/(4*math.pi))
     4260        sig = getCWsig(ins,pos)
     4261        gam = getCWgam(ins,pos)           
     4262        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
     4263    elif 'B' in Parms['Type'][0]:
     4264        for x in ['U','V','W','X','Y','Z','alpha-0','alpha-1','beta-0','beta-1']:
     4265            ins[x] = Parms.get(x,[0.0,0.0])[ind]
     4266        if ifQ:                              #qplot - convert back to 2-theta
     4267            pos = 2.0*asind(pos*getWave(Parms)/(4*math.pi))
     4268        alp = getPinkalpha(ins,pos)
     4269        bet = getPinkbeta(ins,pos)
     4270        sig = getCWsig(ins,pos)
     4271        gam = getCWgam(ins,pos)           
     4272        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]       #default refine intensity 1st
    42134273    return XY
    42144274   

  • trunk/GSASIIplot.py

    r4516 r4519  
    19451945                    G2frame.G2plotNB.status.SetStatusText('q = %9.5f'%q)
    19461946                    return
    1947                 if 'C' in Parms['Type'][0] or 'PKS' in Parms['Type'][0]:
     1947                if 'T' in Parms['Type'][0]: # TOF
     1948                    dT = Parms['difC'][1] * 2 * np.pi * tolerance / q**2
     1949                else: # 'C' or  'B' in Parms['Type'][0] or 'PKS' in Parms['Type'][0]:
    19481950                    wave = G2mth.getWave(Parms)
    19491951                    dT = tolerance*wave*90./(np.pi**2*cosd(xpos/2))
    1950                 else: # TOF
    1951                     dT = Parms['difC'][1] * 2 * np.pi * tolerance / q**2
    19521952            elif plottype in ['SASD','REFD']:
    19531953                q = xpos
    … …  
    19761976                q = 2.*np.pi/dsp
    19771977            if G2frame.Contour: #PWDR only
    1978                 if 'C' in Parms['Type'][0]:
     1978                if 'T' in Parms['Type'][0]:
     1979                    G2frame.G2plotNB.status.SetStatusText('TOF =%9.3f d =%9.5f q = %9.5f pattern ID =%5d'%(xpos,dsp,q,int(ypos)),1)
     1980                else:
    19791981                    G2frame.G2plotNB.status.SetStatusText('2-theta =%9.3f d =%9.5f q = %9.5f pattern ID =%5d'%(xpos,dsp,q,int(ypos)),1)
     1982            else:
     1983                if 'T' in Parms['Type'][0]:
     1984                    if Page.plotStyle['sqrtPlot']:
     1985                        G2frame.G2plotNB.status.SetStatusText('TOF =%9.3f d =%9.5f q =%9.5f sqrt(Intensity) =%9.2f'%(xpos,dsp,q,ypos),1)
     1986                    else:
     1987                        G2frame.G2plotNB.status.SetStatusText('TOF =%9.3f d =%9.5f q =%9.5f Intensity =%9.2f'%(xpos,dsp,q,ypos),1)
    19801988                else:
    1981                     G2frame.G2plotNB.status.SetStatusText('TOF =%9.3f d =%9.5f q = %9.5f pattern ID =%5d'%(xpos,dsp,q,int(ypos)),1)
    1982             else:
    1983                 if 'C' in Parms['Type'][0]:
    19841989                    if 'PWDR' in plottype:
    19851990                        if Page.plotStyle['sqrtPlot']:
    … …  
    19911996                    elif plottype == 'REFD':
    19921997                        G2frame.G2plotNB.status.SetStatusText('q =%12.5g Reflectivity =%12.5g d =%9.1f'%(q,ypos,dsp),1)
    1993                 else:
    1994                     if Page.plotStyle['sqrtPlot']:
    1995                         G2frame.G2plotNB.status.SetStatusText('TOF =%9.3f d =%9.5f q =%9.5f sqrt(Intensity) =%9.2f'%(xpos,dsp,q,ypos),1)
    1996                     else:
    1997                         G2frame.G2plotNB.status.SetStatusText('TOF =%9.3f d =%9.5f q =%9.5f Intensity =%9.2f'%(xpos,dsp,q,ypos),1)
    19981998            s = ''
    19991999            if G2frame.PickId:
    … …  
    21602160            if ind.all() != [0] and ObsLine[0].get_label() in str(pick):    #picked a data point, add a new peak
    21612161                data = G2frame.GPXtree.GetItemPyData(G2frame.PickId)
    2162                 XY = G2mth.setPeakparms(Parms,Parms2,xy[0],xy[1])
     2162                XY = G2mth.setPeakparms(Parms,Parms2,xy[0],xy[1],useFit=True)
    21632163                data['peaks'].append(XY)
    21642164                data['sigDict'] = {}    #now invalid
    … …  
    28912891    else:
    28922892        Plot.set_xlabel(xLabel,fontsize=16)
    2893     if 'C' in ParmList[0]['Type'][0]:
     2893    if 'T' in ParmList[0]['Type'][0]:
     2894        if Page.plotStyle['sqrtPlot']:
     2895            Plot.set_ylabel(r'$\sqrt{Normalized\ intensity}$',fontsize=16)
     2896        else:
     2897            Plot.set_ylabel(r'$Normalized\ intensity$',fontsize=16)
     2898    else:       #neutron TOF
    28942899        if 'PWDR' in plottype:
    28952900            if Page.plotStyle['sqrtPlot']:
    … …  
    29092914            else:
    29102915                Plot.set_ylabel(r'$Reflectivity$',fontsize=16)               
    2911     else:       #neutron TOF
    2912         if Page.plotStyle['sqrtPlot']:
    2913             Plot.set_ylabel(r'$\sqrt{Normalized\ intensity}$',fontsize=16)
    2914         else:
    2915             Plot.set_ylabel(r'$Normalized\ intensity$',fontsize=16)
    29162916    mpl.rcParams['image.cmap'] = G2frame.ContourColor
    29172917    mcolors = mpl.cm.ScalarMappable()       #wants only default as defined in previous line!!
    … …  
    29592959                magLineList = data[0].get('Magnification',[])
    29602960            if ('C' in ParmList[0]['Type'][0] and Page.plotStyle['dPlot']) or \
     2961                ('B' in ParmList[0]['Type'][0] and Page.plotStyle['dPlot']) or \
    29612962                ('T' in ParmList[0]['Type'][0] and Page.plotStyle['qPlot']): # reversed regions relative to data order
    29622963                tcorner = 1
    … …  
    55595560    Z = []
    55605561    W = []
    5561     if 'C' in Parms['Type'][0]:
     5562    if 'T' in Parms['Type'][0]:   #'T'OF
     5563        Plot.set_title('Instrument and sample peak coefficients')
     5564        Plot.set_xlabel(r'$Q, \AA^{-1}$',fontsize=14)
     5565        Plot.set_ylabel(r'$\alpha, \beta, \Delta Q/Q, \Delta d/d$',fontsize=14)
     5566        Xmin,Xmax = limits[1]
     5567        T = np.linspace(Xmin,Xmax,num=101,endpoint=True)
     5568        Z = np.ones_like(T)
     5569        data = G2mth.setPeakparms(Parms,Parms2,T,Z)
     5570        ds = T/difC
     5571        Q = 2.*np.pi/ds
     5572        A = data[4]
     5573        B = data[6]
     5574        S = 1.17741*np.sqrt(data[8])/T
     5575        G = data[10]/T
     5576        Plot.plot(Q,A,color='r',label='Alpha')
     5577        Plot.plot(Q,B,color='g',label='Beta')
     5578        Plot.plot(Q,S,color='b',label='Gaussian')
     5579        Plot.plot(Q,G,color='m',label='Lorentzian')
     5580
     5581        fit = G2mth.setPeakparms(Parms,Parms2,T,Z)
     5582        ds = T/difC
     5583        Q = 2.*np.pi/ds
     5584        Af = fit[4]
     5585        Bf = fit[6]
     5586        Sf = 1.17741*np.sqrt(fit[8])/T
     5587        Gf = fit[10]/T
     5588        Plot.plot(Q,Af,color='r',dashes=(5,5),label='Alpha fit')
     5589        Plot.plot(Q,Bf,color='g',dashes=(5,5),label='Beta fit')
     5590        Plot.plot(Q,Sf,color='b',dashes=(5,5),label='Gaussian fit')
     5591        Plot.plot(Q,Gf,color='m',dashes=(5,5),label='Lorentzian fit')
     5592       
     5593        Tp = []
     5594        Ap = []
     5595        Bp = []
     5596        Sp = []
     5597        Gp = []
     5598        Wp = []
     5599        Qp = []
     5600        for peak in peaks:
     5601            Tp.append(peak[0])
     5602            Ap.append(peak[4])
     5603            Bp.append(peak[6])
     5604            Qp.append(2.*np.pi*difC/peak[0])
     5605            Sp.append(1.17741*np.sqrt(peak[8])/peak[0])
     5606            Gp.append(peak[10]/peak[0])
     5607           
     5608        if Qp:
     5609            Plot.plot(Qp,Ap,'+',color='r',label='Alpha peak')
     5610            Plot.plot(Qp,Bp,'+',color='g',label='Beta peak')
     5611            Plot.plot(Qp,Sp,'+',color='b',label='Gaussian peak')
     5612            Plot.plot(Qp,Gp,'+',color='m',label='Lorentzian peak')
     5613        Plot.legend(loc='best')
     5614    else:       #'C' & 'B'
    55625615        Plot.figure.suptitle(TreeItemText)
    55635616        Plot.set_title('Instrument and sample peak widths')
    … …  
    56125665        SetupLegendPick(legend,new)
    56135666        Page.canvas.draw()
    5614     else:   #'T'OF
    5615         Plot.set_title('Instrument and sample peak coefficients')
    5616         Plot.set_xlabel(r'$Q, \AA^{-1}$',fontsize=14)
    5617         Plot.set_ylabel(r'$\alpha, \beta, \Delta Q/Q, \Delta d/d$',fontsize=14)
    5618         Xmin,Xmax = limits[1]
    5619         T = np.linspace(Xmin,Xmax,num=101,endpoint=True)
    5620         Z = np.ones_like(T)
    5621         data = G2mth.setPeakparms(Parms,Parms2,T,Z)
    5622         ds = T/difC
    5623         Q = 2.*np.pi/ds
    5624         A = data[4]
    5625         B = data[6]
    5626         S = 1.17741*np.sqrt(data[8])/T
    5627         G = data[10]/T
    5628         Plot.plot(Q,A,color='r',label='Alpha')
    5629         Plot.plot(Q,B,color='g',label='Beta')
    5630         Plot.plot(Q,S,color='b',label='Gaussian')
    5631         Plot.plot(Q,G,color='m',label='Lorentzian')
    5632 
    5633         fit = G2mth.setPeakparms(Parms,Parms2,T,Z)
    5634         ds = T/difC
    5635         Q = 2.*np.pi/ds
    5636         Af = fit[4]
    5637         Bf = fit[6]
    5638         Sf = 1.17741*np.sqrt(fit[8])/T
    5639         Gf = fit[10]/T
    5640         Plot.plot(Q,Af,color='r',dashes=(5,5),label='Alpha fit')
    5641         Plot.plot(Q,Bf,color='g',dashes=(5,5),label='Beta fit')
    5642         Plot.plot(Q,Sf,color='b',dashes=(5,5),label='Gaussian fit')
    5643         Plot.plot(Q,Gf,color='m',dashes=(5,5),label='Lorentzian fit')
    5644        
    5645         Tp = []
    5646         Ap = []
    5647         Bp = []
    5648         Sp = []
    5649         Gp = []
    5650         Wp = []
    5651         Qp = []
    5652         for peak in peaks:
    5653             Tp.append(peak[0])
    5654             Ap.append(peak[4])
    5655             Bp.append(peak[6])
    5656             Qp.append(2.*np.pi*difC/peak[0])
    5657             Sp.append(1.17741*np.sqrt(peak[8])/peak[0])
    5658             Gp.append(peak[10]/peak[0])
    5659            
    5660         if Qp:
    5661             Plot.plot(Qp,Ap,'+',color='r',label='Alpha peak')
    5662             Plot.plot(Qp,Bp,'+',color='g',label='Beta peak')
    5663             Plot.plot(Qp,Sp,'+',color='b',label='Gaussian peak')
    5664             Plot.plot(Qp,Gp,'+',color='m',label='Lorentzian peak')
    5665         Plot.legend(loc='best')
     5667       
    56665668    if xylim and not G2frame.G2plotNB.allowZoomReset:
    56675669        # this restores previous plot limits (but I'm not sure why there are two .push_current calls)

  • trunk/GSASIIpwd.py

    r4518 r4519  
    801801    alpTOF = lambda dsp,alp: alp/dsp
    802802    betTOF = lambda dsp,bet0,bet1,betq: bet0+bet1/dsp**4+betq/dsp**2
    803     if 'C' in Inst['Type'][0]:
    804         s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
    805         g = gam(pos/2.,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
    806         return getgamFW(g,s)/100.  #returns FWHM in deg
    807     else:
     803    if 'T' in Inst['Type'][0]:
    808804        dsp = pos/Inst['difC'][0]
    809805        alp = alpTOF(dsp,Inst['alpha'][0])
    … …  
    812808        g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
    813809        return getgamFW(g,s)+np.log(2.0)*(alp+bet)/(alp*bet)
     810    else:
     811        s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
     812        g = gam(pos/2.,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
     813        return getgamFW(g,s)/100.  #returns FWHM in deg
    814814   
    815815def getgamFW(g,s):
    … …  
    854854    if 'T' in dataType:
    855855        q = 2.*np.pi*parmDict[pfx+'difC']/xdata
    856     elif 'C' in dataType:
     856    else:
    857857        wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
    858858        q = npT2q(xdata,wave)
    … …  
    990990    if 'T' in dataType:
    991991        q = 2.*np.pi*parmDict[hfx+'difC']/xdata
    992     elif 'C' in dataType:
     992    else:
    993993        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
    994994        q = 2.*np.pi*npsind(xdata/2.)/wave
    … …  
    10871087            if 'C' in dataType:
    10881088                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
    1089             else: #'T'OF
     1089            else: #'T' or 'B'
    10901090                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
    10911091            iBeg = np.searchsorted(xdata,pkP-fmin)
    … …  
    13211321            except KeyError:        #no more peaks to process
    13221322                return yb+yc
     1323    elif 'B' in dataType:
     1324        iPeak = 0
     1325        dsp = 1.0 #for now - fix later
     1326        while True:
     1327            try:
     1328                pos = parmDict['pos'+str(iPeak)]
     1329                tth = (pos-parmDict['Zero'])
     1330                intens = parmDict['int'+str(iPeak)]
     1331                alpName = 'alp'+str(iPeak)
     1332                if alpName in varyList:
     1333                    alp = parmDict[alpName]
     1334                else:
     1335                    alp = G2mth.getPinkalpha(parmDict,tth)
     1336                alp = max(0.1,alp)
     1337                betName = 'bet'+str(iPeak)
     1338                if betName in varyList:
     1339                    bet = parmDict[betName]
     1340                else:
     1341                    bet = G2mth.getPinkbeta(parmDict,tth)
     1342                bet = max(0.0001,bet)
     1343                sigName = 'sig'+str(iPeak)
     1344                if sigName in varyList:
     1345                    sig = parmDict[sigName]
     1346                else:
     1347                    sig = G2mth.getCWsig(parmDict,tth)
     1348                sig = max(sig,0.001)          #avoid neg sigma^2
     1349                gamName = 'gam'+str(iPeak)
     1350                if gamName in varyList:
     1351                    gam = parmDict[gamName]
     1352                else:
     1353                    gam = G2mth.getCWgam(parmDict,tth)
     1354                gam = max(gam,0.001)             #avoid neg gamma
     1355                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
     1356                iBeg = np.searchsorted(xdata,pos-fmin)
     1357                iFin = np.searchsorted(xdata,pos+fmin)
     1358                if not iBeg+iFin:       #peak below low limit
     1359                    iPeak += 1
     1360                    continue
     1361                elif not iBeg-iFin:     #peak above high limit
     1362                    return yb+yc
     1363                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])
     1364                iPeak += 1
     1365            except KeyError:        #no more peaks to process
     1366                return yb+yc       
    13231367    else:
    13241368        Pdabc = parmDict['Pdabc']
    … …  
    14801524            except KeyError:        #no more peaks to process
    14811525                break
     1526    elif 'B' in dataType:
     1527        iPeak = 0
     1528        while True:
     1529            try:
     1530                pos = parmDict['pos'+str(iPeak)]
     1531                tth = (pos-parmDict['Zero'])
     1532                intens = parmDict['int'+str(iPeak)]
     1533                alpName = 'alp'+str(iPeak)
     1534                if alpName in varyList:
     1535                    alp = parmDict[alpName]
     1536                else:
     1537                    alp = G2mth.getPinkalpha(parmDict,tth)
     1538                    dada0,dada1 = G2mth.getPinkalphaDeriv(tth)
     1539                alp = max(0.0001,alp)
     1540                betName = 'bet'+str(iPeak)
     1541                if betName in varyList:
     1542                    bet = parmDict[betName]
     1543                else:
     1544                    bet = G2mth.getPinkbeta(parmDict,tth)
     1545                    dbdb0,dbdb1 = G2mth.getPinkbetaDeriv(tth)
     1546                bet = max(0.0001,bet)
     1547                sigName = 'sig'+str(iPeak)
     1548                if sigName in varyList:
     1549                    sig = parmDict[sigName]
     1550                    dsdU = dsdV = dsdW = 0
     1551                else:
     1552                    sig = G2mth.getCWsig(parmDict,tth)
     1553                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
     1554                sig = max(sig,0.001)          #avoid neg sigma
     1555                gamName = 'gam'+str(iPeak)
     1556                if gamName in varyList:
     1557                    gam = parmDict[gamName]
     1558                    dgdX = dgdY = dgdZ = 0
     1559                else:
     1560                    gam = G2mth.getCWgam(parmDict,tth)
     1561                    dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
     1562                gam = max(gam,0.001)             #avoid neg gamma
     1563                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig/1.e4,gam/100.)
     1564                iBeg = np.searchsorted(xdata,pos-fmin)
     1565                iFin = np.searchsorted(xdata,pos+fmin)
     1566                if not iBeg+iFin:       #peak below low limit
     1567                    iPeak += 1
     1568                    continue
     1569                elif not iBeg-iFin:     #peak above high limit
     1570                    break
     1571                dMdpk = np.zeros(shape=(7,len(xdata)))
     1572                dMdipk = getdEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])
     1573                for i in range(1,6):
     1574                    dMdpk[i][iBeg:iFin] += cw[iBeg:iFin]*intens*dMdipk[i]
     1575                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
     1576                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4]/1.e4,'gam':dMdpk[5]/100.}
     1577                for parmName in ['pos','int','alp','bet','sig','gam']:
     1578                    try:
     1579                        idx = varyList.index(parmName+str(iPeak))
     1580                        dMdv[idx] = dervDict[parmName]
     1581                    except ValueError:
     1582                        pass
     1583                if 'U' in varyList:
     1584                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
     1585                if 'V' in varyList:
     1586                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
     1587                if 'W' in varyList:
     1588                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
     1589                if 'X' in varyList:
     1590                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
     1591                if 'Y' in varyList:
     1592                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
     1593                if 'Z' in varyList:
     1594                    dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
     1595                if 'alpha-0' in varyList:
     1596                    dMdv[varyList.index('alpha-0')] += dada0*dervDict['alp']
     1597                if 'alpha-1' in varyList:
     1598                    dMdv[varyList.index('alpha-1')] += dada1*dervDict['alp']
     1599                if 'beta-0' in varyList:
     1600                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
     1601                if 'beta-1' in varyList:
     1602                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
     1603                iPeak += 1
     1604            except KeyError:        #no more peaks to process
     1605                break       
    14821606    else:
    14831607        Pdabc = parmDict['Pdabc']
    … …  
    17201844    return True
    17211845           
    1722 def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
     1846def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[],oneCycle=False,controls=None,wtFactor=1.0,dlg=None):
    17231847    '''Called to perform a peak fit, refining the selected items in the peak
    17241848    table as well as selected items in the background.
    … …  
    17441868    :param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
    17451869      and derivType = controls['deriv type']. If None default values are used.
     1870    :param float wtFactor. weight multiplier; = 1.0 by default
    17461871    :param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
    17471872      Defaults to None, which means no updates are done.
    … …  
    18421967            insVals.append(Inst[parm][1])
    18431968            if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
    1844                 'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
     1969                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1'] and Inst[parm][2]:
    18451970                    insVary.append(parm)
    18461971        instDict = dict(zip(insNames,insVals))
    … …  
    18601985                pos = parmDict['pos'+str(iPeak)]
    18611986                if sigName not in varyList:
    1862                     if 'C' in Inst['Type'][0]:
    1863                         parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
    1864                     else:
     1987                    if 'T' in Inst['Type'][0]:
    18651988                        dsp = G2lat.Pos2dsp(Inst,pos)
    18661989                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
     1990                    else:
     1991                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
    18671992                gamName = 'gam'+str(iPeak)
    18681993                if gamName not in varyList:
    1869                     if 'C' in Inst['Type'][0]:
    1870                         parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
    1871                     else:
     1994                    if 'T' in Inst['Type'][0]:
    18721995                        dsp = G2lat.Pos2dsp(Inst,pos)
    18731996                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
     1997                    else:
     1998                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
    18741999                iPeak += 1
    18752000            except KeyError:
    … …  
    18842009        for parm in Inst:
    18852010            if parm in  ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
    1886                 'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
     2011                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1']:
    18872012                ptlbls += "%s" % (parm.center(12))
    18882013                ptstr += ptfmt % (Inst[parm][1])
    … …  
    19012026        if 'C' in dataType:
    19022027            names = ['pos','int','sig','gam']
    1903         else:
     2028        else:   #'T' and 'B'
    19042029            names = ['pos','int','alp','bet','sig','gam']
    19052030        for i,peak in enumerate(Peaks):
    … …  
    19152040        if 'C' in Inst['Type'][0]:
    19162041            names = ['pos','int','sig','gam']
    1917         else:   #'T'
     2042        else:   #'T' & 'B'
    19182043            names = ['pos','int','alp','bet','sig','gam']
    19192044        for i,peak in enumerate(Peaks):
    … …  
    19252050                if parName in varyList:
    19262051                    peak[2*j] = parmDict[parName]
    1927                 elif 'alpha' in parName:
    1928                     peak[2*j] = parmDict['alpha']/dsp
    1929                 elif 'beta' in parName:
    1930                     peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
     2052                elif 'alp' in parName:
     2053                    if 'T' in Inst['Type'][0]:
     2054                        peak[2*j] = G2mth.getTOFalpha(parmDict,dsp)
     2055                    else: #'B'
     2056                        peak[2*j] = G2mth.getPinkalpha(parmDict,pos)
     2057                elif 'bet' in parName:
     2058                    if 'T' in Inst['Type'][0]:
     2059                        peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
     2060                    else:   #'B'
     2061                        peak[2*j] = G2mth.getPinkbeta(parmDict,pos)
    19312062                elif 'sig' in parName:
    1932                     if 'C' in Inst['Type'][0]:
     2063                    if 'T' in Inst['Type'][0]:
     2064                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
     2065                    else:   #'C' & 'B'
    19332066                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
    1934                     else:
    1935                         peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
    19362067                elif 'gam' in parName:
    1937                     if 'C' in Inst['Type'][0]:
     2068                    if 'T' in Inst['Type'][0]:
     2069                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
     2070                    else:   #'C' & 'B'
    19382071                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
    1939                     else:
    1940                         peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
    19412072                       
    19422073    def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
    … …  
    19442075        if 'C' in dataType:
    19452076            names = ['pos','int','sig','gam']
    1946         else:   #'T'
     2077        else:   #'T' & 'B'
    19472078            names = ['pos','int','alp','bet','sig','gam']           
    19482079        head = 13*' '
    … …  
    19562087        if 'C' in dataType:
    19572088            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
    1958         else:
     2089        else:   #'T' & 'B'
    19592090            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
    19602091        for i,peak in enumerate(Peaks):
    … …  
    20242155        badVary = []
    20252156        result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
    2026                args=(x[xBeg:xFin],(y+fixback)[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
     2157               args=(x[xBeg:xFin],(y+fixback)[xBeg:xFin],wtFactor*w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
    20272158        ncyc = int(result[2]['nfev']/2)
    20282159        runtime = time.time()-begin   
    20292160        chisq = np.sum(result[2]['fvec']**2)
    20302161        Values2Dict(parmDict, varyList, result[0])
    2031         Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]*(y+fixback)[xBeg:xFin]**2))*100.      #to %
     2162        Rvals['Rwp'] = np.sqrt(chisq/np.sum(wtFactor*w[xBeg:xFin]*(y+fixback)[xBeg:xFin]**2))*100.      #to %
    20322163        Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
    20332164        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList)))

  • trunk/GSASIIpwdGUI.py

    r4491 r4519  
    659659        poss = x[indx]
    660660        refs = list(zip(poss,mags))
    661         if 'C' in Inst['Type'][0]:   
    662             refs = G2mth.sortArray(refs,0,reverse=True)     #small 2-Thetas first
    663         else:   #'T'OF
    664             refs = G2mth.sortArray(refs,0,reverse=False)    #big TOFs first
     661        if 'T' in Inst['Type'][0]:   
     662            refs = G2mth.sortArray(refs,0,reverse=True)     #big TOFs first
     663        else:   #'c' or 'B'
     664            refs = G2mth.sortArray(refs,0,reverse=False)    #small 2-Thetas first
    665665        for i,ref1 in enumerate(refs):      #reject picks closer than 1 FWHM
    666666            for ref2 in refs[i+1:]:
    667667                if abs(ref2[0]-ref1[0]) < 2.*G2pwd.getFWHM(ref1[0],inst):
    668668                    del(refs[i])
    669         if 'C' in Inst['Type'][0]:   
     669        if 'T' in Inst['Type'][0]:   
     670            refs = G2mth.sortArray(refs,1,reverse=False)
     671        else:   #'C' or 'B'
    670672            refs = G2mth.sortArray(refs,1,reverse=True)
    671         else:   #'T'OF
    672             refs = G2mth.sortArray(refs,1,reverse=False)
    673673        for pos,mag in refs:
    674674            data['peaks'].append(G2mth.setPeakparms(inst,inst2,pos,mag))
    … …  
    741741        finally:
    742742            dlg.Destroy()
    743        
    744743   
    745744    def OnUnDo(event):
    … …  
    837836                fixback = GetFileBackground(G2frame,data,Pattern)
    838837                peaks['sigDict'],result,sig,Rvals,varyList,parmDict,fullvaryList,badVary = G2pwd.DoPeakFit(FitPgm,peaks['peaks'],
    839                     background,limits,inst,inst2,data,fixback,prevVaryList,oneCycle,controls)
     838                    background,limits,inst,inst2,data,fixback,prevVaryList,oneCycle,controls)   #needs wtFactor after controls?
    840839                if len(result[0]) != len(fullvaryList):
    841840                    dlg.Destroy()
    … …  
    881880        Pattern = G2frame.GPXtree.GetItemPyData(PatternId)
    882881        data = Pattern[1]
     882        wtFactor = Pattern[0]['wtFactor']
    883883        bxye = GetFileBackground(G2frame,data,Pattern)
    884884        dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0,
    … …  
    890890            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
    891891        try:
    892             peaks['sigDict'] = G2pwd.DoPeakFit(FitPgm,peaks['peaks'],background,limits,inst,inst2,data,bxye,[],oneCycle,controls,dlg)[0]
     892            peaks['sigDict'] = G2pwd.DoPeakFit(FitPgm,peaks['peaks'],background,limits,inst,inst2,data,bxye,[],oneCycle,controls,wtFactor,dlg)[0]
    893893        finally:
    894894#            dlg.Destroy()
    … …  
    11041104        colLabels = ['position','refine','intensity','refine','alpha','refine',
    11051105            'beta','refine','sigma','refine','gamma','refine']
    1106         Types = [wg.GRID_VALUE_FLOAT+':10,1',wg.GRID_VALUE_BOOL,
    1107             wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
    1108             wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
    1109             wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL,
    1110             wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL,
    1111             wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL]
     1106        if 'T' in Inst['Type'][0]:
     1107            Types = [wg.GRID_VALUE_FLOAT+':10,1',wg.GRID_VALUE_BOOL,
     1108                wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
     1109                wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
     1110                wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL,
     1111                wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL,
     1112                wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL]
     1113        else:  #'B'
     1114            Types = [wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
     1115                wg.GRID_VALUE_FLOAT+':10,1',wg.GRID_VALUE_BOOL,
     1116                wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL,
     1117                wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL,
     1118                wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL,
     1119                wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL]
    11121120    T = []
    11131121    for peak in data['peaks']:
    … …  
    17881796        for key in keys:
    17891797            if key in ['Type','Bank','U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
    1790                 'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','Polariz.',
     1798                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','Polariz.','alpha-0','alpha-1',
    17911799                'Lam','Azimuth','2-theta','fltPath','difC','difA','difB','Zero','Lam1','Lam2']:
    17921800                good.append(key)
    … …  
    22242232                    if item == 'SH/L':
    22252233                        nDig = (10,5)
     2234                    labelLst.append(item)
     2235                    elemKeysLst.append([item,1])
     2236                    dspLst.append(nDig)
     2237                    refFlgElem.append([item,2])
     2238                    instSizer.Add(wx.StaticText(G2frame.dataWindow,-1,lblWdef(item,nDig[1],insDef[item])),0,WACV)
     2239                    itemVal = G2G.ValidatedTxtCtrl(G2frame.dataWindow,insVal,item,nDig=nDig,typeHint=float,OnLeave=NewProfile)
     2240                    instSizer.Add(itemVal,0,WACV)
     2241                    instSizer.Add(RefineBox(item),0,WACV)
     2242            elif 'B' in insVal['Type']:                                   #pink beam CW (x-rays & neutrons(?))
     2243                labelLst.append('Azimuth angle')
     2244                elemKeysLst.append(['Azimuth',1])
     2245                dspLst.append([10,2])
     2246                refFlgElem.append(None)                   
     2247                instSizer.Add(wx.StaticText(G2frame.dataWindow,-1,' Azimuth: '),0,WACV)
     2248                txt = '%7.2f'%(insVal['Azimuth'])
     2249                instSizer.Add(wx.StaticText(G2frame.dataWindow,-1,txt.strip()),0,WACV)
     2250                instSizer.Add((5,5),0)
     2251                key = 'Lam'
     2252                instSizer.Add(wx.StaticText(G2frame.dataWindow,-1,u' Lam (\xc5): (%10.6f)'%(insDef[key])),0,WACV)
     2253                waveVal = G2G.ValidatedTxtCtrl(G2frame.dataWindow,insVal,key,nDig=(10,6),typeHint=float,OnLeave=AfterChange)
     2254                labelLst.append(u'Lam (\xc5)')
     2255                elemKeysLst.append([key,1])
     2256                dspLst.append([10,6])
     2257                instSizer.Add(waveVal,0,WACV)
     2258                refFlgElem.append([key,2])                   
     2259                instSizer.Add(RefineBox(key),0,WACV)
     2260                for item in ['Zero','Polariz.']:
     2261                    if item in insDef:
     2262                        labelLst.append(item)
     2263                        elemKeysLst.append([item,1])
     2264                        dspLst.append([10,4])
     2265                        instSizer.Add(wx.StaticText(G2frame.dataWindow,-1,lblWdef(item,4,insDef[item])),0,WACV)
     2266                        itemVal = G2G.ValidatedTxtCtrl(G2frame.dataWindow,insVal,item,nDig=(10,4),typeHint=float,OnLeave=AfterChange)
     2267                        instSizer.Add(itemVal,0,WACV)
     2268                        refFlgElem.append([item,2])
     2269                        instSizer.Add(RefineBox(item),0,WACV)
     2270                for item in ['U','V','W','X','Y','Z','alpha-0','alpha-1','beta-0','beta-1']:
     2271                    nDig = (10,3)
    22262272                    labelLst.append(item)
    22272273                    elemKeysLst.append([item,1])
    … …  
    31463192    Inst = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))[0]
    31473193    Limits = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,G2frame.PatternId, 'Limits'))[1]
    3148     if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
     3194    if 'T' in Inst['Type'][0]:
     3195        difC = Inst['difC'][1]
     3196        dmin = G2lat.Pos2dsp(Inst,Limits[0])
     3197    else:   #'C', 'B', or 'PKS'
    31493198        wave = G2mth.getWave(Inst)
    31503199        dmin = G2lat.Pos2dsp(Inst,Limits[1])
    3151     else:
    3152         difC = Inst['difC'][1]
    3153         dmin = G2lat.Pos2dsp(Inst,Limits[0])
    31543200   
    31553201    def SetLattice(controls):
    … …  
    33683414        Inst = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))[0]
    33693415        Limits = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,G2frame.PatternId, 'Limits'))[1]
    3370         if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
     3416        if 'T' in Inst['Type'][0]:
     3417            dmin = G2lat.Pos2dsp(Inst,Limits[0])
     3418        else:
    33713419            dmin = G2lat.Pos2dsp(Inst,Limits[1])
    3372         else:
    3373             dmin = G2lat.Pos2dsp(Inst,Limits[0])
    33743420        cell = controls[6:12]
    33753421        A = G2lat.cell2A(cell)
    … …  
    35783624#        if not controls[13]: controls[13] = SPGlist[controls[5]][0]    #don't know if this is needed?   
    35793625        SGData = G2spc.SpcGroup(controls[13])[1]
    3580         if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
     3626        if 'T' in Inst['Type'][0]:
     3627            if ssopt.get('Use',False):
     3628                vecFlags = [True if x in ssopt['ssSymb'] else False for x in ['a','b','g']]
     3629                SSGData = G2spc.SSpcGroup(SGData,ssopt['ssSymb'])[1]
     3630                G2frame.HKL = G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,ssopt['ModVec'],ssopt['maxH'],A)
     3631                peaks = [G2indx.IndexSSPeaks(peaks[0],G2frame.HKL)[1],peaks[1]]   #put peak fit esds back in peaks
     3632                Lhkl,M20,X20,Aref,Vec,Zero = \
     3633                    G2indx.refinePeaksTSS(peaks[0],difC,Inst,SGData,SSGData,ssopt['maxH'],ibrav,A,ssopt['ModVec'],vecFlags,controls[1],controls[0])
     3634            else:
     3635                G2frame.HKL = G2pwd.getHKLpeak(dmin,SGData,A,Inst)
     3636                peaks = [G2indx.IndexPeaks(peaks[0],G2frame.HKL)[1],peaks[1]]   #put peak fit esds back in peaks
     3637                Lhkl,M20,X20,Aref,Zero = G2indx.refinePeaksT(peaks[0],difC,ibrav,A,controls[1],controls[0])           
     3638        else:   
    35813639            if ssopt.get('Use',False):
    35823640                vecFlags = [True if x in ssopt['ssSymb'] else False for x in ['a','b','g']]
    … …  
    35903648                peaks = [G2indx.IndexPeaks(peaks[0],G2frame.HKL)[1],peaks[1]]   #put peak fit esds back in peaks
    35913649                Lhkl,M20,X20,Aref,Zero = G2indx.refinePeaksZ(peaks[0],wave,ibrav,A,controls[1],controls[0])
    3592         else:   
    3593             if ssopt.get('Use',False):
    3594                 vecFlags = [True if x in ssopt['ssSymb'] else False for x in ['a','b','g']]
    3595                 SSGData = G2spc.SSpcGroup(SGData,ssopt['ssSymb'])[1]
    3596                 G2frame.HKL = G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,ssopt['ModVec'],ssopt['maxH'],A)
    3597                 peaks = [G2indx.IndexSSPeaks(peaks[0],G2frame.HKL)[1],peaks[1]]   #put peak fit esds back in peaks
    3598                 Lhkl,M20,X20,Aref,Vec,Zero = \
    3599                     G2indx.refinePeaksTSS(peaks[0],difC,Inst,SGData,SSGData,ssopt['maxH'],ibrav,A,ssopt['ModVec'],vecFlags,controls[1],controls[0])
    3600             else:
    3601                 G2frame.HKL = G2pwd.getHKLpeak(dmin,SGData,A,Inst)
    3602                 peaks = [G2indx.IndexPeaks(peaks[0],G2frame.HKL)[1],peaks[1]]   #put peak fit esds back in peaks
    3603                 Lhkl,M20,X20,Aref,Zero = G2indx.refinePeaksT(peaks[0],difC,ibrav,A,controls[1],controls[0])           
    36043650        controls[1] = Zero
    36053651        controls[6:12] = G2lat.A2cell(Aref)
    … …  
    47814827                colLabels = ['H','K','L','mul','d','pos','sig','gam','Fosq','Fcsq','phase','Icorr','Prfo','Trans','ExtP','I100','mustrain','Size']
    47824828                Types += 6*[wg.GRID_VALUE_FLOAT+':10,3',]
    4783             elif 'T' in Inst['Type'][0]:
     4829            else:
    47844830                colLabels = ['H','K','L','mul','d','pos','sig','gam','Fosq','Fcsq','phase','Icorr','alp','bet','wave','Prfo','Abs','Ext','I100','mustrain','Size']
    47854831                Types += 9*[wg.GRID_VALUE_FLOAT+':10,3',]

  • trunk/GSASIIstrMath.py

    r4514 r4519  
    14861486    ast = np.sqrt(np.diag(G))
    14871487    GS = G/np.outer(ast,ast)
    1488     gs = nl.inv(GS)
    14891488    uAmat,uBmat = G2lat.Gmat2AB(GS)
    14901489    Mast = twopisq*np.multiply.outer(ast,ast)   
    14911490    SGInv = SGData['SGInv']
    14921491    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
    1493     Ncen = len(SGData['SGCen'])
    14941492    Nops = len(SGMT)*(1+SGData['SGInv'])
    14951493    if SGData['SGGray']:
    … …  
    15181516        if SGData['SGGray']:
    15191517            mXYZ = np.hstack((mXYZ,mXYZ))
     1518# this done in wrong place - should be
    15201519        MmodA,MmodB = G2mth.MagMod(glTau,mXYZ,modQ,MSSdata,SGData,SSGData)  #Ntau,Nops,Natm,Mxyz cos,sim parts sum matches drawing
    15211520        Mmod = MmodA+MmodB
    … …  
    15301529            Mmod += GSdata[nxs,:,:,:]
    15311530
     1531#        Kdata = np.inner(Mmod,uAmat)          #Ntau,Nops,Natm,Mxyz
     1532#        Kmean = np.mean(np.sqrt(np.sum(Kdata**2,axis=0)),axis=0)
     1533#        Kdata /= Kmean     #Ntau,Nops,Natm,Mxyz  Cartesian unit vectors along mag moments
     1534       
    15321535    FF = np.zeros(len(Tdata))
    15331536    if 'NC' in calcControls[hfx+'histType']:
    … …  
    15971600        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:           
    15981601           
    1599             MmodA,MmodB = G2mth.MagMod2(H[3],glTau,mXYZ,modQ,MSSdata,SGData,SSGData)  #Nref,Ntau,Nops,Natm,Mxyz
    16001602            phasem = twopi*np.inner(mXYZ,HP.T).T    #2pi(Q.r)
    16011603            cosm = np.cos(phasem)                   #Nref,nops,natm
    … …  
    16031605            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
    16041606            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Mdata*Fdata*MF/(2*Nops)     #Nref,Natm
    1605             HM = np.inner(uAmat,HP.T)                            #put into cartesian space X||H,Z||H*L; uAmat better than uBmat
    1606             eM = (HM/np.sqrt(np.sum(HM**2,axis=0))).T              # normalize  HP  Nref,hkl=Unit vectors || Q
     1607            HM = np.inner(uAmat,HP.T)                    #put into cartesian space X||H,Z||H*L; uAmat better than uBmat
     1608            eM = (HM/np.sqrt(np.sum(HM**2,axis=0))).T    # normalize  HP  Nref,hkl=Unit vectors || Q
     1609#            eDotK = np.sum(eM[:,nxs,nxs,nxs,:]*Kdata[nxs,:,:,:,:],axis=-1)
     1610#            QC = eM[:,nxs,nxs,nxs,:]*eDotK[:,:,:,:,nxs]-Kdata[nxs,:,:,:,:] #Nref,Ntau,Nop,Natm,xyz; Cartesian
     1611#            Q = np.inner(QC,uBmat.T)    #crystal space; Nref,Ntau,Nop,Natm,xyz .T?
    16071612
    16081613            fam0 = 0.
    … …  
    16111616                fam0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*cosm[:,:,:,nxs]    #Nref,Nops,Natm,Mxyz
    16121617                fbm0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*sinm[:,:,:,nxs]
    1613                            
    1614             fams = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA[i]*cosm[i,nxs,:,:,nxs]+   
    1615                 np.sign(H[3,i])*MmodB[i]*sinm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])/2.          #Nref,Ntau,Nops,Natm,Mxyz
     1618#  calc mag. structure factors; Nref,Ntau,Nops,Natm,Mxyz                            
     1619            fams = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*cosm[i,nxs,:,:,nxs]+   
     1620                H[3,i]*MmodB*sinm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
    16161621                       
    1617             fbms = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA[i]*sinm[i,nxs,:,:,nxs]+   
    1618                 np.sign(H[3,i])*MmodB[i]*cosm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])/2.          #Nref,Ntau,Nops,Natm,Mxyz
    1619            
     1622            fbms = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*sinm[i,nxs,:,:,nxs]+   
     1623                H[3,i]*MmodB*cosm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
     1624
    16201625            if not SGData['SGGray']:
    16211626                fams += fam0[:,nxs,:,:,:]
    16221627                fbms += fbm0[:,nxs,:,:,:]
    1623                                
    1624 # # do sum on ops, atms 1st                       
     1628               
     1629# this is the best, but not exactly right
     1630#sum ops & atms                               
    16251631            fasm = np.sum(np.sum(fams,axis=-2),axis=-2)    #Nref,Ntau,Mxyz; sum ops & atoms
    16261632            fbsm = np.sum(np.sum(fbms,axis=-2),axis=-2)
    1627 # #put into cartesian space
    1628             facm = np.inner(fasm,uBmat.T)       #uBmat better than uAmat
     1633#put into cartesian space
     1634            facm = np.inner(fasm,uBmat.T)
    16291635            fbcm = np.inner(fbsm,uBmat.T)
    1630 # #form e.F dot product
     1636#form e.F dot product
    16311637            eDotFa = np.sum(eM[:,nxs,:]*facm,axis=-1)    #Nref,Ntau       
    16321638            eDotFb = np.sum(eM[:,nxs,:]*fbcm,axis=-1)
    1633 # #intensity    Halpern & Johnson
     1639#intensity Halpern & Johnson
     1640#            fass = np.sum((facm-eM[:,nxs,:]*eDotFa[:,:,nxs])**2,axis=-1)
     1641#            fbss = np.sum((fbcm-eM[:,nxs,:]*eDotFb[:,:,nxs])**2,axis=-1)
    16341642            fass = np.sum(facm**2,axis=-1)-eDotFa**2
    16351643            fbss = np.sum(fbcm**2,axis=-1)-eDotFb**2
    1636            
    1637 # #do integration           
     1644## #do integration           
    16381645            fas = np.sum(fass*glWt[nxs,:],axis=1)
    16391646            fbs = np.sum(fbss*glWt[nxs,:],axis=1)
Note: See TracChangeset for help on using the changeset viewer.