Changeset 357
 Timestamp:
 Aug 30, 2011 9:38:39 AM (11 years ago)
 Location:
 trunk
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIpwd.py
r355 r357 50 50 npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave) 51 51 52 53 def factorize(num): 54 ''' Provide prime number factors for integer num 55 Returns dictionary of prime factors (keys) & power for each (data) 56 ''' 57 factors = {} 58 orig = num 59 60 # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1 61 i, sqi = 2, 4 62 while sqi <= num: 63 while not num%i: 64 num /= i 65 factors[i] = factors.get(i, 0) + 1 66 67 sqi += 2*i + 1 68 i += 1 69 70 if num != 1 and num != orig: 71 factors[num] = factors.get(num, 0) + 1 72 73 if factors: 74 return factors 75 else: 76 return {num:1} #a prime number! 77 78 def makeFFTsizeList(nmin=1,nmax=1023,thresh=15): 79 ''' Provide list of optimal data sizes for FFT calculations 80 Input: 81 nmin: minimum data size >= 1 82 nmax: maximum data size > nmin 83 thresh: maximum prime factor allowed 84 Returns: 85 list of data sizes where the maximum prime factor is < thresh 86 ''' 87 plist = [] 88 nmin = max(1,nmin) 89 nmax = max(nmin+1,nmax) 90 for p in range(nmin,nmax): 91 if max(factorize(p).keys()) < thresh: 92 plist.append(p) 93 return plist 94 52 #GSASII pdf calculation routines 53 95 54 def Transmission(Geometry,Abs,Diam): 96 55 #Calculate sample transmission … … 176 135 # return (Pola*npcosd(Azm)**2+(1.Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \ 177 136 # Pola*npsind(Azm)**2+(1.Pola)*npcosd(Azm)**2 178 return((1.0Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \137 pola = ((1.0Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \ 179 138 (1.0Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2 139 dpdPola = npsind(Azm)**2*npsind(Tth)**2 140 return pola,dpdPola 180 141 181 142 def Oblique(ObCoeff,Tth): … … 407 368 Amu += el['FormulaNo']*el['mu'] 408 369 409 #GSASII peak fitting routine: Finger, Cox & Jephcoat model 370 def CalcPDF(data,inst,xydata): 371 auxPlot = [] 372 import copy 373 import scipy.fftpack as ft 374 #subtract backgrounds  if any 375 xydata['IofQ'] = copy.deepcopy(xydata['Sample']) 376 if data['Sample Bkg.']['Name']: 377 xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+ 378 data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult'] 379 if data['Container']['Name']: 380 xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult'] 381 if data['Container Bkg.']['Name']: 382 xycontainer += (xydata['Container Bkg.'][1][1]+ 383 data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult'] 384 xydata['IofQ'][1][1] += xycontainer 385 #get element data & absorption coeff. 386 ElList = data['ElList'] 387 Abs = G2lat.CellAbsorption(ElList,data['Form Vol']) 388 #Apply angle dependent corrections 389 Tth = xydata['Sample'][1][0] 390 dt = (Tth[1]Tth[0]) 391 xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth) 392 xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])[0] 393 if data['DetType'] == 'Image plate': 394 xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth) 395 XY = xydata['IofQ'][1] 396 #convert to Q 397 hc = 12.397639 398 if 'Lam' in inst: 399 wave = inst['Lam'] 400 else: 401 wave = inst['Lam1'] 402 keV = hc/wave 403 minQ = npT2q(Tth[0],wave) 404 maxQ = npT2q(Tth[1],wave) 405 Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True) 406 dq = Qpoints[1]Qpoints[0] 407 XY[0] = npT2q(XY[0],wave) 408 # Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9! 409 T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0) #OK for scipy 0.8 410 Qdata = T(Qpoints) 411 412 qLimits = data['QScaleLim'] 413 minQ = np.searchsorted(Qpoints,qLimits[0]) 414 maxQ = np.searchsorted(Qpoints,qLimits[1]) 415 newdata = [] 416 xydata['IofQ'][1][0] = Qpoints 417 xydata['IofQ'][1][1] = Qdata 418 for item in xydata['IofQ'][1]: 419 newdata.append(item[:maxQ]) 420 xydata['IofQ'][1] = newdata 421 422 423 xydata['SofQ'] = copy.deepcopy(xydata['IofQ']) 424 FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2) #these are <f^2>,<f>^2,Cf 425 Q = xydata['SofQ'][1][0] 426 ruland = Ruland(data['Ruland'],wave,Q,CF) 427 # auxPlot.append([Q,ruland,'Ruland']) 428 CF *= ruland 429 # auxPlot.append([Q,CF,'CFCorr']) 430 scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ]) 431 xydata['SofQ'][1][1] *= scale 432 xydata['SofQ'][1][1] = CF 433 xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF 434 scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ]) 435 xydata['SofQ'][1][1] *= scale 436 437 xydata['FofQ'] = copy.deepcopy(xydata['SofQ']) 438 xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]1.0) 439 if data['Lorch']: 440 xydata['FofQ'][1][1] *= LorchWeight(Q) 441 442 xydata['GofR'] = copy.deepcopy(xydata['FofQ']) 443 nR = len(xydata['GofR'][1][1]) 444 xydata['GofR'][1][1] = dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR]) 445 xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1] 446 447 448 return auxPlot 449 450 #GSASII peak fitting routines: Finger, Cox & Jephcoat model 451 452 def factorize(num): 453 ''' Provide prime number factors for integer num 454 Returns dictionary of prime factors (keys) & power for each (data) 455 ''' 456 factors = {} 457 orig = num 458 459 # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1 460 i, sqi = 2, 4 461 while sqi <= num: 462 while not num%i: 463 num /= i 464 factors[i] = factors.get(i, 0) + 1 465 466 sqi += 2*i + 1 467 i += 1 468 469 if num != 1 and num != orig: 470 factors[num] = factors.get(num, 0) + 1 471 472 if factors: 473 return factors 474 else: 475 return {num:1} #a prime number! 476 477 def makeFFTsizeList(nmin=1,nmax=1023,thresh=15): 478 ''' Provide list of optimal data sizes for FFT calculations 479 Input: 480 nmin: minimum data size >= 1 481 nmax: maximum data size > nmin 482 thresh: maximum prime factor allowed 483 Returns: 484 list of data sizes where the maximum prime factor is < thresh 485 ''' 486 plist = [] 487 nmin = max(1,nmin) 488 nmax = max(nmin+1,nmax) 489 for p in range(nmin,nmax): 490 if max(factorize(p).keys()) < thresh: 491 plist.append(p) 492 return plist 410 493 411 494 np.seterr(divide='ignore') … … 952 1035 PeaksPrint(parmDict,sigDict,varyList) 953 1036 954 def CalcPDF(data,inst,xydata):955 auxPlot = []956 import copy957 import scipy.fftpack as ft958 #subtract backgrounds  if any959 xydata['IofQ'] = copy.deepcopy(xydata['Sample'])960 if data['Sample Bkg.']['Name']:961 xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+962 data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']963 if data['Container']['Name']:964 xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']965 if data['Container Bkg.']['Name']:966 xycontainer += (xydata['Container Bkg.'][1][1]+967 data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']968 xydata['IofQ'][1][1] += xycontainer969 #get element data & absorption coeff.970 ElList = data['ElList']971 Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])972 #Apply angle dependent corrections973 Tth = xydata['Sample'][1][0]974 dt = (Tth[1]Tth[0])975 xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth)976 xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])977 if data['DetType'] == 'Image plate':978 xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)979 XY = xydata['IofQ'][1]980 #convert to Q981 hc = 12.397639982 if 'Lam' in inst:983 wave = inst['Lam']984 else:985 wave = inst['Lam1']986 keV = hc/wave987 minQ = npT2q(Tth[0],wave)988 maxQ = npT2q(Tth[1],wave)989 Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)990 dq = Qpoints[1]Qpoints[0]991 XY[0] = npT2q(XY[0],wave)992 # Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!993 T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0) #OK for scipy 0.8994 Qdata = T(Qpoints)995 996 qLimits = data['QScaleLim']997 minQ = np.searchsorted(Qpoints,qLimits[0])998 maxQ = np.searchsorted(Qpoints,qLimits[1])999 newdata = []1000 xydata['IofQ'][1][0] = Qpoints1001 xydata['IofQ'][1][1] = Qdata1002 for item in xydata['IofQ'][1]:1003 newdata.append(item[:maxQ])1004 xydata['IofQ'][1] = newdata1005 1006 1007 xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])1008 FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2) #these are <f^2>,<f>^2,Cf1009 Q = xydata['SofQ'][1][0]1010 ruland = Ruland(data['Ruland'],wave,Q,CF)1011 # auxPlot.append([Q,ruland,'Ruland'])1012 CF *= ruland1013 # auxPlot.append([Q,CF,'CFCorr'])1014 scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])1015 xydata['SofQ'][1][1] *= scale1016 xydata['SofQ'][1][1] = CF1017 xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF1018 scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])1019 xydata['SofQ'][1][1] *= scale1020 1021 xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])1022 xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]1.0)1023 if data['Lorch']:1024 xydata['FofQ'][1][1] *= LorchWeight(Q)1025 1026 xydata['GofR'] = copy.deepcopy(xydata['FofQ'])1027 nR = len(xydata['GofR'][1][1])1028 xydata['GofR'][1][1] = dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])1029 xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]1030 1031 1032 return auxPlot1033 1034 1037 #testing data 1035 1038 NeedTestData = True 
trunk/GSASIIstruct.py
r350 r357 57 57 58 58 def ShowControls(Controls): 59 print ' Controls:' 60 if Controls['bandWidth']: 61 print ' Band approximated least squares matrix refinement, width: ',Controls['bandWidth'] 62 else: 63 print ' Full matrix least squares refinement' 64 print ' Maximum number of refinement cycles: ',Controls['Ncycles'] 65 print ' Minimum sum(shift/esd)^2 for convergence: ','%.2f'%(Controls['minSumShftESD']) 66 print ' The Marquardt damping factor: ','%.2f'%(Controls['Marquardt']) 67 print ' Maximum allowed atom shift: ','%.2f'%(Controls['maxShift']),'A' 68 if Controls['restraintWeight']: 69 print ' The weights of the restraint histograms will be modified to normalize their contribution','\n' 70 else: 71 print ' The restraint histogram weights are fixed','\n' 59 print ' Least squares controls:' 60 print ' Derivative type: ',Controls['deriv type'] 61 print ' Minimum deltaM/M for convergence: ','%.2g'%(Controls['min dM/M']) 72 62 73 63 def GetPhaseNames(GPXfile): … … 328 318 if isum == jsum: 329 319 G2mv.StoreEquivalence(varyJ,(varyI,)) 320 def cellVary(pfx,SGData): 321 if SGData['SGLaue'] in ['1',]: 322 return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5'] 323 elif SGData['SGLaue'] in ['2/m',]: 324 if SGData['SGUniq'] == 'a': 325 return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3'] 326 elif SGData['SGUniq'] == 'b': 327 return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4'] 328 else: 329 return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5'] 330 elif SGData['SGLaue'] in ['mmm',]: 331 return [pfx+'A0',pfx+'A1',pfx+'A2'] 332 elif SGData['SGLaue'] in ['4/m','4/mmm']: 333 return [pfx+'A0',pfx+'A2'] 334 elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: 335 return [pfx+'A0',pfx+'A2'] 336 elif SGData['SGLaue'] in ['3R', '3mR']: 337 return [pfx+'A0',pfx+'A3'] 338 elif SGData['SGLaue'] in ['m3m','m3']: 339 return [pfx+'A0'] 340 330 341 331 342 def GetPhaseData(PhaseData): 332 343 333 def cellVary(pfx,SGData):334 if SGData['SGLaue'] in ['1',]:335 return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']336 elif SGData['SGLaue'] in ['2/m',]:337 if SGData['SGUniq'] == 'a':338 return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']339 elif SGData['SGUniq'] == 'b':340 return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']341 else:342 return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']343 elif SGData['SGLaue'] in ['mmm',]:344 return [pfx+'A0',pfx+'A1',pfx+'A2']345 elif SGData['SGLaue'] in ['4/m','4/mmm']:346 return [pfx+'A0',pfx+'A2']347 elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:348 return [pfx+'A0',pfx+'A2']349 elif SGData['SGLaue'] in ['3R', '3mR']:350 return [pfx+'A0',pfx+'A3']351 elif SGData['SGLaue'] in ['m3m','m3']:352 return [pfx+'A0']353 354 344 355 345 print ' Phases:' … … 498 488 for i,refl in enumerate(pawleyRef): 499 489 key = pfx+'PWLref:'+str(i) 500 refl[6] = parmDict[key]490 refl[6] = abs(parmDict[key]) #suppress negative Fsq 501 491 if key in sigDict: 502 492 refl[7] = sigDict[key] … … 692 682 ptlbls += '%12s' % (name) 693 683 ptstr += '%12.6f' % (hapData[4][i]) 694 sigstr += '%12.6f' % (mustrainSig[1][i]) 684 if mustrainSig[1][i]: 685 sigstr += '%12.6f' % (mustrainSig[1][i]) 686 else: 687 sigstr += 12*' ' 695 688 print ptlbls 696 689 print ptstr … … 779 772 def GetSampleParms(hId,Sample): 780 773 sampVary = [] 781 pfx = ':'+str(hId)+':'782 sampDict = { pfx+'Gonio. radius':Sample['Gonio. radius']}774 hfx = ':'+str(hId)+':' 775 sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius']} 783 776 Type = Sample['Type'] 784 777 if 'Bragg' in Type: #BraggBrentano 785 778 for item in ['Scale','Shift','Transparency']: #surface roughness?, diffuse scattering? 786 sampDict[ pfx+item] = Sample[item][0]779 sampDict[hfx+item] = Sample[item][0] 787 780 if Sample[item][1]: 788 sampVary.append( pfx+item)781 sampVary.append(hfx+item) 789 782 elif 'Debye' in Type: #DebyeScherrer 790 783 for item in ['Scale','Absorption','DisplaceX','DisplaceY']: 791 sampDict[ pfx+item] = Sample[item][0]784 sampDict[hfx+item] = Sample[item][0] 792 785 if Sample[item][1]: 793 sampVary.append( pfx+item)786 sampVary.append(hfx+item) 794 787 return Type,sampDict,sampVary 795 788 … … 883 876 884 877 def SetBackgroundParms(pfx,Background,parmDict,sigDict): 885 backVals = Background[3:]886 backSig = [0 for i in range(len (Background[3:]))]887 for i in range(len (backVals)):888 backVals[i] = parmDict[pfx+'Back:'+str(i)]878 lenBack = len(Background[3:]) 879 backSig = [0 for i in range(lenBack)] 880 for i in range(lenBack): 881 Background[3+i] = parmDict[pfx+'Back:'+str(i)] 889 882 if Background[1]: 890 883 backSig[i] = sigDict[pfx+'Back:'+str(i)] … … 972 965 hId = Histogram['hId'] 973 966 pfx = ':'+str(hId)+':' 974 975 967 Background = Histogram['Background'][0] 976 968 backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict) … … 1079 1071 1080 1072 def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse): 1073 costh = cosd(refl[5]/2.) 1081 1074 if calcControls[phfx+'SizeType'] == 'isotropic': 1082 gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*cos d(refl[5]/2.))1075 gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*costh) 1083 1076 elif calcControls[phfx+'SizeType'] == 'uniaxial': 1084 1077 H = np.array(refl[:3]) 1085 1078 P = np.array(calcControls[phfx+'SizeAxis']) 1086 1079 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1087 gam = (1.8*wave/np.pi) *parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1']1088 gam /= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2)*cosd(refl[5]/2.)1080 gam = (1.8*wave/np.pi)/parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1'] 1081 gam *= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2)/costh 1089 1082 else: #ellipsoidal crystallites 1090 1083 H = np.array(refl[:3]) 1091 gam += 1.8*wave/(np.pi*cos d(refl[5])*np.inner(H,np.inner(sizeEllipse,H)))1084 gam += 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H))) 1092 1085 if calcControls[phfx+'MustrainType'] == 'isotropic': 1093 1086 gam += 0.018*parmDict[phfx+'Mustrain:0']*tand(refl[5]/2.)/np.pi … … 1109 1102 def GetIntensityCorr(refl,phfx,hfx,calcControls,parmDict): 1110 1103 Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3] #scale*multiplicity 1111 Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth']) 1104 Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0] 1112 1105 1113 1106 return Icorr … … 1188 1181 elif not iBegiFin: #peak above high limit  done 1189 1182 return yc,yb 1190 yc[iBeg:iFin] += G2pwd.getFCJVoigt3(refl[5],Icorr*refl[8],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here1183 yc[iBeg:iFin] += Icorr*refl[8]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here 1191 1184 if Ka2: 1192 1185 pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) … … 1194 1187 iBeg = np.searchsorted(x,pos2fmin) 1195 1188 iFin = np.searchsorted(x,pos2+fmax) 1196 yc[iBeg:iFin] += G2pwd.getFCJVoigt3(pos2,Icorr*refl[8]*kRatio,refl[6],refl[7],shl,x[iBeg:iFin]) #and here1189 yc[iBeg:iFin] += Icorr*refl[8]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin]) #and here 1197 1190 else: 1198 1191 raise ValueError 1199 1192 return yc,yb 1200 1193 1201 1202 1194 def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): 1195 1196 def GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse): 1197 gamDict = {} 1198 costh = cosd(refl[5]/2.) 1199 tanth = tand(refl[5]/2.) 1200 if calcControls[phfx+'SizeType'] == 'isotropic': 1201 gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*costh) 1202 gamDict[phfx+'Size:0'] = gam/parmDict[phfx+'Size:0'] 1203 elif calcControls[phfx+'SizeType'] == 'uniaxial': 1204 H = np.array(refl[:3]) 1205 P = np.array(calcControls[phfx+'SizeAxis']) 1206 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1207 Si = parmDict[phfx+'Size:0'] 1208 Sa = parmDict[phfx+'Size:1'] 1209 gami = (1.8*wave/np.pi)/(Si*Sa) 1210 sqtrm = np.sqrt((cosP*Sa)**2+(sinP*Si)**2) 1211 gam = gami*sqtrm/costh 1212 gamDict[phfx+'Size:0'] = gami*Si*sinP**2/(sqtrm*costh)gam/Si 1213 gamDict[phfx+'Size:1'] = gami*Sa*cosP**2/(sqtrm*costh)gam/Sa 1214 else: #ellipsoidal crystallites  do numerically? 1215 H = np.array(refl[:3]) 1216 gam = 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H))) 1217 1218 if calcControls[phfx+'MustrainType'] == 'isotropic': 1219 gamDict[phfx+'Mustrain:0'] = 0.018*tanth/np.pi 1220 elif calcControls[phfx+'MustrainType'] == 'uniaxial': 1221 H = np.array(refl[:3]) 1222 P = np.array(calcControls[phfx+'MustrainAxis']) 1223 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1224 Si = parmDict[phfx+'Mustrain:0'] 1225 Sa = parmDict[phfx+'Mustrain:1'] 1226 gami = 0.018*Si*Sa*tanth/np.pi 1227 sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2) 1228 gam = gami/sqtrm 1229 gamDict[phfx+'Mustrain:0'] = gam/Sigami*Si*cosP**2/sqtrm**3 1230 gamDict[phfx+'Mustrain:1'] = gam/Sagami*Sa*sinP**2/sqtrm**3 1231 else: #generalized  P.W. Stephens model 1232 pwrs = calcControls[phfx+'MuPwrs'] 1233 const = 0.018*refl[4]**2*tanth 1234 for i,pwr in enumerate(pwrs): 1235 gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] 1236 return gamDict 1237 1238 def GetIntensityDerv(refl,phfx,hfx,calcControls,parmDict): 1239 Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3] #scale*multiplicity 1240 pola,dpdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth']) 1241 dIdpola = Icorr*dpdPola 1242 Icorr *= pola 1243 dIdsh = Icorr/parmDict[hfx+'Scale'] 1244 dIdsp = Icorr/parmDict[phfx+'Scale'] 1245 1246 return Icorr,dIdsh,dIdsp,dIdpola 1247 1248 def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict): 1249 dpr = 180./np.pi 1250 h,k,l = refl[:3] 1251 dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A) 1252 dst = np.sqrt(dstsq) 1253 pos = refl[5] 1254 const = dpr/np.sqrt(1.0wave*dst/4.0) 1255 dpdw = const*dst 1256 dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l]) 1257 dpdA *= const*wave/(2.0*dst) 1258 dpdZ = 1.0 1259 const = 9.e2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns 1260 if 'Bragg' in calcControls[hfx+'instType']: 1261 dpdSh = 4.*const*cosd(pos/2.0) 1262 dpdTr = 1.e7*const*sind(pos) 1263 return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0. 1264 else: #DebyeScherrer  simple but maybe not right 1265 dpdXd = const*cosd(pos) 1266 dpdYd = const*sind(pos) 1267 return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd 1268 1269 def cellVaryDerv(pfx,SGData,dpdA): 1270 if SGData['SGLaue'] in ['1',]: 1271 return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]], 1272 [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]] 1273 elif SGData['SGLaue'] in ['2/m',]: 1274 if SGData['SGUniq'] == 'a': 1275 return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]] 1276 elif SGData['SGUniq'] == 'b': 1277 return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]] 1278 else: 1279 return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]] 1280 elif SGData['SGLaue'] in ['mmm',]: 1281 return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]] 1282 elif SGData['SGLaue'] in ['4/m','4/mmm']: 1283 return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]] 1284 elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: 1285 return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]] 1286 elif SGData['SGLaue'] in ['3R', '3mR']: 1287 return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]] 1288 elif SGData['SGLaue'] in ['m3m','m3']: 1289 return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]] 1290 1291 lenX = len(x) 1292 hId = Histogram['hId'] 1293 hfx = ':%d:'%(hId) 1294 bakType = calcControls[hfx+'bakType'] 1295 dMdv = np.zeros(shape=(len(varylist),len(x))) 1296 if hfx+'Back:0' in varylist: 1297 dMdb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x) 1298 bBpos =varylist.index(hfx+'Back:0') 1299 dMdv[bBpos:bBpos+len(dMdb)] = dMdb 1300 1301 if 'C' in calcControls[hfx+'histType']: 1302 dx = x[1]x[0] 1303 shl = max(parmDict[hfx+'SH/L'],0.002) 1304 Ka2 = False 1305 if hfx+'Lam1' in parmDict.keys(): 1306 wave = parmDict[hfx+'Lam1'] 1307 Ka2 = True 1308 lamRatio = 360*(parmDict[hfx+'Lam2']parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1']) 1309 kRatio = parmDict[hfx+'I(L2)/I(L1)'] 1310 else: 1311 wave = parmDict[hfx+'Lam'] 1312 else: 1313 print 'TOF Undefined at present' 1314 raise ValueError 1315 for phase in Histogram['Reflection Lists']: 1316 refList = Histogram['Reflection Lists'][phase] 1317 Phase = Phases[phase] 1318 SGData = Phase['General']['SGData'] 1319 pId = Phase['pId'] 1320 pfx = '%d::'%(pId) 1321 phfx = '%d:%d:'%(pId,hId) 1322 A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] 1323 G,g = G2lat.A2Gmat(A) #recip & real metric tensors 1324 sizeEllipse = [] 1325 if calcControls[phfx+'SizeType'] == 'ellipsoidal': 1326 sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)]) 1327 for iref,refl in enumerate(refList): 1328 if 'C' in calcControls[hfx+'histType']: 1329 Icorr,dIdsh,dIdsp,dIdpola = GetIntensityDerv(refl,phfx,hfx,calcControls,parmDict) 1330 hkl = refl[:3] 1331 pos = refl[5] 1332 tanth = tand(pos/2.0) 1333 costh = cosd(pos/2.0) 1334 dsdU = tanth**2 1335 dsdV = tanth 1336 dsdW = 1.0 1337 dgdX = 1.0/costh 1338 dgdY = tanth 1339 Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) 1340 iBeg = np.searchsorted(x,refl[5]fmin) 1341 iFin = np.searchsorted(x,refl[5]+fmax) 1342 dMdpk = np.zeros(shape=(6,len(x))) 1343 dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) 1344 for i in range(1,5): 1345 dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[i] 1346 dMdpk[0][iBeg:iFin] += 100.*dx*Icorr*dMdipk[0] 1347 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]} 1348 if Ka2: 1349 pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) 1350 kdelt = int((pos2refl[5])/dx) 1351 iBeg = min(lenX,iBeg+kdelt) 1352 iFin = min(lenX,iFin+kdelt) 1353 if iBegiFin: 1354 dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,xdata[iBeg:iFin]) 1355 for i in range(1,5): 1356 dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[i] 1357 dMdpk[0][iBeg:iFin] += 100.*dx*kRatio*dMdipk2[0] 1358 dMdpk[5][iBeg:iFin] += 100.*dx*dMdipk2[0] 1359 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*Icorr*refl[8]} 1360 try: 1361 idx = varylist.index(pfx+'PWLref:'+str(iref)) 1362 dMdv[idx] = dervDict['int'] 1363 except ValueError: 1364 pass 1365 dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict) 1366 names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'], 1367 hfx+'U':[dsdU,'sig'],hfx+'V':[dsdV,'sig'],hfx+'W':[dsdW,'sig'], 1368 hfx+'X':[dgdX,'gam'],hfx+'Y':[dgdY,'gam'],hfx+'SH/L':[1.0,'shl'], 1369 hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'], 1370 hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'], 1371 hfx+'DisplaceY':[dpdY,'pos'],} 1372 for name in names: 1373 if name in varylist: 1374 item = names[name] 1375 dMdv[varylist.index(name)] += item[0]*dervDict[item[1]] 1376 cellDervNames = cellVaryDerv(pfx,SGData,dpdA) 1377 for name,dpdA in cellDervNames: 1378 if name in varylist: 1379 dMdv[varylist.index(name)] += dpdA*dervDict['pos'] 1380 gamDict = GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) 1381 for name in gamDict: 1382 if name in varylist: 1383 dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam'] 1384 else: 1385 raise ValueError 1386 1387 return dMdv 1388 1203 1389 def Refine(GPXfile): 1390 1391 def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): 1392 parmdict.update(zip(varylist,values)) 1393 G2mv.Dict2Map(parmDict) 1394 Histograms,Phases = HistoPhases 1395 dMdv = np.empty(0) 1396 for histogram in Histograms: 1397 if 'PWDR' in histogram[:4]: 1398 Histogram = Histograms[histogram] 1399 hId = Histogram['hId'] 1400 hfx = ':%d:'%(hId) 1401 Limits = calcControls[hfx+'Limits'] 1402 x,y,w,yc,yb,yd = Histogram['Data'] 1403 xB = np.searchsorted(x,Limits[0]) 1404 xF = np.searchsorted(x,Limits[1]) 1405 dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF], 1406 varylist,Histogram,Phases,calcControls,pawleyLookup) 1407 if len(dMdv): 1408 dMdv = np.concatenate((dMdv,dMdvh)) 1409 else: 1410 dMdv = dMdvh 1411 return dMdv 1204 1412 1205 1413 def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): … … 1229 1437 yc[xB:xF] *= parmDict[hfx+'Scale'] 1230 1438 yc[xB:xF] += yb[xB:xF] 1231 yd[xB:xF] = y [xB:xF]yc[xB:xF]1439 yd[xB:xF] = yc[xB:xF]y[xB:xF] #ycyo then all dydv have no '' needed 1232 1440 Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF])) 1233 1441 M = np.concatenate((M,np.sqrt(w[xB:xF])*(yd[xB:xF]))) … … 1245 1453 parmDict = {} 1246 1454 calcControls = {} 1247 #Controls = GetControls(GPXfile)1248 #ShowControls(Controls)1455 Controls = GetControls(GPXfile) 1456 ShowControls(Controls) 1249 1457 Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile) 1250 1458 if not Phases: … … 1278 1486 Size = dlg.GetSize() 1279 1487 dlg.SetPosition(wx.Point(screenSize[2]Size[0]305,screenSize[1]+5)) 1488 Ftol = Controls['min dM/M'] 1280 1489 try: 1281 result = so.leastsq(errRefine,values,full_output=True,ftol=0.0001,epsfcn=1.e8, 1282 args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) 1490 if Controls['deriv type'] == 'analytic': 1491 result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,ftol=Ftol,col_deriv=True, 1492 args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) 1493 ncyc = int(result[2]['nfev']/2) 1494 else: #'numeric' 1495 result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e8, 1496 args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) 1497 ncyc = int(result[2]['nfev']/len(varyList)) 1283 1498 finally: 1284 1499 dlg.Destroy() 1285 runtime = time.time()begin 1500 runtime = time.time()begin 1286 1501 chisq = np.sum(result[2]['fvec']**2) 1287 ncyc = int(result[2]['nfev']/len(varyList))1288 1502 Values2Dict(parmDict, varyList, result[0]) 1289 1503 G2mv.Dict2Map(parmDict) … … 1316 1530 SetHistogramData(parmDict,sigDict,Histograms) 1317 1531 SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases) 1532 #for testing purposes!!! 1533 # import cPickle 1534 # file = open('structTestdata.dat','wb') 1535 # cPickle.dump(parmDict,file,1) 1536 # cPickle.dump(varyList,file,1) 1537 # for histogram in Histograms: 1538 # if 'PWDR' in histogram[:4]: 1539 # Histogram = Histograms[histogram] 1540 # cPickle.dump(Histogram,file,1) 1541 # cPickle.dump(Phases,file,1) 1542 # cPickle.dump(calcControls,file,1) 1543 # cPickle.dump(pawleyLookup,file,1) 1544 # file.close() 1318 1545 1319 1546 def main():
Note: See TracChangeset
for help on using the changeset viewer.