# Changeset 803

Ignore:
Timestamp:
Nov 21, 2012 1:30:38 PM (10 years ago)
Message:

GSASIImath.py: map search gives mag in name, & change a looping form in findOffset
GSASIIphsGUI.py: map search gives mag in name
GSASIIplot.py: more curves for instrument parms plot
GSASIIpwd.py: make cw array of channel widths & use them
GSASIIstruct.py: ditto
G2pwd_fxye.py: fix for zero sig

Location:
trunk
Files:
6 edited

Unmodified
Removed
• ## trunk/GSASIImath.py

 r802 Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]] ang0 = np.angle(Fh0,deg=True)/360. for j,H in enumerate(Uniq[1:]): ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-Phi[j+1]) for H,phi in zip(Uniq,Phi)[1:]: ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi) dH = H-hkl dang = ang-ang0 continue DH.append(dH) Dphi.append((dang+0.5) % 1.0) Dphi.append((dang+.5) % 1.0) i += 1 DH = np.array(DH) X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]] XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten())) Mmap = np.reshape(np.sum(((np.dot(XYZ,DH.T)+.5)%1.-Dphi)**2,axis=1),newshape=steps) Mmap = np.reshape(np.sum(((np.dot(XYZ,DH.T)+.5)%1.-Dphi)**2,axis=1),newshape=steps)/len(DH) #    hist,bins = np.histogram(Mmap,bins=1000) #    for i,item in enumerate(hist[:10]): #        print item,bins[i] chisq = np.min(Mmap) DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape)) SQ = 0.25/dsp**2 ff *= G2el.ScatFac(FFtable,SQ)[0] if ref[8] > 0.: if ref[8] > 0.:         #use only +ve Fobs**2 E = np.sqrt(ref[8])/ff else: return Ind def setPeakparms(Parms,Parms2,pos,mag,ifQ=False): def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False): ind = 0 if useFit: ind = 1 ins = {} if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif for x in ['U','V','W','X','Y']: ins[x] = Parms[x][0] ins[x] = Parms[x][ind] if ifQ:                              #qplot - convert back to 2-theta pos = 2.0*asind(pos*wave/(4*math.pi)) if 'Pdabc' in Parms2: for x in ['sig-0','sig-1','X','Y']: ins[x] = Parms[x][0] ins[x] = Parms[x][ind] Pdabc = Parms2['Pdabc'].T alp = np.interp(dsp,Pdabc[0],Pdabc[1]) else: for x in ['alpha','beta-0','beta-1','sig-0','sig-1','X','Y']: ins[x] = Parms[x][0] ins[x] = Parms[x][ind] alp = ins['alpha']/dsp bet = ins['beta-0']+ins['beta-1']/dsp**4
• ## trunk/GSASIIphsGUI.py

 r797 event.StopPropagation() def AtomAdd(x,y,z,El='H'): def AtomAdd(x,y,z,El='H',Name='UNK'): atomData = data['Atoms'] generalData = data['General'] Sytsym,Mult = G2spc.SytSym([x,y,z],SGData) if generalData['Type'] == 'macromolecular': atomData.append([0,'UNK','','UNK',El,'',x,y,z,1,Sytsym,Mult,'I',0.10,0,0,0,0,0,0,atId]) atomData.append([0,Name,'',Name,El,'',x,y,z,1,Sytsym,Mult,'I',0.10,0,0,0,0,0,0,atId]) elif generalData['Type'] == 'nuclear': atomData.append(['UNK',El,'',x,y,z,1,Sytsym,Mult,'I',0.01,0,0,0,0,0,0,atId]) atomData.append([Name,El,'',x,y,z,1,Sytsym,Mult,'I',0.01,0,0,0,0,0,0,atId]) elif generalData['Type'] == 'magnetic': atomData.append(['UNK',El,'',x,y,z,1,Sytsym,Mult,0,'I',0.01,0,0,0,0,0,0,0,0,0,atId]) atomData.append([Name,El,'',x,y,z,1,Sytsym,Mult,0,'I',0.01,0,0,0,0,0,0,0,0,0,atId]) SetupGeneral() if 'Atoms' in data['Drawing']: PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,HistoNames[0]) refData = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))[PhaseName] wx.BeginBusyCursor() try: for iref,ref in enumerate(Refs): try: ref[6] = abs(refData[iref][9]) ref[7] = 1.0 if refData[iref][9] < 0.: ref[6] = abs(refData[iref][9]) ref[7] = 1.0 except IndexError: pass def OnPeaksMove(event): if 'Map Peaks' in data: mapPeaks = data['Map Peaks'] mapPeaks = np.array(data['Map Peaks']) peakMax = np.max(mapPeaks.T[0]) Ind = MapPeaks.GetSelectedRows() for ind in Ind: x,y,z,d = mapPeaks[ind][1:] AtomAdd(x,y,z,'C') mag,x,y,z,d = mapPeaks[ind] AtomAdd(x,y,z,'C',Name='M '+'%d'%(int(100*mag/peakMax))) def OnPeaksClear(event):
• ## trunk/GSASIIplot.py

 r800 if ind.all() != [0]:                                    #picked a data point data = G2frame.PatternTree.GetItemPyData(G2frame.PickId) XY = G2mth.setPeakparms(Parms,Parms2,xy[0],xy[1]) XY = G2mth.setPeakparms(Parms,Parms2,xy[0],xy[1])           #what happens for a q-plot??? data.append(XY) G2pdG.UpdatePeakGrid(G2frame,data) Seen when "Instrument Parameters" chosen from powder pattern data tree ''' sig = lambda Th,U,V,W: 1.17741*math.sqrt(U*tand(Th)**2+V*tand(Th)+W)*math.pi/18000. gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/18000. #    sig = lambda Th,U,V,W: 1.17741*math.sqrt(U*tand(Th)**2+V*tand(Th)+W)*math.pi/18000. #    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/18000. gamFW = lambda s,g: np.exp(np.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.) #    gamFW2 = lambda s,g: math.sqrt(s**2+(0.4654996*g)**2)+.5345004*g  #Ubaldo Bafile - private communication Plot.plot(Q,W,color='b',label='G+L') fit = G2mth.setPeakparms(Parms,Parms2,X,Z,useFit=True) sf = 1.17741*np.sqrt(fit[4])*np.pi/18000. gf = fit[6]*np.pi/18000. Gf = gamFW(gf,sf) Yf = sf/nptand(X/2.) Zf = gf/nptand(X/2.) Wf = Gf/nptand(X/2.) Plot.plot(Q,Yf,color='r',dashes=(5,5),label='Gaussian fit') Plot.plot(Q,Zf,color='g',dashes=(5,5),label='Lorentzian fit') Plot.plot(Q,Wf,color='b',dashes=(5,5),label='G+L fit') X = [] Y = [] Z.append(g/tand(peak[0]/2.)) W.append(G/tand(peak[0]/2.)) Plot.plot(X,Y,'+',color='r',label='G peak') Plot.plot(X,Z,'+',color='g',label='L peak') Plot.plot(X,W,'+',color='b',label='G+L peak') if len(peaks): Plot.plot(X,Y,'+',color='r',label='G peak') Plot.plot(X,Z,'+',color='g',label='L peak') Plot.plot(X,W,'+',color='b',label='G+L peak') Plot.legend(loc='best') Page.canvas.draw() Z = np.ones_like(T) data = G2mth.setPeakparms(Parms,Parms2,T,Z) #       = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0] ds = T/difC Q = 2.*np.pi/ds Plot.plot(Q,S,color='b',label='Gaussian') Plot.plot(Q,G,color='m',label='Lorentzian') fit = G2mth.setPeakparms(Parms,Parms2,T,Z) ds = T/difC Q = 2.*np.pi/ds Af = fit[4] Bf = fit[6] Sf = 1.17741*np.sqrt(fit[8])/T Gf = fit[10]/T Plot.plot(Q,Af,color='r',dashes=(5,5),label='Alpha fit') Plot.plot(Q,Bf,color='g',dashes=(5,5),label='Beta fit') Plot.plot(Q,Sf,color='b',dashes=(5,5),label='Gaussian fit') Plot.plot(Q,Gf,color='m',dashes=(5,5),label='Lorentzian fit') T = [] A = []
• ## trunk/GSASIIpwd.py

 r801 yb = np.zeros_like(xdata) nBak = 0 cw = np.diff(xdata) cw = np.append(cw,cw[-1]) while True: key = pfx+'Back:'+str(nBak) dyddb = np.zeros(shape=(3*parmDict[pfx+'nDebye'],len(xdata))) dydpk = np.zeros(shape=(4*parmDict[pfx+'nPeaks'],len(xdata))) dx = xdata[1]-xdata[0] cw = np.diff(xdata) cw = np.append(cw,cw[-1]) if bakType in ['chebyschev','cosine']: bakPos[0] = xdata[0] bakPos[-1] = xdata[-1] dx = bakPos[1]-bakPos[0] for i,pos in enumerate(bakPos): if i == 0: cqr = np.cos(q*dbR) temp = np.exp(-dbU*q**2) dyddb[3*iD] = ff*sqr*temp/(np.pi*dx) dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(np.pi*dbR*dx) dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2/(np.pi*dx) dyddb[3*iD] = ff*sqr*temp/(np.pi*cw) dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(np.pi*dbR*cw) dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2/(np.pi*cw) iD += 1       #ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR) except KeyError: iFin = np.searchsorted(xdata,pkP+fmax) Df,dFdp,dFds,dFdg,dFdsh = getdFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin]) dydpk[4*iD][iBeg:iFin] += 100.*dx*pkI*dFdp dydpk[4*iD+1][iBeg:iFin] += 100.*dx*Df dydpk[4*iD+2][iBeg:iFin] += 100.*dx*pkI*dFds dydpk[4*iD+3][iBeg:iFin] += 100.*dx*pkI*dFdg dydpk[4*iD][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdp dydpk[4*iD+1][iBeg:iFin] += 100.*cw[iBeg:iFin]*Df dydpk[4*iD+2][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFds dydpk[4*iD+3][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdg iD += 1 except KeyError: yb = getBackground('',parmDict,bakType,xdata) yc = np.zeros_like(yb) cw = np.diff(xdata) cw = np.append(cw,cw[-1]) if 'C' in dataType: dx = xdata[1]-xdata[0] U = parmDict['U'] V = parmDict['V'] Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl) iBeg = np.searchsorted(xdata,pos-fmin) lenX = len(xdata) if not iBeg: iFin = np.searchsorted(xdata,pos+fmin) elif iBeg == lenX: iFin = iBeg else: iFin = min(lenX,iBeg+int((fmin+fmax)/dx)) iFin = np.searchsorted(xdata,pos+fmin) if not iBeg+iFin:       #peak below low limit iPeak += 1 if Ka2: pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th) kdelt = int((pos2-pos)/dx) iBeg = min(lenX,iBeg+kdelt) iFin = min(lenX,iFin+kdelt) iBeg = np.searchsorted(xdata,pos2-fmin) iFin = np.searchsorted(xdata,pos2+fmin) if iBeg-iFin: yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin]) ip = names.index(parm) dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip] cw = np.diff(xdata) cw = np.append(cw,cw[-1]) if 'C' in dataType: dx = xdata[1]-xdata[0] U = parmDict['U'] V = parmDict['V'] Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl) iBeg = np.searchsorted(xdata,pos-fmin) lenX = len(xdata) if not iBeg: iFin = np.searchsorted(xdata,pos+fmin) elif iBeg == lenX: iFin = iBeg else: iFin = min(lenX,iBeg+int((fmin+fmax)/dx)) iFin = np.searchsorted(xdata,pos+fmin) if not iBeg+iFin:       #peak below low limit iPeak += 1 dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin]) for i in range(1,5): dMdpk[i][iBeg:iFin] += 100.*dx*intens*dMdipk[i] dMdpk[0][iBeg:iFin] += 100.*dx*dMdipk[0] dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i] dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0] dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]} if Ka2: pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th) kdelt = int((pos2-pos)/dx) iBeg = min(lenX,iBeg+kdelt) iFin = min(lenX,iFin+kdelt) iBeg = np.searchsorted(xdata,pos2-fmin) iFin = np.searchsorted(xdata,pos2+fmin) if iBeg-iFin: dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin]) for i in range(1,5): dMdpk[i][iBeg:iFin] += 100.*dx*intens*kRatio*dMdipk2[i] dMdpk[0][iBeg:iFin] += 100.*dx*kRatio*dMdipk2[0] dMdpk[5][iBeg:iFin] += 100.*dx*dMdipk2[0] dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i] dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0] dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0] dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens} for parmName in ['pos','int','sig','gam']: dMdpk = np.zeros(shape=(7,len(xdata))) dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin]) cw = np.diff(xdata[iBeg:iFin]) cw = np.append(cw,cw[-1]) for i in range(1,6): dMdpk[i][iBeg:iFin] += intens*cw*dMdipk[i] dMdpk[0][iBeg:iFin] += cw*dMdipk[0] dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i] dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0] dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]} for parmName in ['pos','int','alp','bet','sig','gam']:
• ## trunk/GSASIIstruct.py

 r796 import time import math import random import cPickle import numpy as np pId = int(item.split(':')[0]) pNames.append(item) pVals.append(-negWt[pId]*parmDict[item]**2) #            pVals.append(negWt[pId]*parmDict[item]**2) pVals.append(-negWt[pId]*parmDict[item]) pVals = np.array(pVals) return pNames,pVals if item in pNames: pId = int(item.split(':')[0]) pDerv[i][pNames.index(item)] -= 2.*negWt[pId]*parmDict[item] #            pDerv[i][pNames.index(item)] += 2.*negWt[pId]*parmDict[item] pDerv[i][pNames.index(item)] += negWt[pId] return pDerv try: pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]) #                        parmDict[pInd] = max(parmDict[pInd]/2.,parmDict[pInd]) refl[9] = parmDict[pInd] except KeyError: ip = names.index(parm) dMdv[varylist.index(name)] = dMdpk[4*id+ip] cw = np.diff(x) cw = np.append(cw,cw[-1]) if 'C' in calcControls[hfx+'histType']: dx = x[1]-x[0] shl = max(parmDict[hfx+'SH/L'],0.002) Ka2 = False dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) for i in range(1,5): dMdpk[i] += 100.*dx*refl[13]*refl[9]*dMdipk[i] dMdpk[0] += 100.*dx*refl[13]*refl[9]*dMdipk[0] dMdpk[i] += 100.*cw[iBeg:iFin]*refl[13]*refl[9]*dMdipk[i] dMdpk[0] += 100.*cw[iBeg:iFin]*refl[13]*refl[9]*dMdipk[0] dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])} if Ka2: pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th) kdelt = int((pos2-refl[5])/dx) iBeg2 = min(lenX,iBeg+kdelt) iFin2 = min(lenX,iFin+kdelt) iBeg2 = np.searchsorted(x,pos2-fmin) iFin2 = np.searchsorted(x,pos2+fmax) if iBeg2-iFin2: lenBF2 = iFin2-iBeg2 dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) for i in range(1,5): dMdpk2[i] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i] dMdpk2[0] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0] dMdpk2[5] = 100.*dx*refl[13]*dMdipk2[0] dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[13]*refl[9]*kRatio*dMdipk2[i] dMdpk2[0] = 100.*cw[iBeg2:iFin2]*refl[13]*refl[9]*kRatio*dMdipk2[0] dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[13]*dMdipk2[0] dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]} if Phase['General'].get('doPawley'): pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]) idx = varylist.index(pIdx) #                        parmDict[pIdx] = max(parmDict[pIdx]/2.,parmDict[pIdx]) #                        refl[9] = parmDict[pIdx] dMdpw[iBeg:iFin] = dervDict['int']/refl[9] #                        if parmDict[pIdx] < 0.: #                            dMdpw[iBeg:iFin] = 2.*dervDict['int']/refl[9] if Ka2: dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9] #                            if parmDict[pIdx] < 0.: #                                dMdpw[iBeg2:iFin2] += 2.*dervDict['int']/refl[9] dMdv[idx] = dMdpw except: # ValueError:
• ## trunk/imports/G2pwd_fxye.py

 r795 x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees f = float(vals[1]) if f <= 0.0: s = float(vals[2]) if f <= 0.0 or s <= 0.0: y.append(0.0) w.append(1.0)
Note: See TracChangeset for help on using the changeset viewer.