# Changeset 353

Ignore:
Timestamp:
Aug 24, 2011 4:07:29 PM (11 years ago)
Message:

modify psvfcj.for to use sh/l only
GSASIIgrid.py - new LS controls for derivative type
GSASIIlattice.py - work on docs as per P Jemian
GSASIIpwd.py - major mods to use FCJpsvoight fortran code & derivatives
GSASIIpwdGUI.py - mod to use LS controls

Location:
trunk
Files:
6 edited

Unmodified
Removed

• ## trunk/GSASIIlattice.py

 r342 def sec2HMS(sec): """Convert time in sec to H:M:S string :param sec: time in seconds return: H:M:S string (to nearest 100th second) """ H = int(sec/3600) M = int(sec/60-H*60) def rotdMat(angle,axis=0): '''Prepare rotation matrix for angle in degrees about axis(=0,1,2) Returns numpy 3,3 array ''' """Prepare rotation matrix for angle in degrees about axis(=0,1,2) :param angle: angle in degrees :param axis:  axis (0,1,2 = x,y,z) about which for the rotation :return: rotation matrix - 3x3 numpy array """ if axis == 2: return np.array([[cosd(angle),-sind(angle),0],[sind(angle),cosd(angle),0],[0,0,1]]) def rotdMat4(angle,axis=0): '''Prepare rotation matrix for angle in degrees about axis(=0,1,2) with scaling for OpenGL Returns numpy 4,4 array ''' """Prepare rotation matrix for angle in degrees about axis(=0,1,2) with scaling for OpenGL :param angle: angle in degrees :param axis:  axis (0,1,2 = x,y,z) about which for the rotation :return: rotation matrix - 4x4 numpy array (last row/column for openGL scaling) """ Mat = rotdMat(angle,axis) return np.concatenate((np.concatenate((Mat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0) def fillgmat(cell): '''Compute lattice metric tensor from unit cell constants cell is tuple with a,b,c,alpha, beta, gamma (degrees) returns 3x3 numpy array ''' """Compute lattice metric tensor from unit cell constants :param cell: tuple with a,b,c,alpha, beta, gamma (degrees) :return: 3x3 numpy array """ a,b,c,alp,bet,gam = cell g = np.array([ def cell2Gmat(cell): '''Compute real and reciprocal lattice metric tensor from unit cell constants cell is tuple with a,b,c,alpha, beta, gamma (degrees) returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays) ''' """Compute real and reciprocal lattice metric tensor from unit cell constants :param cell: tuple with a,b,c,alpha, beta, gamma (degrees) :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays) """ g = fillgmat(cell) G = nl.inv(g) def A2Gmat(A): '''Fill reciprocal metric tensor (G) from A returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays) ''' """Fill real & reciprocal metric tensor (G) from A :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23] :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays) """ G = np.zeros(shape=(3,3)) G = [ def Gmat2A(G): 'Extract A from reciprocal metric tensor (G)' """Extract A from reciprocal metric tensor (G) :param G: reciprocal maetric tensor (3x3 numpy array :return: A = [G11,G22,G33,2*G12,2*G13,2*G23] """ return [G[0][0],G[1][1],G[2][2],2.*G[0][1],2.*G[0][2],2.*G[1][2]] def cell2A(cell): """Obtain A = [G11,G22,G33,2*G12,2*G13,2*G23] from lattice parameters :param cell: [a,b,c,alpha,beta,gamma] (degrees) :return: G reciprocal metric tensor as 3x3 numpy array """ G,g = cell2Gmat(cell) return Gmat2A(G) def A2cell(A): '''Compute unit cell constants from A tensor returns tuple with a,b,c,alpha, beta, gamma (degrees) ''' """Compute unit cell constants from A :param A: [G11,G22,G33,2*G12,2*G13,2*G23] G - reciprocal metric tensor :return: a,b,c,alpha, beta, gamma (degrees) - lattice parameters """ G,g = A2Gmat(A) return Gmat2cell(g) def Gmat2cell(g): '''Compute lattice parameters from real metric tensor (g) returns tuple with a,b,c,alpha, beta, gamma (degrees) Alternatively,compute reciprocal lattice parameters from inverse metric tensor (G) returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees) ''' """Compute real/reciprocal lattice parameters from real/reciprocal metric tensor (g/G) The math works the same either way. :param g (or G): real (or reciprocal) metric tensor 3x3 array :return: a,b,c,alpha, beta, gamma (degrees) (or a*,b*,c*,alpha*,beta*,gamma* degrees) """ oldset = np.seterr('raise') a = np.sqrt(max(0,g[0][0])) def invcell2Gmat(invcell): '''Compute real and reciprocal lattice metric tensor from reciprocal """Compute real and reciprocal lattice metric tensor from reciprocal unit cell constants invcell is tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees) returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays) ''' :param invcell: [a*,b*,c*,alpha*, beta*, gamma*] (degrees) :return: reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays) """ G = fillgmat(invcell) g = nl.inv(G) def calc_rVsq(A): 'Compute the square of the reciprocal lattice volume (V* **2) from A' """Compute the square of the reciprocal lattice volume (V* **2) from A' """ G,g = A2Gmat(A) rVsq = nl.det(G) def calc_rV(A): 'Compute the reciprocal lattice volume (V*) from A' """Compute the reciprocal lattice volume (V*) from A """ return np.sqrt(calc_rVsq(A)) def calc_V(A): 'Compute the real lattice volume (V) from A' """Compute the real lattice volume (V) from A """ return 1./calc_rV(A) def A2invcell(A): '''Compute reciprocal unit cell constants from A """Compute reciprocal unit cell constants from A returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees) ''' """ G,g = A2Gmat(A) return Gmat2cell(G) def cell2AB(cell): '''Computes orthogonalization matrix from unit cell constants """Computes orthogonalization matrix from unit cell constants cell is tuple with a,b,c,alpha, beta, gamma (degrees) returns tuple of two 3x3 numpy arrays (A,B) A for crystal to Cartesian transformations A*x = np.inner(A,x) = X B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B*x) = x ''' """ G,g = cell2Gmat(cell) cellstar = Gmat2cell(G) def U6toUij(U6): '''Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23] """Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23] returns ''' """ U = np.zeros(shape=(3,3)) U = [ def Uij2betaij(Uij,G): ''' """ Convert Uij to beta-ij tensors input: returns: beta-ij - numpy array [beta-ij] ''' """ pass def CosSinAngle(U,V,G): ''' calculate sin & cos of angle betwee U & V in generalized coordinates """ calculate sin & cos of angle betwee U & V in generalized coordinates defined by metric tensor G input: return: cos(phi) & sin(phi) ''' """ u = U/nl.norm(U) v = V/nl.norm(V) def criticalEllipse(prob): ''' """ Calculate critical values for probability ellipsoids from probability ''' """ if not ( 0.01 <= prob < 1.0): return 1.54 def CellBlock(nCells): ''' """ Generate block of unit cells n*n*n on a side; [0,0,0] centered, n = 2*nCells+1 currently only works for nCells = 0 or 1 (not >1) ''' """ if nCells: N = 2*nCells+1 # def _combinators(_handle, items, n): ''' factored-out common structure of all following combinators ''' """ factored-out common structure of all following combinators """ if n==0: yield [ ] yield this_one + cc def combinations(items, n): ''' take n distinct items, order matters ''' """ take n distinct items, order matters """ def skipIthItem(items, i): return items[:i] + items[i+1:] return _combinators(skipIthItem, items, n) def uniqueCombinations(items, n): ''' take n distinct items, order is irrelevant ''' """ take n distinct items, order is irrelevant """ def afterIthItem(items, i): return items[i+1:] return _combinators(afterIthItem, items, n) def selections(items, n): ''' take n (not necessarily distinct) items, order matters ''' """ take n (not necessarily distinct) items, order matters """ def keepAllItems(items, i): return items return _combinators(keepAllItems, items, n) def permutations(items): ''' take all items, order matters ''' """ take all items, order matters """ return combinations(items, len(items)) def GetBraviasNum(center,system): '''Determine the Bravais lattice number, as used in GenHBravais center = one of: P, C, I, F, R (see SGLatt from GSASIIspc.SpcGroup) lattice = is cubic, hexagonal, tetragonal, orthorhombic, trigonal (R) monoclinic, triclinic (see SGSys from GSASIIspc.SpcGroup) Returns a number between 0 and 13 or throws an exception if the setting is non-standard ''' """Determine the Bravais lattice number, as used in GenHBravais :param center: one of: 'P', 'C', 'I', 'F', 'R' (see SGLatt from GSASIIspc.SpcGroup) :param system: one of 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal' (for R) 'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup) :return: a number between 0 and 13 or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard) """ if center.upper() == 'F' and system.lower() == 'cubic': return 0 def GenHBravais(dmin,Bravais,A): '''Generate the positionally unique powder diffraction reflections input: dmin is minimum d-space Bravais is 0-13 to indicate lattice type (see GetBraviasNum) A is reciprocal cell tensor (see Gmat2A or cell2A) returns: a list of tuples containing: h,k,l,d-space,-1 ''' # Bravais in range(14) to indicate Bravais lattice: #   0 F cubic #   1 I cubic #   2 P cubic #   3 R hexagonal (trigonal not rhombohedral) #   4 P hexagonal #   5 I tetragonal #   6 P tetragonal #   7 F orthorhombic #   8 I orthorhombic #   9 C orthorhombic #  10 P orthorhombic #  11 C monoclinic #  12 P monoclinic #  13 P triclinic # A - as defined in calc_rDsq # returns HKL = [h,k,l,d,0] sorted so d largest first """Generate the positionally unique powder diffraction reflections :param dmin: minimum d-spacing in A :param Bravais: lattice type (see GetBraviasNum) Bravais is one of:: 0 F cubic 1 I cubic 2 P cubic 3 R hexagonal (trigonal not rhombohedral) 4 P hexagonal 5 I tetragonal 6 P tetragonal 7 F orthorhombic 8 I orthorhombic 9 C orthorhombic 10 P orthorhombic 11 C monoclinic 12 P monoclinic 13 P triclinic :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23] :return: HKL unique d list of [h,k,l,d,-1] sorted with largest d first """ import math if Bravais in [9,11]: def GenHLaue(dmin,SGData,A): '''Generate the crystallographically unique powder diffraction reflections """Generate the crystallographically unique powder diffraction reflections for a lattice and Bravais type Input: dmin - minimum d-spacing SGData - space group dictionary with at least: SGLaue - Laue group symbol = '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m' SGLatt - lattice centering = 'P','A','B','C','I','F' SGUniq - code for unique monoclinic axis = 'a','b','c' A - 6 terms as defined in calc_rDsq Return; HKL = list of [h,k,l,d] sorted with largest d first and is unique :param dmin: minimum d-spacing :param SGData: space group dictionary with at least:: 'SGLaue': Laue group symbol: one of '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m' 'SGLatt': lattice centering: one of 'P','A','B','C','I','F' 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise ' ' :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23] :return: HKL = list of [h,k,l,d] sorted with largest d first and is unique part of reciprocal space ignoring anomalous dispersion ''' """ import math SGLaue = SGData['SGLaue']
• ## trunk/GSASIIpwd.py

 r350 fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx') def getBackground(pfx,parmDict,bakType,xdata): yb = np.zeros_like(xdata) if bakType == 'chebyschev': iBak = 0 while True: key = pfx+'Back:'+str(iBak) try: yb += parmDict[key]*(xdata-xdata[0])**iBak iBak += 1 except KeyError: break return yb def getWidths(pos,sig,gam,shl): widths = [np.sqrt(sig)/100.,gam/200.] return intens*Df(xdata)*DX/dx def getBackground(pfx,parmDict,bakType,xdata): yb = np.zeros_like(xdata) if bakType == 'chebyschev': iBak = 0 while True: key = pfx+'Back:'+str(iBak) try: yb += parmDict[key]*(xdata-xdata[0])**iBak iBak += 1 except KeyError: break return yb def getBackgroundDerv(pfx,parmDict,bakType,xdata): dydb = [] if bakType == 'chebyschev': iBak = 0 while True: if pfx+'Back:'+str(iBak) in parmDict: dydb.append((xdata-xdata[0])**iBak) iBak += 1 else: break return dydb #use old fortran routine def getFCJVoigt3(pos,intens,sig,gam,shl,xdata): def getFCJVoigt3(pos,sig,gam,shl,xdata): Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl) Df /= np.sum(Df) return intens*Df return Df def getdFCJVoigt3(pos,sig,gam,shl,xdata): Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl) #might have to make these numpy arrays? sumDf = np.sum(Df) return Df/sumDf,dFdp,dFds,dFdg,dFdsh def getPeakProfile(parmDict,xdata,varyList,bakType): X = parmDict['X'] Y = parmDict['Y'] shl = max(parmDict['SH/L'],0.0005) shl = max(parmDict['SH/L'],0.002) Ka2 = False if 'Lam1' in parmDict.keys(): elif not iBeg-iFin:     #peak above high limit return yb+yc yc[iBeg:iFin] += getFCJVoigt3(pos,intens,sig,gam,shl,xdata[iBeg:iFin]) yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin]) if Ka2: pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th) iFin = min(lenX,iFin+kdelt) if iBeg-iFin: yc[iBeg:iFin] += getFCJVoigt3(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin]) yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin]) iPeak += 1 except KeyError:        #no more peaks to process return yb+yc def getPeakProfileDerv(parmDict,xdata,varyList,bakType): # needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order dMdv = np.zeros(shape=(len(varyList),len(xdata))) if 'Back:0' in varyList:            #background derivs are in front if present dMdb = getBackgroundDerv('',parmDict,bakType,xdata) dMdv[0:len(dMdb)] = dMdb dx = xdata[1]-xdata[0] U = parmDict['U'] V = parmDict['V'] W = parmDict['W'] X = parmDict['X'] Y = parmDict['Y'] shl = max(parmDict['SH/L'],0.002) Ka2 = False if 'Lam1' in parmDict.keys(): Ka2 = True lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1']) kRatio = parmDict['I(L2)/I(L1)'] iPeak = 0 while True: try: pos = parmDict['pos'+str(iPeak)] intens = parmDict['int'+str(iPeak)] sigName = 'sig'+str(iPeak) tanth = tand(pos/2.0) costh = cosd(pos/2.0) if sigName in varyList: sig = parmDict[sigName] else: sig = U*tanth**2+V*tanth+W dsdU = tanth**2 dsdV = tanth dsdW = 1.0 sig = max(sig,0.001)          #avoid neg sigma gamName = 'gam'+str(iPeak) if gamName in varyList: gam = parmDict[gamName] else: gam = X/costh+Y*tanth dgdX = 1.0/costh dgdY = tanth gam = max(gam,0.001)             #avoid neg gamma Wd,fmin,fmax = getWidths(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)) if not iBeg+iFin:       #peak below low limit iPeak += 1 continue elif not iBeg-iFin:     #peak above high limit break dMdpk = np.zeros(shape=(6,len(xdata))) dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin]) for i in range(5): dMdpk[i][iBeg:iFin] = 100.*dx*intens*dMdipk[i] dMdpk[0][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) if iBeg-iFin: dMdipk = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin]) for i in range(5): dMdpk[i][iBeg:iFin] += 100.*dx*intens*kRatio*dMdipk[i] dMdpk[0][iBeg:iFin] += kRatio*dMdipk[0] dMdpk[5][iBeg:iFin] += dMdipk[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']: try: idx = varyList.index(parmName+str(iPeak)) dMdv[idx] = dervDict[parmName] except ValueError: pass if 'U' in varyList: dMdv[varyList.index('U')] += dsdU*dervDict['sig'] if 'V' in varyList: dMdv[varyList.index('V')] += dsdV*dervDict['sig'] if 'W' in varyList: dMdv[varyList.index('W')] += dsdW*dervDict['sig'] if 'X' in varyList: dMdv[varyList.index('X')] += dgdX*dervDict['gam'] if 'Y' in varyList: dMdv[varyList.index('Y')] += dgdY*dervDict['gam'] if 'SH/L' in varyList: dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here if 'I(L2)/I(L1)' in varyList: dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2'] iPeak += 1 except KeyError:        #no more peaks to process break return dMdv def Dict2Values(parmdict, varylist): '''Use before call to leastsq to setup list of values for the parameters parmdict.update(zip(varylist,values)) def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False): def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False,controls=None): def SetBackgroundParms(Background): insVary = [] for i,flag in enumerate(insFlags): if flag: if flag and insNames[i] in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)']: insVary.append(insNames[i]) instDict = dict(zip(insNames,insVals)) instDict['X'] = max(instDict['X'],0.01) instDict['Y'] = max(instDict['Y'],0.01) instDict['SH/L'] = max(instDict['SH/L'],0.0001) instDict['SH/L'] = max(instDict['SH/L'],0.002) return dataType,instDict,insVary ptstr += 12*' ' print '%s'%(('Peak'+str(i+1)).center(8)),ptstr def devPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg): parmdict.update(zip(varylist,values)) return np.sqrt(weights)*getPeakProfileDerv(parmdict,xdata,varylist,bakType) def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg): parmdict.update(zip(varylist,values)) M = np.sqrt(weights)*(ydata-getPeakProfile(parmdict,xdata,varylist,bakType)) M = np.sqrt(weights)*(getPeakProfile(parmdict,xdata,varylist,bakType)-ydata) Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.) if dlg: return M Ftol = 0.0001 if controls: Ftol = controls['min dM/M'] derivType = controls['deriv type'] else: Ftol = 0.0001 derivType = 'analytic' if oneCycle: Ftol = 1.0 dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5)) try: result = so.leastsq(errPeakProfile,values,full_output=True,epsfcn=1.e-8,ftol=Ftol, args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg)) if derivType == 'analytic': result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True, args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg)) ncyc = int(result[2]['nfev']/2) else: result = so.leastsq(errPeakProfile,values,full_output=True,ftol=Ftol,epsfcn=1.e-8, args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg)) ncyc = int(result[2]['nfev']/len(varyList)) finally: dlg.Destroy() runtime = time.time()-begin chisq = np.sum(result[2]['fvec']**2) ncyc = int(result[2]['nfev']/len(varyList)) Values2Dict(parmDict, varyList, result[0]) Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to % 'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5, } global parmDict2 parmDict2 = { 'pos0':5.7,'int0':10.0,'sig0':0.5,'gam0':1.0, 'U':1.,'V':-1.,'W':0.1,'X':0.0,'Y':2.,'SH/L':0.004, 'Back0':5.,'Back1':-0.02,'Back2':.004, 'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5, } global varyList varyList = [] print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0 def test2(name,delt): if NeedTestData: TestData() varyList = [name,] xdata = np.linspace(5.6,5.8,800) hplot = plotter.add('derivatives test for '+name).gca() hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0]) y0 = getPeakProfile(parmDict2,xdata,varyList,bakType) parmDict2[name] += delt y1 = getPeakProfile(parmDict2,xdata,varyList,bakType) hplot.plot(xdata,(y1-y0)/delt,'r+') global plotter plotter = plot.PlotNotebook() test0() #    test0() test2('I(L2)/I(L1)',.0001) print "OK" plotter.StartEventLoop()
• ## trunk/GSASIIpwdGUI.py

 r345 def OnPeakFit(FitPgm,oneCycle=False): SaveState() print 'Peak Fitting:' controls = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,self.root, 'Controls')) if not controls: controls = {'deriv type':'analytic','min dM/M':0.0001,}     #fill in defaults if needed print 'Peak Fitting with '+controls['deriv type']+' derivatives:' PatternId = self.PatternId PickId = self.PickId wx.BeginBusyCursor() try: G2pwd.DoPeakFit(FitPgm,peaks,background,limits,inst,data,oneCycle) G2pwd.DoPeakFit(FitPgm,peaks,background,limits,inst,data,oneCycle,controls) finally: wx.EndBusyCursor()
• ## trunk/fsource/powsubs/psvfcj.for

 r210 SUBROUTINE PSVFCJ(DTT,TTHETA,SL,HL,SIG,GAM,PRFUNC,DPRDT,SLPART, 1  HLPART,SIGPART,GAMPART) SUBROUTINE PSVFCJ(DTT,TTHETA,SIG,GAM,SHL,PRFUNC, 1  DPRDT,SIGPART,GAMPART,SHLPART) !PURPOSE: Compute function & derivatives for Pseudovoigt profile !   [L.W. Finger, D.E. Cox & A.P. Jephcoat (1994) J. Appl. Cryst.,27,892-900.] ! coded 11/95 by B. H. Toby (NIST). revised version ! parameterized as asym1=S/L asym2=H/L ! parameterized as asym=S/L+H/L REAL*4        DTT                 !delta 2-theta in centidegrees REAL*4        TTHETA              !2-theta in centidegrees REAL*4        SL,HL               !S/L & H/L REAL*4        SIG,GAM REAL*4        SHL                 !S/L + H/L REAL*4        PRFUNC REAL*4        DPRDT REAL*4        SLPART,HLPART REAL*4        SIGPART,GAMPART REAL*4        SHLPART !INCLUDE STATEMENTS: REAL*4        DFDA REAL*4        DGDA REAL*4        DGDB REAL*4        DYDA REAL*4        DYDB REAL*4        SIN2THETA2          ! sin(2theta)**2 REAL*4        COS2THETA           ! cos(2theta) REAL*4        COSDELTA2           ! cos(Delta)**2 REAL*4        A                   ! asym1 [coff(7)] REAL*4        B                   ! asym2 [coff(8)] REAL*4        APB                 ! (A+B) REAL*4        AMB                 ! (A-B) REAL*4        APB2                ! (A+B)**2 REAL*4        TTHETAD             ! Two Theta in degrees REAL*4        SUMWDGDA            !      sum of w dGdA REAL*4        SUMWRDGDA           !      sum of w R dGdA REAL*4        SUMWDGDB            !      sum of w dGdB REAL*4        SUMWRDGDB           !      sum of w R dGdB REAL*4        SUMWGDRD2T          !      sum of w G dRd(2theta) REAL*4        SUMWGDRDSIG         !      sum of w G dRdp(n) REAL*4        SUMWGDRDGAM         !      sum of w G dRdp(n) REAL*4        SUMWGDRDA REAL*4        SUMWGDRDB REAL*4        EMIN                ! 2phi minimum REAL*4        EINFL               ! 2phi of inflection point C Asymmetry terms c A = SL            ! A = S/L in FCJ B = HL            ! B = H/L in FCJ ApB = A+B AmB = A-B A = SHL/2.            ! A = S/L in FCJ ApB = SHL ApB2 = ApB*ApB c C handle the case where there is asymmetry c IF (A .ne. 0.0 .or. B .ne. 0.0) then Einfl = Acosd(SQRT(1.0 + AmB**2)*cos2THETA) ! 2phi(infl) FCJ eq 5 (degrees) IF (A .ne. 0.0) then Einfl = Acosd(cos2THETA) ! 2phi(infl) FCJ eq 5 (degrees) tmp2 = 1.0 + ApB2 tmp = SQRT(tmp2)*cos2THETA c C Treat case where A or B is zero - Set Einfl = 2theta c if (A.eq.0.0 .or. B .eq. 0.0)Einfl = Acosd(cos2THETA) C Treat case where A is zero - Set Einfl = 2theta c if (A.eq.0.0) Einfl = Acosd(cos2THETA) if (abs(tmp) .le. 1.0) then Emin = Acosd(tmp)      ! 2phi(min) FCJ eq 4 (degrees) endif if (tmp1 .gt. 0 .and. abs(tmp) .le. 1.0) then dEmindA = -ApB*cos2THETA/SQRT(tmp1) ! N. B. symm w/r A,B dEmindA = -ApB*cos2THETA/SQRT(tmp1) ELSE dEmindA = 0.0 sumWdGdA = 0. sumWRdGdA = 0. sumWdGdB = 0. sumWRdGdB = 0. sumWGdRd2t = 0. sumWGdRdsig = 0. c if(abs(delta-emin) .gt. abs(einfl-emin))then if ( A.ge.B) then c C N.B. this is the only place where d()/dA <> d()/dB c G = 2.0*B*F*RcosDELTA dGdA = 2.0*B*RcosDELTA*(dFdA + F*tanDELTA*dDELTAdA) dGdB = dGdA + 2.0*F*RcosDELTA else G = 2.0*A*F*RcosDELTA dGdB = 2.0*A*RcosDELTA*(dFdA + F*tanDELTA*dDELTAdA) dGdA = dGdB + 2.0*F*RcosDELTA endif dGdA = 2.0*A*RcosDELTA*(dFdA + F*tanDELTA*dDELTAdA) else                                            ! delta .le. einfl .or. min(A,B) .eq. 0 G = (-1.0 + ApB*F) * RcosDELTA dGdA = RcosDELTA*(F - tanDELTA*dDELTAdA 1             + ApB*F*tanDELTA*dDELTAdA + ApB*dFdA) dGdB = dGdA endif sumWdGdA = sumWdGdA + wp(k+it) * dGdA sumWRdGdA = sumWRdGdA + wp(k+it) * R * dGdA sumWdGdB = sumWdGdB + wp(k+it) * dGdB sumWRdGdB = sumWRdGdB + wp(k+it) * R * dGdB sumWGdRd2t = sumWGdRd2t + WG * dRdT ! N.B. 1/centidegrees sumWGdRdsig = sumWGdRdsig + WG * dRdS dydA = (-(sumWRG*sumWdGdA) + 1    sumWG*(sumWRdGdA- 100.0 * todeg * sumWGdRdA)) * RsumWG2 dydB = (-(sumWRG*sumWdGdB) + 1    sumWG*(sumWRdGdB- 100.0 * todeg * sumWGdRdA)) * RsumWG2 sigpart = sumWGdRdsig * RsumWG gampart = sumWGdRdgam * RsumWG CALL PSVOIGT(DTT,SIG,GAM,R,dRdT,dRdS,dRdG) PRFUNC = R dydA = 0.002 * sign(1.0,TTHETA - DTT) dydB = dydA dydA = 0.004 * sign(1.0,TTHETA - DTT) sigpart = dRdS gampart = dRdG DPRDT = -dRdT END IF SLPART = DYDA HLPART = DYDB SHLPART = DYDA RETURN
• ## trunk/fsource/pypowder.for

 r352 C TTHETA in degrees C SPH is S/L + H/L C RETURNS FUNCTION ONLY Cf2py intent(in) NPTS Cf2py intent(in) DTT REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1) SL = SPH/2.0 FW = (2.355*SQRT(SIG)+GAM)/100.0 FMIN = 10.0*(-FW-SPH*COSD(TTHETA)) FMAX = 15.0*FW DO I=0,NPTS-1 CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SL,SL,SIG,GAM, 1    PRFUNC(I),DPRDT,SLPART,HLPART,SIGPART,GAMPART) CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH, 1    PRFUNC(I),DPRDT,SIGPART,GAMPART,SLPART) END DO RETURN END SUBROUTINE PYDPSVFCJ(NPTS,DTT,TTHETA,SIG,GAM,SPH,PRFUNC, 1  DPRDT,SIGPART,GAMPART,SLPART) C DTT in degrees C TTHETA in degrees C SPH is S/L + H/L C RETURNS FUNCTION & DERIVATIVES Cf2py intent(in) NPTS Cf2py intent(in) DTT cf2py depend(NPTS) DTT Cf2py intent(in) TTHETA Cf2py intent(in) SIG Cf2py intent(in) GAM Cf2py intent(in) SPH Cf2py intent(out) PRFUNC Cf2py depend(NPTS) PRFUNC Cf2py intent(out) DPRDT Cf2py depend(NPTS) DPRDT Cf2py intent(out) SIGPART Cf2py depend(NPTS) SIGPART Cf2py intent(out) GAMPART Cf2py depend(NPTS) GAMPART Cf2py intent(out) SLPART Cf2py depend(NPTS) SLPART REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1), 1  GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1) FW = (2.355*SQRT(SIG)+GAM)/100.0 FMIN = 10.0*(-FW-SPH*COSD(TTHETA)) FMAX = 15.0*FW DO I=0,NPTS-1 CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH, 1    PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I),SLPART(I)) DPRDT(I) = DPRDT(I)*100. END DO RETURN END
Note: See TracChangeset for help on using the changeset viewer.