Changeset 357 for trunk/GSASIIpwd.py
 Timestamp:
 Aug 30, 2011 9:38:39 AM (11 years ago)
 File:

 1 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
Note: See TracChangeset
for help on using the changeset viewer.