Changeset 345 for trunk/GSASIIpwd.py
 Timestamp:
 Aug 18, 2011 8:54:16 PM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIpwd.py
r343 r345 410 410 np.seterr(divide='ignore') 411 411 412 # Voigt function = convolution (norm,cauchy) 413 414 class voigt_gen(st.rv_continuous): 415 ''' 416 Voigt function 417 Parameters 418  419 ''' 420 def _pdf(self,x,s,g): 421 z = np.empty([2,len(x)]) 422 z[0] = st.norm.pdf(x,scale=s) 423 z[1] = st.cauchy.pdf(x,scale=g) 424 Z = fft(z) 425 return fftshift(ifft(Z.prod(axis=0))).real 426 # def _cdf(self, x): 427 # def _ppf(self, q): 428 # def _sf(self, x): 429 # def _isf(self, q): 430 # def _stats(self): 431 # def _entropy(self): 432 433 voigt = voigt_gen(name='voigt',shapes='ts,g') 434 435 # FingerCox_Jephcoat D(2phi,2th) function for S/L = H/L 412 # Normal distribution 413 414 # loc = mu, scale = std 415 _norm_pdf_C = 1./math.sqrt(2*math.pi) 416 class norm_gen(st.rv_continuous): 417 418 def pdf(self,x,*args,**kwds): 419 loc,scale=kwds['loc'],kwds['scale'] 420 x = (xloc)/scale 421 return np.exp(x**2/2.0) * _norm_pdf_C / scale 422 423 norm = norm_gen(name='norm',longname='A normal',extradoc=""" 424 425 Normal distribution 426 427 The location (loc) keyword specifies the mean. 428 The scale (scale) keyword specifies the standard deviation. 429 430 normal.pdf(x) = exp(x**2/2)/sqrt(2*pi) 431 """) 432 433 ## Cauchy 434 435 # median = loc 436 437 class cauchy_gen(st.rv_continuous): 438 439 def pdf(self,x,*args,**kwds): 440 loc,scale=kwds['loc'],kwds['scale'] 441 x = (xloc)/scale 442 return 1.0/np.pi/(1.0+x*x) / scale 443 444 cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc=""" 445 446 Cauchy distribution 447 448 cauchy.pdf(x) = 1/(pi*(1+x**2)) 449 450 This is the t distribution with one degree of freedom. 451 """) 452 453 454 #GSASII peak fitting routine: Finger, Cox & Jephcoat model 455 436 456 437 457 class fcjde_gen(st.rv_continuous): … … 457 477 def _pdf(self,x,t,s,dx): 458 478 T = dx*x+t 459 ax = npcosd(T)**2 479 ax2 = abs(npcosd(T)) 480 ax = ax2**2 460 481 bx = npcosd(t)**2 461 482 bx = np.where(ax>bx,bx,ax) 462 fx = np.where(ax>bx,(np.sqrt(bx/(axbx))1./s)/ np.sqrt(ax),0.0)483 fx = np.where(ax>bx,(np.sqrt(bx/(axbx))1./s)/ax2,0.0) 463 484 fx = np.where(fx > 0.,fx,0.0) 464 485 return fx 465 # def _cdf(self, x): 466 # def _ppf(self, q): 467 # def _sf(self, x): 468 # def _isf(self, q): 469 # def _stats(self): 470 # def _entropy(self): 486 487 def pdf(self,x,*args,**kwds): 488 loc=kwds['loc'] 489 return self._pdf(xloc,*args) 471 490 472 491 fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx') … … 485 504 return yb 486 505 487 def calcPeakFFT(x,fxns,widths,pos,args): 488 dx = x[1]x[0] 489 z = np.empty([len(fxns),len(x)]) 490 for i,(fxn,width,arg) in enumerate(zip(fxns,widths,args)): 491 z[i] = fxn.pdf(x,loc=pos,scale=width,*arg)*dx 492 Z = fft(z) 493 if len(fxns)%2: 494 return ifft(Z.prod(axis=0)).real 495 else: 496 return fftshift(ifft(Z.prod(axis=0))).real 497 498 def getFCJVoigt(pos,intens,sig,gam,shl,xdata): 499 506 def getWidths(pos,sig,gam,shl): 507 widths = [np.sqrt(sig)/100.,gam/200.] 508 fwhm = 2.355*widths[0]+2.*widths[1] 509 fmin = 10.*(fwhm+shl*abs(npcosd(pos))) 510 fmax = 15.0*fwhm 511 if pos > 90: 512 fmin,fmax = [fmax,fmin] 513 return widths,fmin,fmax 514 515 def getFCJVoigt(pos,intens,sig,gam,shl,xdata): 500 516 DX = xdata[1]xdata[0] 501 517 widths,fmin,fmax = getWidths(pos,sig,gam,shl) 502 x = np.linspace(posfmin,pos+fmin,1024) 503 if pos > 90: 504 fmin,fmax = [fmax,fmin] 518 x = np.linspace(posfmin,pos+fmin,256) 505 519 dx = x[1]x[0] 506 arg = [pos,shl/10.,dx,] 507 Df = fcjde.pdf(x,*arg,loc=pos,scale=1.0) 508 if len(np.nonzero(Df)[0])>5: 509 fxns = [st.norm,st.cauchy,fcjde] 510 args = [[],[],arg] 520 Norm = norm.pdf(x,loc=pos,scale=widths[0]) 521 Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1]) 522 arg = [pos,shl/57.2958,dx,] 523 FCJ = fcjde.pdf(x,*arg,loc=pos) 524 if len(np.nonzero(FCJ)[0])>5: 525 z = np.column_stack([Norm,Cauchy,FCJ]).T 526 Z = fft(z) 527 Df = ifft(Z.prod(axis=0)).real 511 528 else: 512 fxns = [st.norm,st.cauchy]513 args = [[],[]]514 Df = calcPeakFFT(x,fxns,widths,pos,args)529 z = np.column_stack([Norm,Cauchy]).T 530 Z = fft(z) 531 Df = fftshift(ifft(Z.prod(axis=0))).real 515 532 Df /= np.sum(Df) 516 533 Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0) 517 534 return intens*Df(xdata)*DX/dx 518 535 519 536 def getPeakProfile(parmDict,xdata,varyList,bakType): 520 537 521 538 yb = getBackground('',parmDict,bakType,xdata) 522 539 yc = np.zeros_like(yb) 540 dx = xdata[1]xdata[0] 523 541 U = parmDict['U'] 524 542 V = parmDict['V'] … … 526 544 X = parmDict['X'] 527 545 Y = parmDict['Y'] 528 shl = max(parmDict['SH/L'],0.00 1)546 shl = max(parmDict['SH/L'],0.0005) 529 547 Ka2 = False 530 548 if 'Lam1' in parmDict.keys(): 531 549 Ka2 = True 532 lamRatio = parmDict['Lam2']/parmDict['Lam1']550 lamRatio = 360*(parmDict['Lam2']parmDict['Lam1'])/(np.pi*parmDict['Lam1']) 533 551 kRatio = parmDict['I(L2)/I(L1)'] 534 552 iPeak = 0 … … 548 566 else: 549 567 gam = X/cosd(pos/2.0)+Y*tand(pos/2.0) 550 gam = max(gam,0.0 1) #avoid neg gamma568 gam = max(gam,0.001) #avoid neg gamma 551 569 Wd,fmin,fmax = getWidths(pos,sig,gam,shl) 552 if pos > 90:553 fmin,fmax = [fmax,fmin]554 570 iBeg = np.searchsorted(xdata,posfmin) 555 iFin = np.searchsorted(xdata,pos+fmax) 571 lenX = len(xdata) 572 if not iBeg: 573 iFin = np.searchsorted(xdata,pos+fmin) 574 elif iBeg == lenX: 575 iFin = iBeg 576 else: 577 iFin = min(lenX,iBeg+int((fmin+fmax)/dx)) 556 578 if not iBeg+iFin: #peak below low limit 557 579 iPeak += 1 … … 561 583 yc[iBeg:iFin] += getFCJVoigt(pos,intens,sig,gam,shl,xdata[iBeg:iFin]) 562 584 if Ka2: 563 pos2 = 2.0*asind(lamRatio*sind(pos/2.0)) 564 Wd,fmin,fmax = getWidths(pos2,sig,gam,shl) 565 if pos > 90: 566 fmin,fmax = [fmax,fmin] 567 iBeg = np.searchsorted(xdata,pos2fmin) 568 iFin = np.searchsorted(xdata,pos2+fmax) 569 yc[iBeg:iFin] += getFCJVoigt(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin]) 585 pos2 = pos+lamRatio*tand(pos/2.0) # + 360/pi * Dlam/lam * tan(th) 586 kdelt = int((pos2pos)/dx) 587 iBeg = min(lenX,iBeg+kdelt) 588 iFin = min(lenX,iFin+kdelt) 589 if iBegiFin: 590 yc[iBeg:iFin] += getFCJVoigt(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin]) 570 591 iPeak += 1 571 592 except KeyError: #no more peaks to process 572 593 return yb+yc 573 574 def getWidths(pos,sig,gam,shl): 575 if shl: 576 widths = [np.sqrt(sig)/100.,gam/200.,.1] 577 else: 578 widths = [np.sqrt(sig)/100.,gam/200.] 579 fwhm = 2.355*widths[0]+2.*widths[1] 580 fmin = 10.*(fwhm+shl*abs(npcosd(pos))) 581 fmax = 15.0*fwhm 582 return widths,fmin,fmax 583 594 584 595 def Dict2Values(parmdict, varylist): 585 596 '''Use before call to leastsq to setup list of values for the parameters … … 635 646 insVary.append(insNames[i]) 636 647 instDict = dict(zip(insNames,insVals)) 637 instDict['X'] = max(instDict['X'],0. 1)638 instDict['Y'] = max(instDict['Y'],0. 1)639 instDict['SH/L'] = max(instDict['SH/L'],0.00 1)648 instDict['X'] = max(instDict['X'],0.01) 649 instDict['Y'] = max(instDict['Y'],0.01) 650 instDict['SH/L'] = max(instDict['SH/L'],0.0001) 640 651 return dataType,instDict,insVary 641 652 … … 735 746 return M 736 747 737 Ftol = 0.0 1748 Ftol = 0.0001 738 749 if oneCycle: 739 750 Ftol = 1.0 … … 759 770 dlg.SetPosition(wx.Point(screenSize[2]Size[0]305,screenSize[1]+5)) 760 771 try: 761 result = so.leastsq(errPeakProfile,values,full_output=True, ftol=Ftol,772 result = so.leastsq(errPeakProfile,values,full_output=True, #ftol=Ftol, 762 773 args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg)) 763 774 finally:
Note: See TracChangeset
for help on using the changeset viewer.