Changeset 939 for trunk/GSASIIpwd.py
 Timestamp:
 Jun 2, 2013 11:07:35 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIpwd.py
r808 r939 1 1 #/usr/bin/env python 2 2 # * coding: utf8 * 3 #GSASII powder calculation module 3 ''' 4 *GSASII powder calculation module* 5 ================================== 6 7 ''' 4 8 ########### SVN repository information ################### 5 9 # $Date$ … … 54 58 55 59 def Transmission(Geometry,Abs,Diam): 56 #Calculate sample transmission 57 # Geometry: one of 'Cylinder','BraggBrentano','Tilting flat plate in transmission','Fixed flat plate' 58 # Abs: absorption coeff in cm1 59 # Diam: sample thickness/diameter in mm 60 ''' 61 Calculate sample transmission 62 63 :param str Geometry: one of 'Cylinder','BraggBrentano','Tilting flat plate in transmission','Fixed flat plate' 64 :param float Abs: absorption coeff in cm1 65 :param float Diam: sample thickness/diameter in mm 66 ''' 60 67 if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2theta = 0; beam fully illuminates sample 61 68 MuR = Abs*Diam/20.0 … … 84 91 85 92 def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0): 86 #Calculate sample absorption 87 # Geometry: one of 'Cylinder','BraggBrentano','Tilting Flat Plate in transmission','Fixed flat plate' 88 # MuR: absorption coeff * sample thickness/2 or radius 89 # Tth: 2theta scattering angle  can be numpy array 90 # Phi: flat plate tilt angle  future 91 # Psi: flat plate tilt axis  future 93 '''Calculate sample absorption 94 :param str Geometry: one of 'Cylinder','BraggBrentano','Tilting Flat Plate in transmission','Fixed flat plate' 95 :param float MuR: absorption coeff * sample thickness/2 or radius 96 :param Tth: 2theta scattering angle  can be numpy array 97 :param float Phi: flat plate tilt angle  future 98 :param float Psi: flat plate tilt axis  future 99 ''' 92 100 Sth2 = npsind(Tth/2.0)**2 93 101 Cth2 = 1.Sth2 … … 127 135 128 136 def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0): 137 'needs a doc string' 129 138 dA = 0.001 130 139 AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi) … … 137 146 def Polarization(Pola,Tth,Azm=0.0): 138 147 """ Calculate angle dependent xray polarization correction (not scaled correctly!) 139 input: 140 Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized 141 Azm: azimuthal angle e.g. 0.0 in plane of polarization 142 Tth: 2theta scattering angle  can be numpy array 143 which (if either) of these is "right"? 144 return: 145 pola = ((1Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \ 146 (1Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2 147 dpdPola: derivative needed for least squares 148 149 :param Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized 150 :param Azm: azimuthal angle e.g. 0.0 in plane of polarization 151 :param Tth: 2theta scattering angle  can be numpy array 152 which (if either) of these is "right"? 153 :return: (pola, dpdPola) 154 * pola = ((1Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \ 155 (1Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2 156 * dpdPola: derivative needed for least squares 157 148 158 """ 149 159 pola = ((1.0Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \ … … 153 163 154 164 def Oblique(ObCoeff,Tth): 165 'needs a doc string' 155 166 if ObCoeff: 156 167 return (1.ObCoeff)/(1.0np.exp(np.log(ObCoeff)/npcosd(Tth))) … … 159 170 160 171 def Ruland(RulCoff,wave,Q,Compton): 172 'needs a doc string' 161 173 C = 2.9978e8 162 174 D = 1.5e3 … … 168 180 169 181 def LorchWeight(Q): 182 'needs a doc string' 170 183 return np.sin(np.pi*(Q[1]Q)/(2.0*Q[1])) 171 184 172 185 def GetAsfMean(ElList,Sthl2): 173 # Calculate various scattering factor terms for PDF calcs 174 # ElList: element dictionary contains scattering factor coefficients, etc. 175 # Sthl2: numpy array of sin theta/lambda squared values 176 # returns: mean(f^2), mean(f)^2, mean(compton) 186 '''Calculate various scattering factor terms for PDF calcs 187 188 :param dict ElList: element dictionary contains scattering factor coefficients, etc. 189 :param np.array Sthl2: numpy array of sin theta/lambda squared values 190 :returns: mean(f^2), mean(f)^2, mean(compton) 191 ''' 177 192 sumNoAtoms = 0.0 178 193 FF = np.zeros_like(Sthl2) … … 191 206 192 207 def GetNumDensity(ElList,Vol): 208 'needs a doc string' 193 209 sumNoAtoms = 0.0 194 210 for El in ElList: … … 197 213 198 214 def CalcPDF(data,inst,xydata): 215 'needs a doc string' 199 216 auxPlot = [] 200 217 import copy … … 303 320 def makeFFTsizeList(nmin=1,nmax=1023,thresh=15): 304 321 ''' Provide list of optimal data sizes for FFT calculations 305 Input: 306 nmin: minimum data size >= 1 307 nmax: maximum data size > nmin 308 thresh: maximum prime factor allowed 309 Returns: 310 list of data sizes where the maximum prime factor is < thresh 322 323 :param int nmin: minimum data size >= 1 324 :param int nmax: maximum data size > nmin 325 :param int thresh: maximum prime factor allowed 326 :Returns: list of data sizes where the maximum prime factor is < thresh 311 327 ''' 312 328 plist = [] … … 325 341 _norm_pdf_C = 1./math.sqrt(2*math.pi) 326 342 class norm_gen(st.rv_continuous): 327 343 'needs a doc string' 344 328 345 def pdf(self,x,*args,**kwds): 329 346 loc,scale=kwds['loc'],kwds['scale'] … … 346 363 347 364 class cauchy_gen(st.rv_continuous): 365 'needs a doc string' 348 366 349 367 def pdf(self,x,*args,**kwds): … … 369 387 FingerCoxJephcoat D(2phi,2th) function for S/L = H/L 370 388 Ref: J. Appl. Cryst. (1994) 27, 892900. 371 Parameters 372  373 x: array 1 to 1 374 t: 2theta position of peak 375 s: sum(S/L,H/L); S: sample height, H: detector opening, 376 L: sample to detector opening distance 377 dx: 2theta step size in deg 378 Result for fcj.pdf 379  380 T = x*dx+t 381 s = S/L+H/L 382 if x < 0: 383 fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}1)  1/s]/cos(T) 384 if x >= 0: 385 fcj.pdf = 0 389 390 :param x: array 1 to 1 391 :param t: 2theta position of peak 392 :param s: sum(S/L,H/L); S: sample height, H: detector opening, 393 L: sample to detector opening distance 394 :param dx: 2theta step size in deg 395 396 :returns: for fcj.pdf 397 398 * T = x*dx+t 399 * s = S/L+H/L 400 * if x < 0:: 401 402 fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}1)  1/s]/cos(T) 403 404 * if x >= 0: fcj.pdf = 0 386 405 """ 387 406 def _pdf(self,x,t,s,dx): … … 402 421 403 422 def getWidthsCW(pos,sig,gam,shl): 423 'needs a doc string' 404 424 widths = [np.sqrt(sig)/100.,gam/200.] 405 425 fwhm = 2.355*widths[0]+2.*widths[1] … … 411 431 412 432 def getWidthsTOF(pos,alp,bet,sig,gam): 433 'needs a doc string' 413 434 lnf = 3.3 # =log(0.001/2) 414 435 widths = [np.sqrt(sig),gam] … … 419 440 420 441 def getFWHM(TTh,Inst): 442 'needs a doc string' 421 443 sig = lambda Th,U,V,W: 1.17741*math.sqrt(max(0.001,U*tand(Th)**2+V*tand(Th)+W))*math.pi/180. 422 444 gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/180. … … 427 449 428 450 def getFCJVoigt(pos,intens,sig,gam,shl,xdata): 451 'needs a doc string' 429 452 DX = xdata[1]xdata[0] 430 453 widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl) … … 448 471 449 472 def getBackground(pfx,parmDict,bakType,xdata): 473 'needs a doc string' 450 474 yb = np.zeros_like(xdata) 451 475 nBak = 0 … … 529 553 530 554 def getBackgroundDerv(pfx,parmDict,bakType,xdata): 555 'needs a doc string' 531 556 nBak = 0 532 557 while True: … … 626 651 #use old fortran routine 627 652 def getFCJVoigt3(pos,sig,gam,shl,xdata): 653 'needs a doc string' 628 654 629 655 Df = pyd.pypsvfcj(len(xdata),xdatapos,pos,sig,gam,shl) … … 633 659 634 660 def getdFCJVoigt3(pos,sig,gam,shl,xdata): 661 'needs a doc string' 635 662 636 663 Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdatapos,pos,sig,gam,shl) … … 640 667 641 668 def getEpsVoigt(pos,alp,bet,sig,gam,xdata): 669 'needs a doc string' 642 670 Df = pyd.pyepsvoigt(len(xdata),xdatapos,alp,bet,sig,gam) 643 671 Df /= np.sum(Df) … … 645 673 646 674 def getdEpsVoigt(pos,alp,bet,sig,gam,xdata): 675 'needs a doc string' 647 676 Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdatapos,alp,bet,sig,gam) 648 677 sumDf = np.sum(Df) … … 650 679 651 680 def ellipseSize(H,Sij,GB): 681 'needs a doc string' 652 682 HX = np.inner(H.T,GB) 653 683 lenHX = np.sqrt(np.sum(HX**2)) … … 658 688 659 689 def ellipseSizeDerv(H,Sij,GB): 690 'needs a doc string' 660 691 lenR = ellipseSize(H,Sij,GB) 661 692 delt = 0.001 … … 668 699 669 700 def getHKLpeak(dmin,SGData,A): 701 'needs a doc string' 670 702 HKL = G2lat.GenHLaue(dmin,SGData,A) 671 703 HKLs = [] … … 677 709 678 710 def getPeakProfile(dataType,parmDict,xdata,varyList,bakType): 711 'needs a doc string' 679 712 680 713 yb = getBackground('',parmDict,bakType,xdata) … … 795 828 796 829 def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType): 830 'needs a doc string' 797 831 # needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order 798 832 dMdv = np.zeros(shape=(len(varyList),len(xdata))) … … 1006 1040 1007 1041 def SetBackgroundParms(Background): 1042 'needs a doc string' 1008 1043 if len(Background) == 1: # fix up old backgrounds 1009 1044 BackGround.append({'nDebye':0,'debyeTerms':[]}) … … 1047 1082 1048 1083 def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None): 1084 'needs a doc string' 1049 1085 1050 1086 … … 1347 1383 1348 1384 def calcIncident(Iparm,xdata): 1385 'needs a doc string' 1349 1386 1350 1387 def IfunAdv(Iparm,xdata): … … 1403 1440 NeedTestData = True 1404 1441 def TestData(): 1442 'needs a doc string' 1405 1443 # global NeedTestData 1406 1444 NeedTestData = False
Note: See TracChangeset
for help on using the changeset viewer.