Changeset 5271
 Timestamp:
 May 5, 2022 9:29:34 PM (19 months ago)
 Location:
 trunk
 Files:

 6 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIdataGUI.py
r5265 r5271 4628 4628 if self.G2plotNB: 4629 4629 self.G2plotNB.Destroy() 4630 if self.undofile :4630 if self.undofile and os.path.exists(self.undofile): 4631 4631 os.remove(self.undofile) 4632 4632 sys.exit() … … 7944 7944 # import importlib as imp 7945 7945 # imp.reload(G2pdG) 7946 # imp.reload(G2pwd) 7947 # imp.reload(G2plt) 7946 7948 # print('reloading G2pdG') 7947 7949 G2pdG.UpdatePeakGrid(G2frame,data) 
trunk/GSASIImath.py
r5260 r5271 4761 4761 gam = getTOFgamma(ins,dsp) 4762 4762 XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0] 4763 elif 'C' in Parms['Type'][0] :4763 elif 'C' in Parms['Type'][0] or 'LF' in Parms['Type'][0]: 4764 4764 for x in ['U','V','W','X','Y','Z']: 4765 4765 ins[x] = Parms.get(x,[0.0,0.0])[ind] 
trunk/GSASIIplot.py
r5268 r5271 3722 3722 G2frame.dataWindow.movePeak.Enable(len(selectedPeaks) == 1) # allow peak move from table when one peak is selected 3723 3723 for i,item in enumerate(data['peaks']): 3724 if type(item) is dict: continue 3724 3725 if i in selectedPeaks: 3725 3726 Ni = N+1 
trunk/GSASIIpwd.py
r5269 r5271 52 52 print ('pydiffax is not available for this platform') 53 53 import GSASIIfiles as G2fil 54 #import LaueFringe as LF55 54 56 55 # trig functions in degrees … … 864 863 alpPink = lambda pos,alp0,alp1: alp0+alp1*tand(pos/2.) 865 864 betPink = lambda pos,bet0,bet1: bet0+bet1*tand(pos/2.) 866 if 'T' in Inst['Type'][0]: 865 if 'LF' in Inst['Type'][0]: 866 return 3 867 elif 'T' in Inst['Type'][0]: 867 868 dsp = pos/Inst['difC'][1] 868 869 alp = alpTOF(dsp,Inst['alpha'][1]) … … 1409 1410 yb = getBackground('',parmDict,bakType,dataType,xdata,fixback)[0] 1410 1411 yc = np.zeros_like(yb) 1411 # if parmDict.get('LaueFringe',False): 1412 # if 'Lam1' in parmDict.keys(): 1413 # lam = parmDict['Lam1'] 1414 # lam2 = parmDict['Lam2'] 1415 # Ka2 = True 1416 # lamRatio = 360*(lam2lam)/(np.pi*lam) 1417 # kRatio = parmDict['I(L2)/I(L1)'] 1418 # else: 1419 # lam = parmDict['Lam'] 1420 # Ka2 = False 1421 # shol = 0 1422 # # loop over peaks 1423 # iPeak = 1 1424 # cells = parmDict['ncell'] 1425 # while True: 1426 # iPeak += 1 1427 # try: 1428 # pos = parmDict['pos'+str(iPeak)] 1429 # #tth = (posparmDict['Zero']) 1430 # intens = parmDict['int'+str(iPeak)] 1431 # damp = parmDict['damp'+str(iPeak)] 1432 # asym = parmDict['asym'+str(iPeak)] 1433 # sig = parmDict['sig'+str(iPeak)] 1434 # gam = parmDict['gam'+str(iPeak)] 1435 # fmin = 4 # for now make peaks 4 degrees wide 1436 # fmin = min(0.9*abs(xdata[1]  xdata[0]),fmin) # unless the data range is smaller 1437 # iBeg = np.searchsorted(xdata,posfmin/2) 1438 # iFin = np.searchsorted(xdata,pos+fmin/2) 1439 # if not iBeg+iFin: # skip peak below low limit 1440 # continue 1441 # elif not iBegiFin: # got peak above high limit (peaks sorted, so we can stop) 1442 # break 1443 # LF.LaueFringePeakCalc(xdata,yc,lam,pos,intens,damp,asym,sig,gam,shol,cells,fmin) 1444 # if Ka2: 1445 # pos2 = pos+lamRatio*tand(pos/2.0) # + 360/pi * Dlam/lam * tan(th) 1446 # iBeg = np.searchsorted(xdata,pos2fmin) 1447 # iFin = np.searchsorted(xdata,pos2+fmin) 1448 # if iBegiFin: 1449 # LF.LaueFringePeakCalc(xdata,yc,lam2,pos2,intens*kRatio,damp,asym,sig,gam,shol,cells,fmin) 1450 # except KeyError: #no more peaks to process 1451 # return yb+yc 1452 # elif 'C' in dataType: 1453 if 'C' in dataType: 1412 if 'LF' in dataType: 1413 if 'Lam1' in parmDict.keys(): 1414 lam = parmDict['Lam1'] 1415 lam2 = parmDict['Lam2'] 1416 Ka2 = True 1417 lamRatio = 360*(lam2lam)/(np.pi*lam) 1418 kRatio = parmDict['I(L2)/I(L1)'] 1419 else: 1420 lam = parmDict['Lam'] 1421 Ka2 = False 1422 shol = 0 1423 # loop over peaks 1424 iPeak = 1 1425 try: 1426 ncells = parmDict['ncell'] 1427 clat = parmDict['clat'] 1428 except KeyError: # no Laue info must be bkg fit 1429 print('Laue Fit: no params, assuming bkg fit') 1430 return yb 1431 while True: 1432 iPeak += 1 1433 try: 1434 #Qcen = 2 * np.pi * lam * (iPeak+1) / parmDict['clat'] 1435 l = parmDict['l'+str(iPeak)] 1436 pos = 360 * np.arcsin(0.5 * lam * l / parmDict['clat']) / np.pi 1437 parmDict['pos'+str(iPeak)] = pos 1438 #tth = (posparmDict['Zero']) 1439 intens = parmDict['int'+str(iPeak)] 1440 damp = parmDict['damp'+str(iPeak)] 1441 asym = parmDict['asym'+str(iPeak)] 1442 sig = parmDict['sig'+str(iPeak)] 1443 gam = parmDict['gam'+str(iPeak)] 1444 fmin = 8 # for now make peaks 8 degrees wide 1445 fmin = min(0.9*abs(xdata[1]  xdata[0]),fmin) # unless the data range is smaller 1446 iBeg = np.searchsorted(xdata,posfmin/2) 1447 iFin = np.searchsorted(xdata,pos+fmin/2) 1448 if not iBeg+iFin: # skip peak below low limit 1449 continue 1450 elif not iBegiFin: # got peak above high limit (peaks sorted, so we can stop) 1451 break 1452 #LF.plotme(fmin,lam,pos,intens,sig,gam,shol,ncells,clat,damp,asym) 1453 LaueFringePeakCalc(xdata,yc,lam,pos,intens,sig,gam,shol,ncells,clat,damp,asym,fmin) 1454 if Ka2: 1455 pos2 = pos+lamRatio*tand(pos/2.0) # + 360/pi * Dlam/lam * tan(th) 1456 iBeg = np.searchsorted(xdata,pos2fmin) 1457 iFin = np.searchsorted(xdata,pos2+fmin) 1458 if iBegiFin: 1459 LaueFringePeakCalc(xdata,yc,lam2,pos2,intens*kRatio,sig,gam,shol,ncells,clat,damp,asym,fmin) 1460 except KeyError: #no more peaks to process 1461 return yb+yc 1462 elif 'C' in dataType: 1454 1463 shl = max(parmDict['SH/L'],0.002) 1455 1464 Ka2 = False … … 1647 1656 ip = names.index(parm) 1648 1657 dMdv[varyList.index(name)] = dMdpk[4*int(Id)+ip] 1649 # if parmDict.get('LaueFringe',False): 1650 # for i,name in enumerate(varyList): 1651 # if not np.all(dMdv[i] == 0): continue 1652 # deltaParmDict = parmDict.copy() 1653 # delta = max(parmDict[name]/1e5,0.001) 1654 # deltaParmDict[name] += delta 1655 # intArrP = getPeakProfile(dataType,deltaParmDict,xdata,fixback,varyList,bakType) 1656 # deltaParmDict[name] = 2*delta 1657 # intArrM = getPeakProfile(dataType,deltaParmDict,xdata,fixback,varyList,bakType) 1658 # dMdv[i] = 0.5 * (intArrP  intArrM) / delta 1659 # return dMdv 1658 if 'LF' in dataType: 1659 for i,name in enumerate(varyList): 1660 if not np.all(dMdv[i] == 0): continue 1661 deltaParmDict = parmDict.copy() 1662 delta = max(parmDict[name]/1e5,0.001) 1663 deltaParmDict[name] += delta 1664 #print('num. deriv for',name,'val',deltaParmDict[name],'delta',delta) 1665 intArrP = getPeakProfile(dataType,deltaParmDict,xdata,fixback,varyList,bakType) 1666 deltaParmDict[name] = 2*delta 1667 intArrM = getPeakProfile(dataType,deltaParmDict,xdata,fixback,varyList,bakType) 1668 dMdv[i] = 0.5 * (intArrP  intArrM) / delta 1669 return dMdv 1660 1670 if 'C' in dataType: 1661 1671 shl = max(parmDict['SH/L'],0.002) … … 2112 2122 return True 2113 2123 2114 def getHeaderInfo( FitPgm,dataType):2124 def getHeaderInfo(dataType): 2115 2125 '''Provide parameter name, label name and formatting information for the 2116 2126 contents of the peak table and where used in DoPeakFit … … 2118 2128 names = ['pos','int'] 2119 2129 lnames = ['position','intensity'] 2120 if FitPgm == 'LaueFringe':2121 names += ['damp','asym','sig','gam']2122 lnames += ['satellite\ndamping','satellite\nasym','sigma','gamma']2123 fmt = ["%10. 5f","%10.2f","%10.3f","%10.3f","%10.3f","%10.3f"]2130 if 'LF' in dataType: 2131 names = ['int','sig','gam','damp','asym','l'] 2132 lnames = ['intensity','sigma','gamma','satellite\ndamping','satellite\nasym','00l'] 2133 fmt = ["%10.2f","%10.3f","%10.3f","%10.3f","%10.3f","%4d"] 2124 2134 elif 'C' in dataType: 2125 2135 names += ['sig','gam'] … … 2146 2156 2147 2157 :param str FitPgm: type of fit to perform. At present this is ignored. 2148 :param str FitPgm: type of fit to perform.2149 This should be 'LSQ' for a standard leastsquares fit2150 with pseudoVoigt peaks or 'LaueFringe' for onspecular thinfilm fitting.2151 2158 :param list Peaks: a list of peaks. Each peak entry is a list with paired values: 2152 2159 The number of pairs depends on the data type (see :func:`getHeaderInfo`). … … 2333 2340 peakDict = {} 2334 2341 peakVary = [] 2335 names,_,_ = getHeaderInfo(FitPgm,dataType) 2342 names,_,_ = getHeaderInfo(dataType) 2343 off = 0 2344 if 'LF' in dataType: off = 2 2336 2345 for i,peak in enumerate(Peaks): 2337 2346 if type(peak) is dict: … … 2340 2349 for j,name in enumerate(names): 2341 2350 parName = name+str(i) 2342 peakDict[parName] = peak[ 2*j]2343 if peak[ 2*j+1]:2351 peakDict[parName] = peak[off+2*j] 2352 if peak[off+2*j+1]: 2344 2353 peakVary.append(parName) 2345 2354 return peakDict,peakVary 2346 2355 2347 2356 def GetPeaksParms(Inst,parmDict,Peaks,varyList): 2348 names,_,_ = getHeaderInfo(FitPgm,Inst['Type'][0]) 2357 off = 0 2358 if 'LF' in Inst['Type'][0]: 2359 off = 2 2360 if 'clat' in varyList: 2361 Peaks[1]['clat'] = parmDict['clat'] 2362 names,_,_ = getHeaderInfo(Inst['Type'][0]) 2349 2363 for i,peak in enumerate(Peaks): 2350 2364 if type(peak) is dict: 2351 2365 continue 2366 for j in range(len(names)): 2367 parName = names[j]+str(i) 2368 if parName in varyList or not peakInstPrmMode: 2369 peak[2*j+off] = parmDict[parName] 2352 2370 pos = parmDict['pos'+str(i)] 2371 if 'LF' in Inst['Type'][0]: peak[0] = pos 2353 2372 if 'difC' in Inst: 2354 2373 dsp = pos/Inst['difC'][1] 2355 2374 for j in range(len(names)): 2356 parName = names[j]+str(i) 2357 if parName in varyList or not peakInstPrmMode: 2358 peak[2*j] = parmDict[parName] 2359 elif 'alp' in parName: 2375 if 'alp' in parName: 2360 2376 if 'T' in Inst['Type'][0]: 2361 2377 peak[2*j] = G2mth.getTOFalpha(parmDict,dsp) … … 2381 2397 2382 2398 def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW): 2399 if 'clat' in varyList: 2400 print('c = {:.6f} esd {:.6f}'.format(parmDict['clat'],sigDict['clat'])) 2383 2401 print ('Peak coefficients:') 2384 names,fmt,_ = getHeaderInfo( FitPgm,dataType)2402 names,fmt,_ = getHeaderInfo(dataType) 2385 2403 head = 13*' ' 2386 2404 for name in names: … … 2447 2465 xBeg = np.searchsorted(x,Limits[0]) 2448 2466 xFin = np.searchsorted(x,Limits[1])+1 2467 # find out what is varied 2449 2468 bakType,bakDict,bakVary = SetBackgroundParms(Background) 2450 2469 dataType,insDict,insVary = SetInstParms(Inst) … … 2460 2479 else: 2461 2480 varyList = bakVary+insVary+peakVary 2481 if 'LF' in Inst['Type'][0] and Peaks: 2482 if Peaks[1].get('clatref'): varyList += ['clat'] 2462 2483 fullvaryList = varyList[:] 2463 2484 if not peakInstPrmMode: … … 2467 2488 raise Exception('Instrumental profile terms cannot be varied '+ 2468 2489 'after setPeakInstPrmMode(False) is used') 2469 parmDict['LaueFringe'] = False 2470 # if FitPgm == 'LaueFringe': 2471 # #====================================================================== 2472 # print('Debug: reload LaueFringe') # TODO: remove me 2473 # import imp 2474 # imp.reload(LF) 2475 # # TODO: remove ^^^^ 2476 # #====================================================================== 2477 # for v in ('U','V','W','X','Y','Z','alpha','alpha0','alpha1', 2478 # 'beta0','beta1','betaq','sig0','sig1','sig2','sigq',): 2479 # if v in varyList: 2480 # raise Exception('Instrumental profile terms cannot be varied '+ 2481 # 'in Laue Fringe fits') 2482 # parmDict['LaueFringe'] = True 2490 if 'LF' in Inst['Type'][0]: 2491 warn = [] 2492 for v in ('U','V','W','X','Y','Z','alpha','alpha0','alpha1', 2493 'beta0','beta1','betaq','sig0','sig1','sig2','sigq',): 2494 if v in varyList: 2495 warn.append(v) 2496 del varyList[varyList.index(v)] 2497 if warn: 2498 print('Instrumental profile terms cannot be varied '+ 2499 'in Laue Fringe fits:',warn) 2483 2500 2484 2501 while not noFit: … … 2487 2504 Rvals = {} 2488 2505 badVary = [] 2489 result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True, 2506 try: 2507 result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True, 2490 2508 args=(x[xBeg:xFin],y[xBeg:xFin],fixback[xBeg:xFin],wtFactor*w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg)) 2509 except Exception as msg: 2510 if GSASIIpath.GetConfigValue('debug'): 2511 print('peak fit failure\n',msg) 2512 import traceback 2513 print (traceback.format_exc()) 2514 else: 2515 print('peak fit failure') 2516 return 2491 2517 ncyc = int(result[2]['nfev']/2) 2492 2518 runtime = time.time()begin … … 2522 2548 yd[xBeg:xFin] = y[xBeg:xFin]yc[xBeg:xFin] 2523 2549 if noFit: 2550 GetPeaksParms(Inst,parmDict,Peaks,varyList) 2524 2551 return 2525 2552 sigDict = dict(zip(varyList,sig)) … … 5030 5057 newRefs = np.array(newRefs) 5031 5058 return True,newRefs 5032 5059 5060 #===Laue Fringe code =================================================================== 5061 import NIST_profile as FP 5062 5063 class profileObj(FP.FP_profile): 5064 def conv_Lauefringe(self): 5065 """Compute the FT of the Laue Fringe function""" 5066 5067 me=self.get_function_name() #the name of this convolver,as a string 5068 wave = self.param_dicts['conv_global']['dominant_wavelength']*1.e10 # in A 5069 pos = np.rad2deg(self.param_dicts["conv_global"]["twotheta0"]) # peak position as 2theta in deg 5070 posQ = np.pi * 4 * np.sin(self.param_dicts["conv_global"]["twotheta0"]/2) / wave # peak position as Q 5071 ttwid = self.twotheta_window_fullwidth_deg 5072 ncell = self.param_dicts[me]['Ncells'] 5073 co2 = self.param_dicts[me]['clat'] / 2. 5074 damp = self.param_dicts[me]['damp'] 5075 asym = self.param_dicts[me]['asym'] 5076 ttlist = np.linspace(posttwid/2,pos+ttwid/2,len(self._epsb2)) 5077 Qs = np.pi * 4 * np.sin(np.deg2rad(ttlist/2)) / wave 5078 w = np.exp(10**((dampasym) * (Qs  posQ)**2)) 5079 w2 = np.exp(10**((damp+asym) * (Qs  posQ)**2)) 5080 w[len(w)//2:] = w2[len(w)//2:] 5081 weqdiv = w * np.sin(Qs * ncell * co2)**2 / (np.sin(Qs * co2)**2) 5082 conv = FP.best_rfft(weqdiv) 5083 conv[1::2] *= 1 #flip center 5084 return conv 5085 5086 def conv_Lorentzian(self): 5087 """Compute the FT of a Lorentz function where gamma is the FWHM""" 5088 ttwid = self.twotheta_window_fullwidth_deg 5089 me=self.get_function_name() #the name of this convolver,as a string 5090 g2gam = self.param_dicts[me]['g2gam'] # gsasii gamma in centidegrees 5091 gamma = g2gam/100 # deg 5092 ttlist = np.linspace(ttwid/2,ttwid/2,len(self._epsb2)) 5093 eqdiv = (0.5 * gamma / np.pi) / (gamma**2/4. + ttlist**2) 5094 conv = FP.best_rfft(eqdiv) 5095 conv[1::2] *= 1 #flip center 5096 return conv 5097 5098 def conv_Gaussian(self): 5099 """Compute the FT of a Gaussian where sigma**2 is the variance""" 5100 ttwid = self.twotheta_window_fullwidth_deg 5101 me=self.get_function_name() #the name of this convolver,as a string 5102 g2sig2 = self.param_dicts[me]['g2sig2'] # gsasii sigma**2 in centidegr**2 5103 sigma = math.sqrt(g2sig2)/100. 5104 ttlist = np.linspace(ttwid/2,ttwid/2,len(self._epsb2)) 5105 eqdiv = np.exp(0.5*ttlist**2/sigma**2) / math.sqrt(2*np.pi*sigma**2) 5106 conv = FP.best_rfft(eqdiv) 5107 conv[1::2] *= 1 #flip center 5108 return conv 5109 5110 def LaueFringePeakCalc(ttArr,intArr,lam,peakpos,intens,sigma2,gamma,shol,ncells,clat,damp,asym,calcwid,plot=False): 5111 '''Compute the peakshape for a Laue Fringe peak convoluted with a Gaussian, Lorentzian & 5112 an axial divergence asymmetry correction. 5113 5114 :param np.array ttArr: Array of twotheta values (in degrees) 5115 :param np.array intArr: Array of intensity values (peaks are added to this) 5116 :param float lam: wavelength in Angstrom 5117 :param float peakpos: peak position in twotheta (deg.) 5118 :param float intens: intensity factor for peak 5119 :param float sigma2: Gaussian variance (in centidegrees**2) ** 5120 :param float gamma: Lorenzian FWHM (in centidegrees) ** 5121 :param float shol: FCJ (S + H)/L where S=samplehalf height, H=slit halfheight, L=radius ** 5122 :param float ncells: number of unit cells in specular direction ** 5123 :param float clat: c lattice parameter ** 5124 :param float damp: 5125 :param float asym: 5126 :param float calcwid: twotheta (deg.) width for cutoff of peak computation. 5127 Defaults to 5 5128 :param bool plot: for debugging, shows contributions to peak 5129 5130 ** If term is <= zero, item is removed from convolution 5131 ''' 5132 def LaueFringePeakPlot(ttArr,intArr): 5133 import matplotlib.pyplot as plt 5134 refColors = ['xkcd:blue','xkcd:red','xkcd:green','xkcd:cyan','xkcd:magenta','xkcd:black', 5135 'xkcd:pink','xkcd:brown','xkcd:teal','xkcd:orange','xkcd:grey','xkcd:violet',] 5136 fig, ax = plt.subplots() 5137 ax.set(title='Peak convolution functions @ 2theta={:.3f}'.format(peakpos), 5138 xlabel=r'$\Delta 2\theta, deg$', 5139 ylabel=r'Intensity (arbitrary)') 5140 ax.set_yscale("log",nonpositive='mask') 5141 ttmin = ttmax = 0 5142 for i,conv in enumerate(convList): 5143 f = NISTpk.convolver_funcs[conv]() 5144 if f is None: continue 5145 FFT = FP.best_irfft(f) 5146 if f[1].real > 0: FFT = np.roll(FFT,int(len(FFT)/2.)) 5147 FFT /= FFT.max() 5148 if i == 0: 5149 tt = np.linspace(NISTpk.twotheta_window_fullwidth_deg/2, 5150 NISTpk.twotheta_window_fullwidth_deg/2,len(FFT)) 5151 ttmin = min(ttmin,tt[np.argmax(FFT>.005)]) 5152 ttmax = max(ttmax,tt[::1][np.argmax(FFT[::1]>.005)]) 5153 color = refColors[i%len(refColors)] 5154 ax.plot(tt,FFT,color,label=conv[5:]) 5155 color = refColors[(i+1)%len(refColors)] 5156 ax.plot(ttArrpeakpos,intArr/max(intArr),color,label='Convolution') 5157 ax.set_xlim((ttmin,ttmax)) 5158 ax.legend(loc='best') 5159 plt.show() 5160 # hardcoded constants 5161 diffRadius = 220 # diffractometer radius in mm; needed for axial divergence, etc, but should not matter 5162 axial_factor = 1.5 # fudge factor to bring sh/l broadening to ~ agree with FPA 5163 equatorial_divergence_deg = 0.5 # not sure exactly what this impacts 5164 NISTparms = { 5165 "": { 5166 'equatorial_divergence_deg' : equatorial_divergence_deg, 5167 'dominant_wavelength' : 1.e10 * lam, 5168 'diffractometer_radius' : 1e3* diffRadius, # diffractometer radius in m 5169 'oversampling' : 8, 5170 }, 5171 "emission": { 5172 'emiss_wavelengths' : 1.e10 * np.array([lam]), 5173 'emiss_intensities' : np.array([1.]), 5174 'emiss_gauss_widths' : 1.e10 * 1.e3 * np.array([0.001]), 5175 'emiss_lor_widths' : 1.e10 * 1.e3 * np.array([0.001]), 5176 'crystallite_size_gauss' : 1.e9 * 1e6, 5177 'crystallite_size_lor' : 1.e9 * 1e6 5178 }, 5179 "axial": { 5180 'axDiv':"full", 5181 'slit_length_source' : 1e3 * diffRadius * shol * axial_factor, 5182 'slit_length_target' : 1e3 * diffRadius * shol * 1.00001 * axial_factor, # != 'slit_length_source' 5183 'length_sample' : 1e3 * diffRadius * shol * axial_factor, 5184 'n_integral_points' : 10, 5185 'angI_deg' : 2.5, 5186 'angD_deg': 2.5, 5187 }, 5188 'Gaussian': {'g2sig2': sigma2}, 5189 'Lorentzian': {'g2gam': gamma}, 5190 'Lauefringe': {'Ncells': ncells, 'clat':clat, 'damp': damp, 'asym': asym}, 5191 } 5192 NISTpk=profileObj(anglemode="twotheta", 5193 output_gaussian_smoother_bins_sigma=1.0, 5194 oversampling=NISTparms.get('oversampling',10)) 5195 NISTpk.debug_cache=False 5196 for key in NISTparms: #set parameters for each convolver 5197 if key: 5198 NISTpk.set_parameters(convolver=key,**NISTparms[key]) 5199 else: 5200 NISTpk.set_parameters(**NISTparms[key]) 5201 # find closest point to peak location (which may be outside limits of the array) 5202 center_bin_idx=min(ttArr.searchsorted(peakpos),len(ttArr)1) 5203 step = (ttArr[1]ttArr[0])/(len(ttArr)1) 5204 NISTpk.set_optimized_window(twotheta_exact_bin_spacing_deg=step, 5205 twotheta_window_center_deg=ttArr[center_bin_idx], 5206 twotheta_approx_window_fullwidth_deg=calcwid, 5207 ) 5208 NISTpk.set_parameters(twotheta0_deg=peakpos) 5209 convList = ['conv_emission'] 5210 if ncells: convList += ['conv_Lauefringe'] 5211 if sigma2 > 0: convList += ['conv_Gaussian'] 5212 if gamma > 0: convList += ['conv_Lorentzian'] 5213 if shol > 0: convList += ['conv_axial'] 5214 5215 # global deriv 5216 # if deriv: 5217 # peakObj = NISTpk.compute_line_profile(convolver_names=convList,compute_derivative=True) 5218 # else: 5219 # peakObj = NISTpk.compute_line_profile(convolver_names=convList) 5220 peakObj = NISTpk.compute_line_profile(convolver_names=convList) 5221 5222 pkPts = len(peakObj.peak) 5223 pkMax = peakObj.peak.max() 5224 startInd = center_bin_idx(pkPts//2) 5225 istart = None 5226 pstart = None 5227 iend = None 5228 pend = None 5229 # adjust data range if peak calc begins below data range or ends above data range 5230 # but range of peak calc should not extend past both ends of ttArr 5231 if startInd < 0: 5232 iend = startInd+pkPts 5233 pstart = startInd 5234 elif startInd > len(intArr): 5235 return 5236 elif startInd+pkPts >= len(intArr): 5237 offset = pkPts  len( intArr[startInd:] ) 5238 istart = startInd 5239 iend = startInd+pkPtsoffset 5240 pend = offset 5241 else: 5242 istart = startInd 5243 iend = startInd+pkPts 5244 intArr[istart:iend] += intens * peakObj.peak[pstart:pend]/pkMax 5245 if plot: 5246 LaueFringePeakPlot(ttArr[istart:iend], (intens * peakObj.peak[pstart:pend]/pkMax)) 5247 5248 def LaueSatellite(peakpos,wave,c,ncell,j=[4,3,2,1,0,1,2,3,4]): 5249 '''Returns the locations of the Laue satellite positions relative 5250 to the peak position 5251 5252 :param float peakpos: the peak position in degrees 2theta 5253 :param float ncell: Laue fringe parameter, number of unit cells in layer 5254 :param list j: the satellite order, where j=1 is the first satellite 5255 on the lower 2theta side and j=1 is the first satellite on the high 5256 2theta side. j=0 gives the peak position 5257 ''' 5258 Qpos = 4 * np.pi * np.sin(peakpos * np.pi / 360) / wave 5259 dQvals = (2 * np.array(j) + np.sign(j)) * np.pi / (c * ncell) 5260 return np.arcsin((Qpos+dQvals)*wave/(4*np.pi)) * (360 / np.pi) 5261 5033 5262 #### testing data 5034 5263 NeedTestData = True 
trunk/GSASIIpwdGUI.py
r5265 r5271 834 834 if not G2frame.GSASprojectfile: #force a save of the gpx file so SaveState can write in the same directory 835 835 G2frame.OnFileSaveas(event) 836 FitPgm = 'LSQ' 837 # if data.get('FitPgm',0) == 1: 838 # FitPgm = 'LaueFringe' 839 wx.CallAfter(OnPeakFit,FitPgm) 836 wx.CallAfter(OnPeakFit) 840 837 841 838 def OnOneCycle(event): … … 844 841 reflGrid.HideCellEditControl() 845 842 reflGrid.DisableCellEditControl() 846 FitPgm = 'LSQ' 847 # if data.get('FitPgm',0) == 1: 848 # FitPgm = 'LaueFringe' 849 wx.CallAfter(OnPeakFit,FitPgm,oneCycle=True) 843 wx.CallAfter(OnPeakFit,oneCycle=True) 850 844 851 845 def OnSeqPeakFit(event): … … 876 870 print ('Peak Fitting with '+controls['deriv type']+' derivatives:') 877 871 oneCycle = False 878 FitPgm = 'LSQ'879 # if data.get('FitPgm',0) == 1:880 # FitPgm = 'LaueFringe'881 872 prevVaryList = [] 882 873 peaks = None … … 903 894 data = Pattern[1] 904 895 fixback = GetFileBackground(G2frame,data,background,scale=False) 905 peaks['sigDict'],result,sig,Rvals,varyList,parmDict,fullvaryList,badVary = G2pwd.DoPeakFit( FitPgm,peaks['peaks'],896 peaks['sigDict'],result,sig,Rvals,varyList,parmDict,fullvaryList,badVary = G2pwd.DoPeakFit(None,peaks['peaks'], 906 897 background,limits,inst,inst2,data,fixback,prevVaryList,oneCycle,controls) #needs wtFactor after controls? 907 898 if len(result[0]) != len(fullvaryList): … … 933 924 G2plt.PlotPatterns(G2frame,plotType='PWDR') 934 925 935 def OnPeakFit( FitPgm,oneCycle=False,noFit=False):926 def OnPeakFit(oneCycle=False,noFit=False): 936 927 'Do peak fitting by least squares' 937 928 SaveState() … … 939 930 if not controls: 940 931 controls = {'deriv type':'analytic','min dM/M':0.001,} #fill in defaults if needed 941 print ('Peak Fitting with '+controls['deriv type']+' derivatives:')932 #print ('Peak Fitting with '+controls['deriv type']+' derivatives:') 942 933 PatternId = G2frame.PatternId 943 934 peaks = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,PatternId, 'Peak List')) … … 952 943 wtFactor = Pattern[0]['wtFactor'] 953 944 bxye = GetFileBackground(G2frame,data,background,scale=False) 954 overallInfo = [] 955 # peaks['LaueFringe'] = peaks.get('LaueFringe',{}) 956 # peaks['LaueFringe']['satellites'] = [] 957 # import LaueFringe as LF 958 # lines = peaks['LaueFringe'].get('Show') 959 # if 'Laue' in FitPgm: 960 # overallInfo = [{'ncell':peaks['LaueFringe']['ncell']}] # add overall info 961 # if lines: 962 # for i in peaks['peaks']: 963 # pks = list(range(lines,0)) + list(range(1,lines+1)) 964 # peaks['LaueFringe']['satellites'].extend( 965 # LF.LaueSatellite(i[0],peaks['LaueFringe']['ncell'],pks)) 966 # #====================================================================== 967 # print('Debug: reload G2pwd') # TODO: remove me 968 # import imp 969 # imp.reload(G2pwd) 970 # imp.reload(G2plt) 971 # # TODO: remove ^^^^ 972 #====================================================================== 945 overallInfo = {} 946 peaks['LaueFringe'] = peaks.get('LaueFringe',{}) 947 peaks['LaueFringe']['satellites'] = [] 948 lines = peaks['LaueFringe'].get('Show') 949 if 'LF' in inst['Type'][0]: 950 wave = G2mth.getWave(inst) 951 overallInfo = {'ncell':peaks['LaueFringe']['ncell'], 952 'clat': peaks['LaueFringe']['clat'], 953 'clatref': peaks['LaueFringe']['clatref'] 954 } # add overall info 955 if lines: 956 for i in peaks['peaks']: 957 pks = list(range(lines,0)) + list(range(1,lines+1)) 958 peaks['LaueFringe']['satellites'].extend( 959 G2pwd.LaueSatellite(i[0],wave,peaks['LaueFringe']['clat'],peaks['LaueFringe']['ncell'],pks)) 960 peaksplus = peaks['peaks'] + [overallInfo] 973 961 if noFit: 974 results = G2pwd.DoPeakFit( FitPgm,peaks['peaks']+overallInfo,background,limits,inst,inst2,data,bxye,[],oneCycle,controls,wtFactor,noFit=True)962 results = G2pwd.DoPeakFit(None,peaksplus,background,limits,inst,inst2,data,bxye,[],oneCycle,controls,wtFactor,noFit=True) 975 963 G2plt.PlotPatterns(G2frame,plotType='PWDR') 976 964 return 977 dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0, 978 parent=G2frame, 979 style = wx.PD_ELAPSED_TIMEwx.PD_AUTO_HIDEwx.PD_REMAINING_TIMEwx.PD_CAN_ABORT) 980 results = G2pwd.DoPeakFit(FitPgm,peaks['peaks']+overallInfo,background,limits,inst,inst2,data,bxye,[],oneCycle,controls,wtFactor,dlg) 965 try: 966 dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0,parent=G2frame, 967 style = wx.PD_ELAPSED_TIMEwx.PD_AUTO_HIDEwx.PD_REMAINING_TIMEwx.PD_CAN_ABORT) 968 results = G2pwd.DoPeakFit(None,peaksplus,background,limits,inst,inst2,data,bxye,[],oneCycle,controls,wtFactor,dlg) 969 finally: 970 dlg.Destroy() 971 if results is None: 972 return 981 973 peaks['sigDict'] = results[0] 974 if 'LF' in inst['Type'][0] and 'clat' in overallInfo: 975 peaks['LaueFringe']['clat'] = overallInfo['clat'] 982 976 text = 'Peak fit: Rwp=%.2f%% Nobs= %d Nparm= %d Npeaks= %d'%(results[3]['Rwp'],results[1][2]['fjac'].shape[1],len(results[0]),len(peaks['peaks'])) 983 977 newpeaks = copy.copy(peaks) … … 1134 1128 wx.CallAfter(UpdatePeakGrid,G2frame,data) 1135 1129 1136 def updateMe(*args):1137 'Redraw the peak listings after the mode changes'1138 wx.CallAfter(UpdatePeakGrid,G2frame,data)1139 wx.CallLater(100,RefreshPeakGrid,None)1140 1141 1130 def RefreshPeakGrid(event): 1142 1131 'recompute & plot the peaks any time a value in the table is edited' 1143 FitPgm = 'LSQ'1144 # if data.get('FitPgm',0) == 1:1145 # FitPgm = 'LaueFringe' 1146 OnPeakFit( FitPgm,noFit=True)1132 if 'LF' in Inst['Type'][0]: 1133 for i in range(len(data['LFpeaks'])): 1134 data['peaks'][i][2:] = data['LFpeaks'][i] 1135 OnPeakFit(noFit=True) 1147 1136 1148 1137 def OnSetPeakWidMode(event): … … 1155 1144 #### beginning of UpdatePeakGrid 1156 1145 #====================================================================== 1157 data['FitPgm'] = data.get('FitPgm',0)1158 1146 G2frame.GetStatusBar().SetStatusText('Global refine: select refine column & press Y or N',1) 1159 1147 G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.PeakMenu) … … 1171 1159 G2frame.Bind(wx.EVT_MENU, OnResetSigGam, id=G2G.wxID_RESETSIGGAM) 1172 1160 G2frame.Bind(wx.EVT_MENU, OnSetPeakWidMode, id=G2G.wxID_SETUNVARIEDWIDTHS) 1161 OnSetPeakWidMode(None) 1173 1162 if data['peaks']: 1174 1163 G2frame.dataWindow.AutoSearch.Enable(False) … … 1190 1179 colLabels = [] 1191 1180 Types = [] 1192 for _,f,l in zip(*G2pwd.getHeaderInfo( 1193 ('LSQ','LaueFringe')[data.get('FitPgm',0)],Inst['Type'][0])): 1181 for _,f,l in zip(*G2pwd.getHeaderInfo(Inst['Type'][0])): 1194 1182 colLabels.append(l) 1195 colLabels.append('refine') 1196 Types.append(wg.GRID_VALUE_FLOAT + 1197 f.replace('%',':').replace('f','').replace('.',',')) 1198 Types.append(wg.GRID_VALUE_BOOL) 1183 if "d" in f: 1184 Types.append(wg.GRID_VALUE_NUMBER) 1185 else: 1186 Types.append(wg.GRID_VALUE_FLOAT + f.replace('%',':').replace('f','').replace('.',',')) 1187 colLabels.append('refine') 1188 Types.append(wg.GRID_VALUE_BOOL) 1199 1189 T = [] 1200 1190 for peak in data['peaks']: … … 1207 1197 for key in T: X.append(D[key]) 1208 1198 data['peaks'] = X 1209 # pad peak table, if needed 1210 for j in range(len(data['peaks'])): 1211 for i in range(len(data['peaks'][j]),len(colLabels)): 1212 if colLabels[i] == 'refine': 1213 data['peaks'][j].append(False) 1214 else: 1215 data['peaks'][j].append(1.0) 1199 1216 1200 G2frame.dataWindow.ClearData() 1217 1201 mainSizer = wx.BoxSizer(wx.VERTICAL) 1218 1202 G2frame.GPXtree.SetItemPyData(G2frame.PickId,data) 1219 G2frame.PeakTable = G2G.Table(data['peaks'],rowLabels=rowLabels,colLabels=colLabels,types=Types) 1203 if 'LF' in Inst['Type'][0]: 1204 data['LFpeaks'] = [] 1205 for i in range(len(data['peaks'])): 1206 if len(data['peaks'][i]) == 8: # need to extend the entry 1207 if 'Lam' in Inst: 1208 lam = Inst['Lam'][0] 1209 else: 1210 lam = (Inst['Lam1'][0] + 1211 Inst['I(L2)/I(L1)'][0] * Inst['Lam2'][0]) / ( 1212 1 + Inst['I(L2)/I(L1)'][0]) 1213 if i == 0: 1214 pos = data['peaks'][i][0] 1215 data['LaueFringe']['clat'] = 0.5 * lam / np.sin(pos*np.pi/360) 1216 l = 1 1217 else: 1218 l = data['LaueFringe']['clat'] / (0.5 * lam / np.sin(data['peaks'][i][0]*np.pi/360)) 1219 l = int(l + 0.5) 1220 data['peaks'][i][0] = 360 * np.arcsin(0.5 * lam / (data['LaueFringe']['clat']/l)) / np.pi 1221 data['peaks'][i] += [0., False, 0., False, l, None] 1222 data['LFpeaks'].append(data['peaks'][i][2:]) 1223 G2frame.PeakTable = G2G.Table(data['LFpeaks'],rowLabels=rowLabels,colLabels=colLabels,types=Types) 1224 else: 1225 G2frame.PeakTable = G2G.Table(data['peaks'],rowLabels=rowLabels,colLabels=colLabels,types=Types) 1220 1226 #G2frame.SetLabel(G2frame.GetLabel().split('')[0]+'  '+'Peak List') 1221 1227 G2frame.dataWindow.currentGrids = [] … … 1234 1240 mainSizer.Add(topSizer,0,wx.EXPAND) 1235 1241 G2G.HorizontalLine(mainSizer,G2frame.dataWindow) 1236 # if 'C' in Inst['Type'][0]:1237 # topSizer = wx.BoxSizer(wx.HORIZONTAL)1238 # topSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Fitting mode: '),0,WACV)1239 # pkType = G2G.G2ChoiceButton(G2frame.dataWindow,('Powder diffraction','Laue fringes'),1240 # indLoc=data,indKey='FitPgm',1241 # onChoice=updateMe)1242 # topSizer.Add(pkType,0,WACV)1243 # if data['FitPgm']:1244 # topSizer.Add((15,1))1245 # if 'LaueFringe' not in data:1246 # data['LaueFringe'] = {'ncell':20}1247 # elif 'ncell' not in data['LaueFringe']:1248 # data['LaueFringe']['ncell'] = 201249 # siz = G2G.G2SpinWidget(G2frame.dataWindow,data['LaueFringe'] ,'ncell',1250 # 'Laue ncell',1251 # onChange=RefreshPeakGrid,onChangeArgs=[None])1252 # topSizer.Add(siz,0,WACV)1253 # data['LaueFringe']['Show'] = data['LaueFringe'].get('Show',0)1254 # #ch = G2G.G2CheckBox(G2frame.dataWindow,' Show',data['LaueFringe'],'Show',1255 # # OnChange=RefreshPeakGrid)1256 #topSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Show '),0,WACV)1257 #ch = G2G.EnumSelector(G2frame.dataWindow,data['LaueFringe'],'Show',1258 #['None','1','2','3','4','5','6'],list(range(7)),1259 #OnChange=RefreshPeakGrid)1260 #topSizer.Add(ch,0,WACV)1261 #topSizer.Add(wx.StaticText(G2frame.dataWindow,label=' satellites'),0,WACV)1262 #mainSizer.Add(topSizer)1242 if 'LF' in Inst['Type'][0]: 1243 mainSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Laue Fringe fitting')) 1244 topSizer = wx.BoxSizer(wx.HORIZONTAL) 1245 topSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Overall parms: c='),0,WACV) 1246 data['LaueFringe'] = data.get('LaueFringe',{}) 1247 data['LaueFringe']['ncell'] = data['LaueFringe'].get('ncell',20) 1248 data['LaueFringe']['clat'] = data['LaueFringe'].get('clat',9.0) 1249 data['LaueFringe']['clatref'] = data['LaueFringe'].get('clatref',False) 1250 data['LaueFringe']['Show'] = data['LaueFringe'].get('Show',0) 1251 cVal = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['LaueFringe'],'clat', 1252 typeHint=float,nDig=(10,4), 1253 OnLeave=lambda *arg,**kw:RefreshPeakGrid(None)) 1254 topSizer.Add(cVal,0,WACV) 1255 cRef = G2G.G2CheckBox(G2frame.dataWindow,'ref',data['LaueFringe'],'clatref') 1256 topSizer.Add(cRef,0,WACV) 1257 topSizer.Add((15,1)) 1258 siz = G2G.G2SpinWidget(G2frame.dataWindow,data['LaueFringe'] ,'ncell', 1259 'Laue ncell', 1260 onChange=RefreshPeakGrid,onChangeArgs=[None]) 1261 topSizer.Add(siz,0,WACV) 1262 topSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Show '),0,WACV) 1263 ch = G2G.EnumSelector(G2frame.dataWindow,data['LaueFringe'],'Show', 1264 ['None','1','2','3','4','5','6'],list(range(7)), 1265 OnChange=RefreshPeakGrid) 1266 topSizer.Add(ch,0,WACV) 1267 topSizer.Add(wx.StaticText(G2frame.dataWindow,label=' satellites'),0,WACV) 1268 mainSizer.Add(topSizer) 1263 1269 mainSizer.Add(reflGrid,0,wx.ALL,10) 1264 1270 G2frame.dataWindow.SetSizer(mainSizer) 
trunk/imports/G2pwd_csv.py
r5268 r5271 80 80 try: 81 81 positions[j] = int(v.strip().split('=')[1]) 82 if j == 1: positions[2] = 1 82 83 except: 83 84 print('Error parsing: "'+S+'"')
Note: See TracChangeset
for help on using the changeset viewer.