Changeset 1459
- Timestamp:
- Aug 8, 2014 2:59:09 PM (9 years ago)
- Location:
- trunk
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASII.py
r1453 r1459 979 979 return [G2IO.makeInstDict(names,data,codes),{}] 980 980 elif 'T' in DataType: 981 names = ['Type',' 2-theta','difC','difA','difB','Zero','alpha','beta-0','beta-1',982 'beta-q','sig-0','sig-1','sig-q','X', 'Y','Azimuth']983 codes = [0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0]981 names = ['Type','fltPath','2-theta','difC','difA', 'difB','Zero','alpha','beta-0','beta-1', 982 'beta-q','sig-0','sig-1','sig-q','X', 'Y','Azimuth',] 983 codes = [0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,] 984 984 azm = 0. 985 985 if 'INS 1DETAZM' in Iparm: 986 986 azm = float(Iparm['INS 1DETAZM']) 987 s = Iparm['INS FPATH1'].split() 988 fltPath0 = G2IO.sfloat(s[0]) 987 989 s = Iparm['INS 1BNKPAR'].split() 990 fltPath1 = G2IO.sfloat(s[0]) 991 data.extend([fltPath0+fltPath1,]) #Flight path source-sample-detector 988 992 data.extend([G2IO.sfloat(s[1]),]) #2-theta for bank 989 993 s = Iparm['INS 1 ICONS'].split() … … 1889 1893 self.Offset = [0.0,0.0] 1890 1894 self.delOffset = .02 1891 self.refOffset = -1 00.01895 self.refOffset = -1.0 1892 1896 self.refDelt = .01 1893 1897 self.Weight = False … … 2791 2795 if name2 == 'Peak List': 2792 2796 peaks = self.PatternTree.GetItemPyData(item2)['peaks'] 2793 file.write("%s \n" % (name+' Peak List')) 2797 file.write("%s \n" % (name+' Peak List')) 2798 if len(peaks[0]) == 8: 2799 file.write('%10s %12s %10s %10s\n'%('pos','int','sig','gam')) 2800 else: 2801 file.write('%10s %12s %10s %10s %10s %10s\n'%('pos','int','alp','bet','sig','gam')) 2794 2802 for peak in peaks: 2795 file.write("%10.5f %12.2f %10.3f %10.3f \n" % \ 2796 (peak[0],peak[2],peak[4],peak[6])) 2803 if len(peak) == 8: #CW 2804 file.write("%10.5f %12.2f %10.3f %10.3f \n" % \ 2805 (peak[0],peak[2],peak[4],peak[6])) 2806 else: #TOF - more cols 2807 file.write("%10.5f %12.2f %10.3f %10.3f %10.3f %10.3f\n" % \ 2808 (peak[0],peak[2],peak[4],peak[6],peak[8],peak[10])) 2797 2809 item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2) 2798 2810 item, cookie = self.PatternTree.GetNextChild(self.root, cookie) -
trunk/GSASIIIO.py
r1446 r1459 1694 1694 defaultIparm_lbl.append('10m TOF backscattering bank') 1695 1695 defaultIparms.append({ 1696 'INS FPATH1':' 9.00', 1696 1697 'INS HTYPE ':'PNT', 1697 1698 'INS 1 ICONS':' 5000.00 0.00 0.00', … … 1703 1704 defaultIparm_lbl.append('10m TOF 90deg bank') 1704 1705 defaultIparms.append({ 1706 'INS FPATH1':' 9.00', 1705 1707 'INS HTYPE ':'PNT', 1706 1708 'INS 1 ICONS':' 3500.00 0.00 0.00', … … 1712 1714 defaultIparm_lbl.append('63m POWGEN 90deg bank') 1713 1715 defaultIparms.append({ 1716 'INS FPATH1':' 60.00', 1714 1717 'INS HTYPE ':'PNT', 1715 1718 'INS 1 ICONS':' 22585.80 0.00 0.00', 1716 'INS 1BNKPAR':' 1.000090.000',1719 'INS 1BNKPAR':' 3.169 90.000', 1717 1720 'INS 1PRCF1 ':' 1 8 0.01000', 1718 1721 'INS 1PRCF11':' 0.000000E+00 1.000000E+00 3.000000E-02 4.000000E-03', -
trunk/GSASIIplot.py
r1451 r1459 874 874 if G2frame.SqrtPlot: 875 875 G2frame.delOffset = .002 876 G2frame.refOffset = -1 0.0876 G2frame.refOffset = -1.0 877 877 G2frame.refDelt = .001 878 878 else: 879 879 G2frame.delOffset = .02 880 G2frame.refOffset = -1 00.0880 G2frame.refOffset = -1.0 881 881 G2frame.refDelt = .01 882 882 newPlot = True … … 1886 1886 else: 1887 1887 Plot.set_ylabel(r'$\mathsf{\Delta}T/T$',fontsize=14) 1888 colors=['b','g','r','c','m','k']1889 1888 for ixy,xy in enumerate(XY): 1890 1889 X,Y = xy … … 1899 1898 Plot.errorbar(X,Y,ecolor='k',yerr=E) 1900 1899 Plot.plot(X,Y,'kx',picker=3) 1900 Plot.axhline(0.,color='r',linestyle='--') 1901 1901 if not newPlot: 1902 1902 Page.toolbar.push_current() -
trunk/GSASIIpwd.py
r1443 r1459 561 561 bakInt = si.interp1d(bakPos,bakVals,'linear') 562 562 yb = bakInt(xdata) 563 if 'difC' in parmDict:563 if pfx+'difC' in parmDict: 564 564 ff = 1. 565 565 else: … … 602 602 return yb 603 603 604 def getBackgroundDerv( pfx,parmDict,bakType,xdata):604 def getBackgroundDerv(hfx,parmDict,bakType,xdata): 605 605 'needs a doc string' 606 606 nBak = 0 607 607 while True: 608 key = pfx+'Back:'+str(nBak)608 key = hfx+'Back:'+str(nBak) 609 609 if key in parmDict: 610 610 nBak += 1 … … 612 612 break 613 613 dydb = np.zeros(shape=(nBak,len(xdata))) 614 dyddb = np.zeros(shape=(3*parmDict[ pfx+'nDebye'],len(xdata)))615 dydpk = np.zeros(shape=(4*parmDict[ pfx+'nPeaks'],len(xdata)))614 dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata))) 615 dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata))) 616 616 cw = np.diff(xdata) 617 617 cw = np.append(cw,cw[-1]) … … 650 650 np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.), 651 651 np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.)) 652 if 'difC' in parmDict:652 if hfx+'difC' in parmDict: 653 653 ff = 1. 654 654 else: 655 655 try: 656 wave = parmDict[ pfx+'Lam']656 wave = parmDict[hfx+'Lam'] 657 657 except KeyError: 658 wave = parmDict[ pfx+'Lam1']658 wave = parmDict[hfx+'Lam1'] 659 659 q = 4.0*np.pi*npsind(xdata/2.0)/wave 660 660 SQ = (q/(4*np.pi))**2 … … 664 664 while True: 665 665 try: 666 dbA = parmDict[pfx+'DebyeA:'+str(iD)] 667 dbR = parmDict[pfx+'DebyeR:'+str(iD)] 668 dbU = parmDict[pfx+'DebyeU:'+str(iD)] 666 if hfx+'difC' in parmDict: 667 q = 2*np.pi*parmDict[hfx+'difC']/xdata 668 dbA = parmDict[hfx+'DebyeA:'+str(iD)] 669 dbR = parmDict[hfx+'DebyeR:'+str(iD)] 670 dbU = parmDict[hfx+'DebyeU:'+str(iD)] 669 671 sqr = np.sin(q*dbR)/(q*dbR) 670 672 cqr = np.cos(q*dbR) … … 679 681 while True: 680 682 try: 681 pkP = parmDict[ pfx+'BkPkpos;'+str(iD)]682 pkI = parmDict[ pfx+'BkPkint;'+str(iD)]683 pkS = parmDict[ pfx+'BkPksig;'+str(iD)]684 pkG = parmDict[ pfx+'BkPkgam;'+str(iD)]683 pkP = parmDict[hfx+'BkPkpos;'+str(iD)] 684 pkI = parmDict[hfx+'BkPkint;'+str(iD)] 685 pkS = parmDict[hfx+'BkPksig;'+str(iD)] 686 pkG = parmDict[hfx+'BkPkgam;'+str(iD)] 685 687 shl = 0.002 686 688 Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl) … … 1172 1174 peakDsp = [] 1173 1175 peakWt = [] 1174 for peak in IndexPeaks:1176 for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]): 1175 1177 if peak[2] and peak[3]: 1176 1178 peakPos.append(peak[0]) 1177 1179 peakDsp.append(peak[8]) 1178 peakWt.append(1/ peak[1])1180 peakWt.append(1/sig**2) 1179 1181 peakPos = np.array(peakPos) 1180 1182 peakDsp = np.array(peakDsp) -
trunk/GSASIIpwdGUI.py
r1453 r1459 87 87 'SlitLen':0.0, #Slit length - in Q(A-1) 88 88 } 89 def SetupSampleLabels(histName,dataType ):89 def SetupSampleLabels(histName,dataType,histType): 90 90 '''Setup a list of labels and number formatting for use in 91 91 labeling sample parameters. … … 95 95 parms = [] 96 96 parms.append(['Scale','Histogram scale factor: ',[10,4]]) 97 parms.append(['Gonio. radius','Goniometer radius (mm): ',[10,3]]) 97 if 'C' in histType: 98 parms.append(['Gonio. radius','Goniometer radius (mm): ',[10,3]]) 98 99 if 'PWDR' in histName: 99 100 if dataType == 'Debye-Scherrer': … … 101 102 ['DisplaceY',u'Sample Y displ. || to beam (\xb5m): ',[10,3]], 102 103 ['Absorption',u'Sample absorption (\xb5\xb7r): ',[10,4]],] 104 if 'T' in histType: 105 parms[-1] = ['Absorption',u'Sample absorption (\xb5\xb7r/l): ',[10,4]] 103 106 elif dataType == 'Bragg-Brentano': 104 107 parms += [['Shift',u'Sample displacement(\xb5m): ',[10,4]], … … 927 930 if key in ['Type','U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha', 928 931 'beta-0','beta-1','beta-q','sig-0','sig-1','sig-q','Polariz.', 929 'Lam','Azimuth','2-theta',' difC','difA','difB','Zero','Lam1','Lam2']:932 'Lam','Azimuth','2-theta','fltPath','difC','difA','difB','Zero','Lam1','Lam2']: 930 933 good.append(key) 931 934 return good … … 964 967 G2frame.ErrorDialog('Can not calibrate','Index Peak List not indexed') 965 968 return 966 G2pwd.DoCalibInst(IndexPeaks [0],data)969 G2pwd.DoCalibInst(IndexPeaks,data) 967 970 UpdateInstrumentGrid(G2frame,data) 968 971 XY = [] … … 1255 1258 else: #time of flight (neutrons) 1256 1259 subSizer = wx.BoxSizer(wx.HORIZONTAL) 1260 subSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,' Fligth path: '),0,WACV) 1261 txt = '%8.3f'%(insVal['fltPath']) 1262 subSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,txt.strip()),0,WACV) 1263 labelLst.append('flight path') 1264 elemKeysLst.append(['fltpath',1]) 1265 dspLst.append([10,2]) 1266 refFlgElem.append(None) 1257 1267 subSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,' 2-theta: '),0,WACV) 1258 1268 txt = '%7.2f'%(insVal['2-theta']) … … 1273 1283 subSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,' alpha, beta: fixed by table'),0,WACV) 1274 1284 else: 1275 Items = ['difC','difA','difB','Zero','alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-q','X','Y'] 1285 Items = ['difC','difA','difB','Zero','alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-q','','X','Y'] 1286 mainSizer.Add((5,5),0) 1276 1287 mainSizer.Add(subSizer) 1288 mainSizer.Add((5,5),0) 1277 1289 for item in Items: 1290 if item == '': 1291 instSizer.Add((5,5),0) 1292 instSizer.Add((5,5),0) 1293 instSizer.Add((5,5),0) 1294 continue 1278 1295 nDig = (10,3) 1279 1296 fmt = '%10.3f' … … 1292 1309 refFlgElem.append([item,2]) 1293 1310 instSizer.Add(RefineBox(item),0,WACV) 1294 # if not ifHisto and item in ['difC','difA','difB','Zero',]:1295 # refFlgElem.append(None)1296 # instSizer.Add((5,5),0)1297 # else:1298 # refFlgElem.append([item,2])1299 # instSizer.Add(RefineBox(item),0,WACV)1300 1311 elif 'S' in insVal['Type']: #single crystal data 1301 1312 if 'C' in insVal['Type']: #constant wavelength … … 1527 1538 # Assemble a list of item labels 1528 1539 TextTable = {key:label for key,label,dig in 1529 SetupSampleLabels(hst,data.get('Type') )1540 SetupSampleLabels(hst,data.get('Type'),Inst['Type'][0]) 1530 1541 } 1531 1542 # get flexible labels … … 1648 1659 #reload(G2gd) 1649 1660 ###################################################################### 1661 Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId( 1662 G2frame,G2frame.PatternId, 'Instrument Parameters'))[0] 1650 1663 histName = G2frame.PatternTree.GetItemText(G2frame.PatternId) 1651 1664 if G2frame.dataDisplay: … … 1698 1711 #patch end 1699 1712 labelLst,elemKeysLst,dspLst,refFlgElem = [],[],[],[] 1700 parms = SetupSampleLabels(histName,data.get('Type') )1713 parms = SetupSampleLabels(histName,data.get('Type'),Inst['Type'][0]) 1701 1714 mainSizer = wx.BoxSizer(wx.VERTICAL) 1702 1715 topSizer = wx.BoxSizer(wx.HORIZONTAL) … … 1724 1737 nameSizer.Add(wx.StaticText(G2frame.dataDisplay,wx.ID_ANY,' Diffractometer type: '), 1725 1738 0,WACV) 1726 choices = ['Debye-Scherrer','Bragg-Brentano',] 1739 if 'T' in Inst['Type'][0]: 1740 choices = ['Debye-Scherrer',] 1741 else: 1742 choices = ['Debye-Scherrer','Bragg-Brentano',] 1727 1743 histoType = wx.ComboBox(G2frame.dataDisplay,wx.ID_ANY,value=data['Type'],choices=choices, 1728 1744 style=wx.CB_READONLY|wx.CB_DROPDOWN) … … 2123 2139 ibrav = bravaisSymb.index(controls[5]) 2124 2140 SGData = G2spc.SpcGroup(controls[13])[1] 2125 dmin = G2indx.getDmin(peaks )-0.0052141 dmin = G2indx.getDmin(peaks[0])-0.005 2126 2142 G2frame.HKL = G2pwd.getHKLpeak(dmin,SGData,A) 2127 2143 G2indx.IndexPeaks(peaks[0],G2frame.HKL) … … 2475 2491 if G2frame.dataDisplay: 2476 2492 G2frame.dataFrame.Clear() 2493 Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))[0] 2477 2494 rowLabels = [] 2478 2495 if HKLF: … … 2494 2511 except TypeError: 2495 2512 refList = np.array([refl[:11] for refl in data[G2frame.RefList]]) 2496 I100 = refList.T[8]*np.array([refl[1 3] for refl in data[G2frame.RefList]])2513 I100 = refList.T[8]*np.array([refl[11] for refl in data[G2frame.RefList]]) 2497 2514 Imax = np.max(I100) 2498 2515 if Imax: 2499 2516 I100 *= 100.0/Imax 2500 refs = np.vstack((refList.T[:11],I100)).T 2517 if 'C' in Inst['Type'][0]: 2518 refs = np.vstack((refList.T[:11],I100)).T 2519 elif 'T' in Inst['Type'][0]: 2520 refs = np.vstack((refList.T[:15],I100)).T 2521 2501 2522 for i in range(len(refs)): rowLabels.append(str(i)) 2523 Types = 4*[wg.GRID_VALUE_LONG,]+4*[wg.GRID_VALUE_FLOAT+':10,4',]+ \ 2524 2*[wg.GRID_VALUE_FLOAT+':10,2',]+[wg.GRID_VALUE_FLOAT+':10,3',]+ \ 2525 [wg.GRID_VALUE_FLOAT+':10,3',] 2502 2526 if HKLF: 2503 2527 colLabels = ['H','K','L','mul','d','Fosq','sig','Fcsq','FoTsq','FcTsq','phase','Ext',] 2504 2528 else: 2505 colLabels = ['H','K','L','mul','d','pos','sig','gam','Fosq','Fcsq','phase','I100',] 2506 Types = 4*[wg.GRID_VALUE_LONG,]+4*[wg.GRID_VALUE_FLOAT+':10,4',]+ \ 2507 2*[wg.GRID_VALUE_FLOAT+':10,2',]+[wg.GRID_VALUE_FLOAT+':10,3',]+ \ 2508 [wg.GRID_VALUE_FLOAT+':10,3',] 2529 if 'C' in Inst['Type'][0]: 2530 colLabels = ['H','K','L','mul','d','pos','sig','gam','Fosq','Fcsq','phase','Icorr','I100',] 2531 Types += [wg.GRID_VALUE_FLOAT+':10,3',] 2532 elif 'T' in Inst['Type'][0]: 2533 colLabels = ['H','K','L','mul','d','pos','sig','gam','Fosq','Fcsq','phase','Icorr','alp','bet','wave','I100',] 2534 Types += 4*[wg.GRID_VALUE_FLOAT+':10,3',] 2535 2509 2536 G2frame.PeakTable = G2gd.Table(refs,rowLabels=rowLabels,colLabels=colLabels,types=Types) 2510 2537 G2frame.dataFrame.SetLabel('Reflection List for '+phaseName) -
trunk/GSASIIstrIO.py
r1415 r1459 239 239 hId = Histogram['hId'] 240 240 hfx = ':%d:'%(hId) 241 keV = controlDict[hfx+'keV'] 242 for El in FFtables: 243 Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0]) 244 FP,FPP,Mu = G2el.FPcalc(Orbs, keV) 245 FFtables[El][hfx+'FP'] = FP 246 FFtables[El][hfx+'FPP'] = FPP 241 if 'X' in controlDict[hfx+'histType']: 242 keV = controlDict[hfx+'keV'] 243 for El in FFtables: 244 Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0]) 245 FP,FPP,Mu = G2el.FPcalc(Orbs, keV) 246 FFtables[El][hfx+'FP'] = FP 247 FFtables[El][hfx+'FPP'] = FPP 247 248 248 249 def GetPhaseNames(GPXfile): … … 1718 1719 limits = Histogram['Limits'][1] 1719 1720 inst = Histogram['Instrument Parameters'][0] 1720 Zero = inst['Zero'][ 1]1721 Zero = inst['Zero'][0] 1721 1722 if 'C' in inst['Type'][1]: 1722 1723 try: … … 1725 1726 wave = inst['Lam1'][1] 1726 1727 dmin = wave/(2.0*sind(limits[1]/2.0)) 1728 elif 'T' in inst['Type'][0]: 1729 dmin = limits[0]/inst['difC'][1] 1727 1730 pfx = str(pId)+':'+str(hId)+':' 1728 1731 for item in ['Scale','Extinction']: … … 1810 1813 pos = 2.0*asind(wave/(2.0*d))+Zero 1811 1814 if limits[0] < pos < limits[1]: 1812 refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0,0.0,0.0])1815 refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0]) 1813 1816 Uniq.append(uniq) 1814 1817 Phi.append(phi) 1815 else: 1816 raise ValueError 1818 elif 'T' in inst['Type'][0]: 1819 pos = inst['difC'][1]*d+inst['difA'][1]*d**2+inst['difB'][1]*d**3+Zero 1820 if limits[0] < pos < limits[1]: 1821 wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0]) 1822 refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,wave]) 1823 Uniq.append(uniq) 1824 Phi.append(phi) 1817 1825 Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{}} 1818 1826 elif 'HKLF' in histogram: … … 2218 2226 if len(Inst[item]) > 2 and Inst[item][2]: 2219 2227 insVary.append(insName) 2220 # instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001) 2221 # instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001) 2222 instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005) 2228 if 'C' in dataType: 2229 instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005) 2223 2230 return dataType,instDict,insVary 2224 2231 … … 2281 2288 def PrintInstParms(Inst): 2282 2289 print >>pFile,'\n Instrument Parameters:' 2283 ptlbls = ' name :'2284 ptstr = ' value :'2285 varstr = ' refine:'2286 2290 insKeys = Inst.keys() 2287 2291 insKeys.sort() 2288 for item in insKeys: 2289 if item not in ['Type','Source']: 2290 ptlbls += '%12s' % (item) 2291 ptstr += '%12.6f' % (Inst[item][1]) 2292 if item in ['Lam1','Lam2','Azimuth']: 2293 varstr += 12*' ' 2294 else: 2295 varstr += '%12s' % (str(bool(Inst[item][2]))) 2296 print >>pFile,ptlbls 2297 print >>pFile,ptstr 2298 print >>pFile,varstr 2292 iBeg = 0 2293 Ok = True 2294 while Ok: 2295 ptlbls = ' name :' 2296 ptstr = ' value :' 2297 varstr = ' refine:' 2298 iFin = min(iBeg+9,len(insKeys)) 2299 for item in insKeys[iBeg:iFin]: 2300 if item not in ['Type','Source']: 2301 ptlbls += '%12s' % (item) 2302 ptstr += '%12.6f' % (Inst[item][1]) 2303 if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]: 2304 varstr += 12*' ' 2305 else: 2306 varstr += '%12s' % (str(bool(Inst[item][2]))) 2307 print >>pFile,ptlbls 2308 print >>pFile,ptstr 2309 print >>pFile,varstr 2310 iBeg = iFin 2311 if iBeg == len(insKeys): 2312 Ok = False 2313 else: 2314 print >>pFile,'\n' 2299 2315 2300 2316 def PrintSampleParms(Sample): … … 2347 2363 Type,instDict,insVary = GetInstParms(hId,Inst) 2348 2364 controlDict[pfx+'histType'] = Type 2349 if pfx+'Lam1' in instDict: 2350 controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1'] 2351 else: 2352 controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam'] 2365 if 'XC' in Type: 2366 if pfx+'Lam1' in instDict: 2367 controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1'] 2368 else: 2369 controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam'] 2353 2370 histDict.update(instDict) 2354 2371 histVary += insVary … … 2383 2400 Inst = Histogram['Instrument Parameters'][0] 2384 2401 controlDict[pfx+'histType'] = Inst['Type'][0] 2385 histDict[pfx+'Lam'] = Inst['Lam'][1] 2386 controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam'] 2402 if 'X' in Inst['Type'][0]: 2403 histDict[pfx+'Lam'] = Inst['Lam'][1] 2404 controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam'] 2387 2405 return histVary,histDict,controlDict 2388 2406 … … 2498 2516 2499 2517 def PrintInstParmsSig(Inst,instSig): 2500 ptlbls = ' names :'2501 ptstr = ' value :'2502 sigstr = ' sig :'2503 2518 refine = False 2504 2519 insKeys = instSig.keys() 2505 2520 insKeys.sort() 2506 for name in insKeys: 2507 if name not in ['Type','Lam1','Lam2','Azimuth','Source']: 2508 ptlbls += '%12s' % (name) 2509 ptstr += '%12.6f' % (Inst[name][1]) 2510 if instSig[name]: 2511 refine = True 2512 sigstr += '%12.6f' % (instSig[name]) 2513 else: 2514 sigstr += 12*' ' 2515 if refine: 2516 print >>pFile,'\n Instrument Parameters:' 2517 print >>pFile,ptlbls 2518 print >>pFile,ptstr 2519 print >>pFile,sigstr 2521 iBeg = 0 2522 Ok = True 2523 while Ok: 2524 ptlbls = ' names :' 2525 ptstr = ' value :' 2526 sigstr = ' sig :' 2527 iFin = min(iBeg+9,len(insKeys)) 2528 for name in insKeys[iBeg:iFin]: 2529 if name not in ['Type','Lam1','Lam2','Azimuth','Source','fltPath']: 2530 ptlbls += '%12s' % (name) 2531 ptstr += '%12.6f' % (Inst[name][1]) 2532 if instSig[name]: 2533 refine = True 2534 sigstr += '%12.6f' % (instSig[name]) 2535 else: 2536 sigstr += 12*' ' 2537 if refine: 2538 print >>pFile,'\n Instrument Parameters:' 2539 print >>pFile,ptlbls 2540 print >>pFile,ptstr 2541 print >>pFile,sigstr 2542 iBeg = iFin 2543 if iBeg == len(insKeys): 2544 Ok = False 2520 2545 2521 2546 def PrintSampleParmsSig(Sample,sampleSig): -
trunk/GSASIIstrMath.py
r1456 r1459 573 573 for iref,refl in enumerate(refDict['RefList']): 574 574 if 'NT' in calcControls[hfx+'histType']: 575 FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[1 2])575 FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14]) 576 576 fbs = np.array([0,0]) 577 577 H = refl[:3] … … 665 665 SQfactor = 4.0*SQ*twopisq 666 666 if 'T' in calcControls[hfx+'histType']: 667 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12]) 667 if 'P' in calcControls[hfx+'histType']: 668 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14]) 669 else: 670 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12]) 668 671 FP = np.repeat(FP.T,len(SGT),axis=0) 669 672 FPP = np.repeat(FPP.T,len(SGT),axis=0) … … 710 713 mSize = len(Mdata) 711 714 FF = np.zeros(len(Tdata)) 712 if 'N ' in calcControls[hfx+'histType']:715 if 'NC' in calcControls[hfx+'histType']: 713 716 FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam']) 714 el se:717 elif 'X' in calcControls[hfx+'histType']: 715 718 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 716 719 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) … … 724 727 dFdbab = np.zeros((nRef,2)) 725 728 for iref,refl in enumerate(refDict['RefList']): 729 if 'T' in calcControls[hfx+'histType']: 730 FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12]) 726 731 H = np.array(refl[:3]) 727 732 SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2 … … 731 736 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 732 737 FF = refDict['FF']['FF'][iref].T[Tindx] 733 # FF = [refDict['FF'][iref][El] for El in Tdata]734 738 Uniq = np.inner(H,SGMT) 735 739 Phi = np.inner(H,SGT) … … 911 915 'Spherical harmonics texture' 912 916 IFCoup = 'Bragg' in calcControls[hfx+'instType'] 917 if 'T' in calcControls[hfx+'histType']: 918 tth = parmDict[hfx+'2-theta'] 919 else: 920 tth = refl[5] 913 921 odfCor = 1.0 914 922 H = refl[:3] … … 917 925 Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']] 918 926 phi,beta = G2lat.CrsAng(H,cell,SGData) 919 psi,gam,x,x = G2lat.SamAng( refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.927 psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs. 920 928 SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder']) 921 929 for item in SHnames: … … 929 937 def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict): 930 938 'Spherical harmonics texture derivatives' 939 if 'T' in calcControls[hfx+'histType']: 940 tth = parmDict[hfx+'2-theta'] 941 else: 942 tth = refl[5] 931 943 FORPI = 4.0*np.pi 932 944 IFCoup = 'Bragg' in calcControls[hfx+'instType'] … … 939 951 Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']] 940 952 phi,beta = G2lat.CrsAng(H,cell,SGData) 941 psi,gam,dPSdA,dGMdA = G2lat.SamAng( refl[5]/2.,Gangls,Sangls,IFCoup)953 psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) 942 954 SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder']) 943 955 for item in SHnames: … … 954 966 def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict): 955 967 'spherical harmonics preferred orientation (cylindrical symmetry only)' 968 if 'T' in calcControls[hfx+'histType']: 969 tth = parmDict[hfx+'2-theta'] 970 else: 971 tth = refl[5] 956 972 odfCor = 1.0 957 973 H = refl[:3] … … 965 981 IFCoup = False 966 982 phi,beta = G2lat.CrsAng(H,cell,SGData) 967 psi,gam,x,x = G2lat.SamAng( refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.983 psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs. 968 984 SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False) 969 985 for item in SHnames: … … 975 991 def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict): 976 992 'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)' 993 if 'T' in calcControls[hfx+'histType']: 994 tth = parmDict[hfx+'2-theta'] 995 else: 996 tth = refl[5] 977 997 FORPI = 12.5663706143592 978 998 odfCor = 1.0 … … 988 1008 IFCoup = False 989 1009 phi,beta = G2lat.CrsAng(H,cell,SGData) 990 psi,gam,x,x = G2lat.SamAng( refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.1010 psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs. 991 1011 SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False) 992 1012 for item in SHnames: … … 1039 1059 'Needs a doc string' 1040 1060 if 'Debye' in calcControls[hfx+'instType']: 1041 return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0) 1061 if 'T' in calcControls[hfx+'histType']: 1062 return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14],parmDict[hfx+'2-theta'],0,0) 1063 else: 1064 return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0) 1042 1065 else: 1043 1066 return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5]) … … 1046 1069 'Needs a doc string' 1047 1070 if 'Debye' in calcControls[hfx+'instType']: 1048 return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0) 1071 if 'T' in calcControls[hfx+'histType']: 1072 return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14],parmDict[hfx+'2-theta'],0,0) 1073 else: 1074 return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0) 1049 1075 else: 1050 1076 return G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5]) … … 1059 1085 Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict) 1060 1086 Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict) 1061 re fl[11] = Icorr1087 return Icorr 1062 1088 1063 1089 def GetIntensityDerv(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict): … … 1084 1110 return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb 1085 1111 1086 def GetSampleSigGam(refl,wave,G,GB, phfx,calcControls,parmDict):1112 def GetSampleSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict): 1087 1113 'Needs a doc string' 1088 costh = cosd(refl[5]/2.) 1089 #crystallite size 1090 if calcControls[phfx+'SizeType'] == 'isotropic': 1091 Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh) 1092 elif calcControls[phfx+'SizeType'] == 'uniaxial': 1093 H = np.array(refl[:3]) 1094 P = np.array(calcControls[phfx+'SizeAxis']) 1095 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1096 Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh) 1097 Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2) 1098 else: #ellipsoidal crystallites 1099 Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] 1100 H = np.array(refl[:3]) 1101 lenR = G2pwd.ellipseSize(H,Sij,GB) 1102 Sgam = 1.8*wave/(np.pi*costh*lenR) 1103 #microstrain 1104 if calcControls[phfx+'MustrainType'] == 'isotropic': 1105 Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi 1106 elif calcControls[phfx+'MustrainType'] == 'uniaxial': 1107 H = np.array(refl[:3]) 1108 P = np.array(calcControls[phfx+'MustrainAxis']) 1109 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1110 Si = parmDict[phfx+'Mustrain;i'] 1111 Sa = parmDict[phfx+'Mustrain;a'] 1112 Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2)) 1113 else: #generalized - P.W. Stephens model 1114 pwrs = calcControls[phfx+'MuPwrs'] 1115 sum = 0 1116 for i,pwr in enumerate(pwrs): 1117 sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] 1118 Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum 1114 if 'C' in calcControls[hfx+'histType']: 1115 costh = cosd(refl[5]/2.) 1116 #crystallite size 1117 if calcControls[phfx+'SizeType'] == 'isotropic': 1118 Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh) 1119 elif calcControls[phfx+'SizeType'] == 'uniaxial': 1120 H = np.array(refl[:3]) 1121 P = np.array(calcControls[phfx+'SizeAxis']) 1122 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1123 Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh) 1124 Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2) 1125 else: #ellipsoidal crystallites 1126 Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] 1127 H = np.array(refl[:3]) 1128 lenR = G2pwd.ellipseSize(H,Sij,GB) 1129 Sgam = 1.8*wave/(np.pi*costh*lenR) 1130 #microstrain 1131 if calcControls[phfx+'MustrainType'] == 'isotropic': 1132 Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi 1133 elif calcControls[phfx+'MustrainType'] == 'uniaxial': 1134 H = np.array(refl[:3]) 1135 P = np.array(calcControls[phfx+'MustrainAxis']) 1136 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1137 Si = parmDict[phfx+'Mustrain;i'] 1138 Sa = parmDict[phfx+'Mustrain;a'] 1139 Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2)) 1140 else: #generalized - P.W. Stephens model 1141 pwrs = calcControls[phfx+'MuPwrs'] 1142 sum = 0 1143 for i,pwr in enumerate(pwrs): 1144 sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] 1145 Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum 1146 elif 'T' in calcControls[hfx+'histType']: 1147 #crystallite size 1148 if calcControls[phfx+'SizeType'] == 'isotropic': 1149 Sgam = 1.e-4*parmDict[hfx+'difC']*parmDict[phfx+'Size;i'] 1150 elif calcControls[phfx+'SizeType'] == 'uniaxial': 1151 H = np.array(refl[:3]) 1152 P = np.array(calcControls[phfx+'SizeAxis']) 1153 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1154 Sgam = 1.e-4*parmDict[hfx+'difC']*(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']) 1155 Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2) 1156 else: #ellipsoidal crystallites 1157 Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] 1158 H = np.array(refl[:3]) 1159 lenR = G2pwd.ellipseSize(H,Sij,GB) 1160 Sgam = 1.e-4*parmDict[hfx+'difC']*(refl[4]**2*lenR) 1161 #microstrain 1162 if calcControls[phfx+'MustrainType'] == 'isotropic': 1163 Mgam = 1.e-6*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;i'] 1164 elif calcControls[phfx+'MustrainType'] == 'uniaxial': 1165 H = np.array(refl[:3]) 1166 P = np.array(calcControls[phfx+'MustrainAxis']) 1167 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1168 Si = parmDict[phfx+'Mustrain;i'] 1169 Sa = parmDict[phfx+'Mustrain;a'] 1170 Mgam = 1.e-6*parmDict[hfx+'difC']*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2) 1171 else: #generalized - P.W. Stephens model 1172 pwrs = calcControls[phfx+'MuPwrs'] 1173 sum = 0 1174 for i,pwr in enumerate(pwrs): 1175 sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] 1176 Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]**2*sum 1177 1119 1178 gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx'] 1120 1179 sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2 … … 1122 1181 return sig,gam 1123 1182 1124 def GetSampleSigGamDerv(refl,wave,G,GB, phfx,calcControls,parmDict):1183 def GetSampleSigGamDerv(refl,wave,G,GB,hfx,phfx,calcControls,parmDict): 1125 1184 'Needs a doc string' 1126 1185 gamDict = {} 1127 1186 sigDict = {} 1128 costh = cosd(refl[5]/2.) 1129 tanth = tand(refl[5]/2.) 1130 #crystallite size derivatives 1131 if calcControls[phfx+'SizeType'] == 'isotropic': 1132 Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh) 1133 gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh) 1134 sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2) 1135 elif calcControls[phfx+'SizeType'] == 'uniaxial': 1136 H = np.array(refl[:3]) 1137 P = np.array(calcControls[phfx+'SizeAxis']) 1138 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1139 Si = parmDict[phfx+'Size;i'] 1140 Sa = parmDict[phfx+'Size;a'] 1141 gami = (1.8*wave/np.pi)/(Si*Sa) 1142 sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2) 1143 Sgam = gami*sqtrm 1144 gam = Sgam/costh 1145 dsi = (gami*Si*cosP**2/(sqtrm*costh)-gam/Si) 1146 dsa = (gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa) 1147 gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx'] 1148 gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx'] 1149 sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2 1150 sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2 1151 else: #ellipsoidal crystallites 1152 const = 1.8*wave/(np.pi*costh) 1153 Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] 1154 H = np.array(refl[:3]) 1155 lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB) 1156 Sgam = 1.8*wave/(np.pi*costh*lenR) 1157 for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]): 1158 gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx'] 1159 sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2 1160 gamDict[phfx+'Size;mx'] = Sgam 1161 sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2 1162 1163 #microstrain derivatives 1164 if calcControls[phfx+'MustrainType'] == 'isotropic': 1165 Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi 1166 gamDict[phfx+'Mustrain;i'] = 0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi 1167 sigDict[phfx+'Mustrain;i'] = 0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2) 1168 elif calcControls[phfx+'MustrainType'] == 'uniaxial': 1169 H = np.array(refl[:3]) 1170 P = np.array(calcControls[phfx+'MustrainAxis']) 1171 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1172 Si = parmDict[phfx+'Mustrain;i'] 1173 Sa = parmDict[phfx+'Mustrain;a'] 1174 gami = 0.018*Si*Sa*tanth/np.pi 1175 sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2) 1176 Mgam = gami/sqtrm 1177 dsi = -gami*Si*cosP**2/sqtrm**3 1178 dsa = -gami*Sa*sinP**2/sqtrm**3 1179 gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx'] 1180 gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx'] 1181 sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2 1182 sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2 1183 else: #generalized - P.W. Stephens model 1184 pwrs = calcControls[phfx+'MuPwrs'] 1185 const = 0.018*refl[4]**2*tanth 1186 sum = 0 1187 for i,pwr in enumerate(pwrs): 1188 term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] 1189 sum += parmDict[phfx+'Mustrain:'+str(i)]*term 1190 gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain;mx'] 1191 sigDict[phfx+'Mustrain:'+str(i)] = \ 1192 2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2 1193 Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum 1194 for i in range(len(pwrs)): 1195 sigDict[phfx+'Mustrain:'+str(i)] *= Mgam 1196 gamDict[phfx+'Mustrain;mx'] = Mgam 1197 sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2 1187 if 'C' in calcControls[hfx+'histType']: 1188 costh = cosd(refl[5]/2.) 1189 tanth = tand(refl[5]/2.) 1190 #crystallite size derivatives 1191 if calcControls[phfx+'SizeType'] == 'isotropic': 1192 Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh) 1193 gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh) 1194 sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2) 1195 elif calcControls[phfx+'SizeType'] == 'uniaxial': 1196 H = np.array(refl[:3]) 1197 P = np.array(calcControls[phfx+'SizeAxis']) 1198 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1199 Si = parmDict[phfx+'Size;i'] 1200 Sa = parmDict[phfx+'Size;a'] 1201 gami = (1.8*wave/np.pi)/(Si*Sa) 1202 sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2) 1203 Sgam = gami*sqtrm 1204 gam = Sgam/costh 1205 dsi = (gami*Si*cosP**2/(sqtrm*costh)-gam/Si) 1206 dsa = (gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa) 1207 gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx'] 1208 gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx'] 1209 sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2 1210 sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2 1211 else: #ellipsoidal crystallites 1212 const = 1.8*wave/(np.pi*costh) 1213 Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] 1214 H = np.array(refl[:3]) 1215 lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB) 1216 Sgam = const/lenR 1217 for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]): 1218 gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx'] 1219 sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2 1220 gamDict[phfx+'Size;mx'] = Sgam 1221 sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2 1222 1223 #microstrain derivatives 1224 if calcControls[phfx+'MustrainType'] == 'isotropic': 1225 Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi 1226 gamDict[phfx+'Mustrain;i'] = 0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi 1227 sigDict[phfx+'Mustrain;i'] = 0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2) 1228 elif calcControls[phfx+'MustrainType'] == 'uniaxial': 1229 H = np.array(refl[:3]) 1230 P = np.array(calcControls[phfx+'MustrainAxis']) 1231 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1232 Si = parmDict[phfx+'Mustrain;i'] 1233 Sa = parmDict[phfx+'Mustrain;a'] 1234 gami = 0.018*Si*Sa*tanth/np.pi 1235 sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2) 1236 Mgam = gami/sqtrm 1237 dsi = -gami*Si*cosP**2/sqtrm**3 1238 dsa = -gami*Sa*sinP**2/sqtrm**3 1239 gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx'] 1240 gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx'] 1241 sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2 1242 sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2 1243 else: #generalized - P.W. Stephens model 1244 pwrs = calcControls[phfx+'MuPwrs'] 1245 const = 0.018*refl[4]**2*tanth 1246 sum = 0 1247 for i,pwr in enumerate(pwrs): 1248 term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] 1249 sum += parmDict[phfx+'Mustrain:'+str(i)]*term 1250 gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain;mx'] 1251 sigDict[phfx+'Mustrain:'+str(i)] = \ 1252 2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2 1253 Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum 1254 for i in range(len(pwrs)): 1255 sigDict[phfx+'Mustrain:'+str(i)] *= Mgam 1256 gamDict[phfx+'Mustrain;mx'] = Mgam 1257 sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2 1258 else: #'T'OF 1259 if calcControls[phfx+'SizeType'] == 'isotropic': 1260 Sgam = 1.e-4*parmDict[hfx+'difC']*parmDict[phfx+'Size;i'] 1261 gamDict[phfx+'Size;i'] = 1.e-4*parmDict[hfx+'difC']*parmDict[phfx+'Size;mx'] 1262 sigDict[phfx+'Size;i'] = 1.e-4*parmDict[hfx+'difC']*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2 1263 elif calcControls[phfx+'SizeType'] == 'uniaxial': 1264 const = 1.e-4*parmDict[hfx+'difC'] 1265 H = np.array(refl[:3]) 1266 P = np.array(calcControls[phfx+'SizeAxis']) 1267 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1268 Si = parmDict[phfx+'Size;i'] 1269 Sa = parmDict[phfx+'Size;a'] 1270 gami = const*(Si*Sa) 1271 sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2) 1272 Sgam = gami*sqtrm 1273 dsi = gami*Si*cosP**2/sqtrm-gam/Si 1274 dsa = gami*Sa*sinP**2/sqtrm-gam/Sa 1275 gamDict[phfx+'Size;i'] = const*parmDict[phfx+'Size;mx']*Sa 1276 gamDict[phfx+'Size;a'] = const*parmDict[phfx+'Size;mx']*Si 1277 sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2 1278 sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2 1279 else: #ellipsoidal crystallites 1280 const = 1.e-4*parmDict[hfx+'difC'] 1281 Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] 1282 H = np.array(refl[:3]) 1283 lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB) 1284 Sgam = const/lenR 1285 for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]): 1286 gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx'] 1287 sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2 1288 gamDict[phfx+'Size;mx'] = Sgam 1289 sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2 1290 1291 #microstrain derivatives 1292 if calcControls[phfx+'MustrainType'] == 'isotropic': 1293 Mgam = 1.e-6*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;i'] 1294 gamDict[phfx+'Mustrain;i'] = 1.e-6*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx'] 1295 sigDict[phfx+'Mustrain;i'] = 2.e-6*parmDict[hfx+'difC']*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2 1296 elif calcControls[phfx+'MustrainType'] == 'uniaxial': 1297 H = np.array(refl[:3]) 1298 P = np.array(calcControls[phfx+'MustrainAxis']) 1299 cosP,sinP = G2lat.CosSinAngle(H,P,G) 1300 Si = parmDict[phfx+'Mustrain;i'] 1301 Sa = parmDict[phfx+'Mustrain;a'] 1302 gami = 1.e-6*parmDict[hfx+'difC']*Si*Sa 1303 sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2) 1304 Mgam = gami/sqtrm 1305 dsi = -gami*Si*cosP**2/sqtrm**3 1306 dsa = -gami*Sa*sinP**2/sqtrm**3 1307 gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx'] 1308 gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx'] 1309 sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2 1310 sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2 1311 else: #generalized - P.W. Stephens model 1312 pwrs = calcControls[phfx+'MuPwrs'] 1313 const = 1.e-6*parmDict[hfx+'difC']*refl[4]**2 1314 sum = 0 1315 for i,pwr in enumerate(pwrs): 1316 term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] 1317 sum += parmDict[phfx+'Mustrain:'+str(i)]*term 1318 gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain;mx'] 1319 sigDict[phfx+'Mustrain:'+str(i)] = \ 1320 2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2 1321 Mgam = const*sum 1322 for i in range(len(pwrs)): 1323 sigDict[phfx+'Mustrain:'+str(i)] *= Mgam 1324 gamDict[phfx+'Mustrain;mx'] = Mgam 1325 sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2 1326 1198 1327 return sigDict,gamDict 1199 1328 … … 1205 1334 1206 1335 refl[4] = d 1207 pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero'] 1208 const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns 1209 if 'Bragg' in calcControls[hfx+'instType']: 1210 pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \ 1211 parmDict[hfx+'Transparency']*sind(pos)*100.0) #trans(=1/mueff) in cm 1212 else: #Debye-Scherrer - simple but maybe not right 1213 pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos)) 1336 if 'C' in calcControls[hfx+'histType']: 1337 pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero'] 1338 const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns 1339 if 'Bragg' in calcControls[hfx+'instType']: 1340 pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \ 1341 parmDict[hfx+'Transparency']*sind(pos)*100.0) #trans(=1/mueff) in cm 1342 else: #Debye-Scherrer - simple but maybe not right 1343 pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos)) 1344 elif 'T' in calcControls[hfx+'histType']: 1345 pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']*d**3+parmDict[hfx+'Zero'] 1346 #do I need sample position effects - maybe? 1214 1347 return pos 1215 1348 … … 1220 1353 dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A) 1221 1354 dst = np.sqrt(dstsq) 1222 pos = refl[5]-parmDict[hfx+'Zero'] 1223 const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0) 1224 dpdw = const*dst 1225 dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l]) 1226 dpdA *= const*wave/(2.0*dst) 1227 dpdZ = 1.0 1228 const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns 1229 if 'Bragg' in calcControls[hfx+'instType']: 1230 dpdSh = -4.*const*cosd(pos/2.0) 1231 dpdTr = -const*sind(pos)*100.0 1232 return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0. 1233 else: #Debye-Scherrer - simple but maybe not right 1234 dpdXd = -const*cosd(pos) 1235 dpdYd = -const*sind(pos) 1236 return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd 1355 dsp = 1./dst 1356 if 'C' in calcControls[hfx+'histType']: 1357 pos = refl[5]-parmDict[hfx+'Zero'] 1358 const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0) 1359 dpdw = const*dst 1360 dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l]) 1361 dpdA *= const*wave/(2.0*dst) 1362 dpdZ = 1.0 1363 const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns 1364 if 'Bragg' in calcControls[hfx+'instType']: 1365 dpdSh = -4.*const*cosd(pos/2.0) 1366 dpdTr = -const*sind(pos)*100.0 1367 return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0. 1368 else: #Debye-Scherrer - simple but maybe not right 1369 dpdXd = -const*cosd(pos) 1370 dpdYd = -const*sind(pos) 1371 return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd 1372 elif 'T' in calcControls[hfx+'histType']: 1373 dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l]) 1374 dpdZ = 1.0 1375 dpdDC = dsp 1376 dpdDA = dsp**2 1377 dpdDB = dsp**3 1378 return dpdA,dpdZ,dpdDC,dpdDA,dpdDB 1237 1379 1238 1380 def GetHStrainShift(refl,SGData,phfx,parmDict): … … 1307 1449 hfx = ':%d:'%(hId) 1308 1450 Limits = calcControls[hfx+'Limits'] 1309 shl = max(parmDict[hfx+'SH/L'],0.0005) 1310 Ka2 = False 1311 kRatio = 0.0 1312 if hfx+'Lam1' in parmDict.keys(): 1313 Ka2 = True 1314 lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1']) 1315 kRatio = parmDict[hfx+'I(L2)/I(L1)'] 1451 if 'C' in calcControls[hfx+'histType']: 1452 shl = max(parmDict[hfx+'SH/L'],0.0005) 1453 Ka2 = False 1454 kRatio = 0.0 1455 if hfx+'Lam1' in parmDict.keys(): 1456 Ka2 = True 1457 lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1']) 1458 kRatio = parmDict[hfx+'I(L2)/I(L1)'] 1316 1459 x,y,w,yc,yb,yd = Histogram['Data'] 1317 1460 xB = np.searchsorted(x,Limits[0]) … … 1338 1481 iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF)) 1339 1482 iFin2 = iFin 1340 yp[iBeg:iFin] = refl[11]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl, x[iBeg:iFin]) #>90% of time spent here1483 yp[iBeg:iFin] = refl[11]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #>90% of time spent here 1341 1484 if Ka2: 1342 1485 pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) … … 1348 1491 elif not iBeg2-iFin2: #peak above high limit - done 1349 1492 break 1350 yp[iBeg2:iFin2] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl, x[iBeg2:iFin2]) #and here1493 yp[iBeg2:iFin2] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2])) #and here 1351 1494 refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11]*(1.+kRatio)),0.0)) 1352 1495 elif 'T' in calcControls[hfx+'histType']: 1353 print 'TOF Undefined at present' 1354 raise Exception #no TOF yet 1496 yp = np.zeros_like(yb) 1497 Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7]) 1498 iBeg = max(xB,np.searchsorted(x,refl[5]-fmin)) 1499 iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF)) 1500 yp[iBeg:iFin] = refl[11]*refl[9]*G2pwd.getEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin])) #>90% of time spent here 1501 refl[8] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11],0.0)) 1355 1502 Fo = np.sqrt(np.abs(refl[8])) 1356 1503 Fc = np.sqrt(np.abs(refl[9])) … … 1370 1517 'Needs a doc string' 1371 1518 1372 def GetReflSigGam (refl,wave,G,GB,hfx,phfx,calcControls,parmDict):1519 def GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict): 1373 1520 U = parmDict[hfx+'U'] 1374 1521 V = parmDict[hfx+'V'] … … 1377 1524 Y = parmDict[hfx+'Y'] 1378 1525 tanPos = tand(refl[5]/2.0) 1379 Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB, phfx,calcControls,parmDict)1526 Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) 1380 1527 sig = U*tanPos**2+V*tanPos+W+Ssig #save peak sigma 1381 1528 sig = max(0.001,sig) … … 1383 1530 gam = max(0.001,gam) 1384 1531 return sig,gam 1532 1533 def GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict): 1534 sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4]**2+parmDict[hfx+'sig-q']*refl[4] 1535 gam = parmDict[hfx+'X']*refl[4]+parmDict[hfx+'Y']*refl[4]**2 1536 Ssig,Sgam = GetSampleSigGam(refl,0.0,G,GB,hfx,phfx,calcControls,parmDict) 1537 sig += Ssig #save peak sigma 1538 sig = max(0.001,sig) 1539 gam += Sgam #save peak gamma 1540 gam = max(0.001,gam) 1541 return sig,gam 1542 1543 def GetReflAlpBet(refl,hfx,parmDict): 1544 alp = parmDict[hfx+'alpha']/refl[4] 1545 bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4]**4+parmDict[hfx+'beta-q']/refl[4] 1546 return alp,bet 1385 1547 1386 1548 hId = Histogram['hId'] … … 1400 1562 else: 1401 1563 wave = parmDict[hfx+'Lam'] 1402 else:1403 print 'TOF Undefined at present - might be nothing need be done here'1404 1564 for phase in Histogram['Reflection Lists']: 1405 1565 refDict = Histogram['Reflection Lists'][phase] … … 1427 1587 Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.)) #Lorentz correction 1428 1588 refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift 1429 refl[6:8] = GetReflSigGam (refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam1430 GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[11]1431 refl[11] *= Vst*Lorenz1589 refl[6:8] = GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict) #peak sig & gam 1590 refl[11] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)*Vst*Lorenz 1591 1432 1592 if Phase['General'].get('doPawley'): 1433 1593 try: … … 1456 1616 yc[iBeg:iFin] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #and here 1457 1617 elif 'T' in calcControls[hfx+'histType']: 1458 print 'TOF Undefined at present' 1459 raise Exception #no TOF yet 1618 h,k,l = refl[:3] 1619 Uniq = np.inner(refl[:3],SGMT) 1620 refl[5] = GetReflPos(refl,0.0,G,hfx,calcControls,parmDict) #corrected reflection position 1621 Lorenz = sind(parmDict[hfx+'2-theta']/2)*refl[4]**4 #TOF Lorentz correction 1622 refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift 1623 refl[6:8] = GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict) #peak sig & gam 1624 refl[12:14] = GetReflAlpBet(refl,hfx,parmDict) 1625 refl[11] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)*Vst*Lorenz 1626 if Phase['General'].get('doPawley'): 1627 try: 1628 pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]) 1629 refl[9] = parmDict[pInd] 1630 except KeyError: 1631 # print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l) 1632 continue 1633 Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7]) 1634 iBeg = np.searchsorted(x,refl[5]-fmin) 1635 iFin = np.searchsorted(x,refl[5]+fmax) 1636 if not iBeg+iFin: #peak below low limit - skip peak 1637 continue 1638 elif not iBeg-iFin: #peak above high limit - done 1639 break 1640 yc[iBeg:iFin] += refl[11]*refl[9]*G2pwd.getEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin])) 1460 1641 # print 'profile calc time: %.3fs'%(time.time()-time0) 1461 1642 return yc,yb … … 1518 1699 cw = np.diff(x) 1519 1700 cw = np.append(cw,cw[-1]) 1701 Ka2 = False #also for TOF! 1520 1702 if 'C' in calcControls[hfx+'histType']: 1521 1703 shl = max(parmDict[hfx+'SH/L'],0.002) 1522 Ka2 = False1523 1704 if hfx+'Lam1' in parmDict.keys(): 1524 1705 wave = parmDict[hfx+'Lam1'] … … 1528 1709 else: 1529 1710 wave = parmDict[hfx+'Lam'] 1530 else:1531 print 'TOF Undefined at present'1532 raise ValueError1533 1711 for phase in Histogram['Reflection Lists']: 1534 1712 refDict = Histogram['Reflection Lists'][phase] … … 1549 1727 time0 = time.time() 1550 1728 for iref,refl in enumerate(refDict['RefList']): 1729 h,k,l = refl[:3] 1730 Uniq = np.inner(refl[:3],SGMT) 1731 dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) 1551 1732 if 'C' in calcControls[hfx+'histType']: #CW powder 1552 h,k,l = refl[:3]1553 Uniq = np.inner(refl[:3],SGMT)1554 dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)1555 1733 Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl) 1556 iBeg = np.searchsorted(x,refl[5]-fmin) 1557 iFin = np.searchsorted(x,refl[5]+fmax) 1558 if not iBeg+iFin: #peak below low limit - skip peak 1559 continue 1560 elif not iBeg-iFin: #peak above high limit - done 1561 break 1562 pos = refl[5] 1734 else: #'T'OF 1735 Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7]) 1736 iBeg = np.searchsorted(x,refl[5]-fmin) 1737 iFin = np.searchsorted(x,refl[5]+fmax) 1738 if not iBeg+iFin: #peak below low limit - skip peak 1739 continue 1740 elif not iBeg-iFin: #peak above high limit - done 1741 break 1742 pos = refl[5] 1743 if 'C' in calcControls[hfx+'histType']: 1563 1744 tanth = tand(pos/2.0) 1564 1745 costh = cosd(pos/2.0) … … 1581 1762 dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11]*dMdipk2[0] 1582 1763 dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]} 1583 if Phase['General'].get('doPawley'): 1584 dMdpw = np.zeros(len(x)) 1585 try: 1586 pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]) 1587 idx = varylist.index(pIdx) 1588 dMdpw[iBeg:iFin] = dervDict['int']/refl[9] 1589 if Ka2: 1590 dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9] 1591 dMdv[idx] = dMdpw 1592 except: # ValueError: 1593 pass 1764 else: #'T'OF 1765 lenBF = iFin-iBeg 1766 dMdpk = np.zeros(shape=(6,lenBF)) 1767 dMdipk = G2pwd.getdEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin])) 1768 for i in range(5): 1769 dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11]*refl[9]*dMdipk[i] 1770 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]} 1771 if Phase['General'].get('doPawley'): 1772 dMdpw = np.zeros(len(x)) 1773 try: 1774 pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]) 1775 idx = varylist.index(pIdx) 1776 dMdpw[iBeg:iFin] = dervDict['int']/refl[9] 1777 if Ka2: #not for TOF either 1778 dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9] 1779 dMdv[idx] = dMdpw 1780 except: # ValueError: 1781 pass 1782 if 'C' in calcControls[hfx+'histType']: 1594 1783 dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict) 1595 1784 names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'], … … 1604 1793 else: 1605 1794 names.update({hfx+'Absorption':[dFdAb,'int'],}) 1606 for name in names: 1607 item = names[name] 1608 if name in varylist: 1609 dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]] 1795 else: #'T'OF 1796 dpdA,dpdZ,dpdDC,dpdDA,dpdDB = GetReflPosDerv(refl,0.0,A,hfx,calcControls,parmDict) 1797 names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'], 1798 hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'], 1799 hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4],'gam'],hfx+'Y':[refl[4]**2,'gam'], 1800 hfx+'alpha':[1./refl[4],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4]**4,'bet'], 1801 hfx+'beta-q':[1./refl[4],'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4]**2,'sig'], 1802 hfx+'sig-q':[refl[4],'sig'],hfx+'Absorption':[dFdAb,'int'],} 1803 for name in names: 1804 item = names[name] 1805 if name in varylist: 1806 dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]] 1807 if Ka2: 1808 dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]] 1809 elif name in dependentVars: 1810 depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]] 1811 if Ka2: 1812 depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]] 1813 for iPO in dIdPO: 1814 if iPO in varylist: 1815 dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int'] 1816 if Ka2: 1817 dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int'] 1818 elif iPO in dependentVars: 1819 depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int'] 1820 if Ka2: 1821 depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int'] 1822 for i,name in enumerate(['omega','chi','phi']): 1823 aname = pfx+'SH '+name 1824 if aname in varylist: 1825 dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int'] 1826 if Ka2: 1827 dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int'] 1828 elif aname in dependentVars: 1829 depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int'] 1830 if Ka2: 1831 depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int'] 1832 for iSH in dFdODF: 1833 if iSH in varylist: 1834 dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int'] 1835 if Ka2: 1836 dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int'] 1837 elif iSH in dependentVars: 1838 depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int'] 1839 if Ka2: 1840 depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int'] 1841 cellDervNames = cellVaryDerv(pfx,SGData,dpdA) 1842 for name,dpdA in cellDervNames: 1843 if name in varylist: 1844 dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos'] 1845 if Ka2: 1846 dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos'] 1847 elif name in dependentVars: 1848 depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos'] 1849 if Ka2: 1850 depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos'] 1851 dDijDict = GetHStrainShiftDerv(refl,SGData,phfx) 1852 for name in dDijDict: 1853 if name in varylist: 1854 dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos'] 1855 if Ka2: 1856 dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos'] 1857 elif name in dependentVars: 1858 depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos'] 1859 if Ka2: 1860 depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos'] 1861 if 'C' in calcControls[hfx+'histType']: 1862 sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) 1863 else: #'T'OF 1864 sigDict,gamDict = GetSampleSigGamDerv(refl,0.0,G,GB,hfx,phfx,calcControls,parmDict) 1865 for name in gamDict: 1866 if name in varylist: 1867 dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam'] 1868 if Ka2: 1869 dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam'] 1870 elif name in dependentVars: 1871 depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam'] 1872 if Ka2: 1873 depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam'] 1874 for name in sigDict: 1875 if name in varylist: 1876 dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig'] 1877 if Ka2: 1878 dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig'] 1879 elif name in dependentVars: 1880 depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig'] 1881 if Ka2: 1882 depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig'] 1883 for name in ['BabA','BabU']: 1884 if refl[9]: 1885 if phfx+name in varylist: 1886 dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9] 1610 1887 if Ka2: 1611 dMdv[varylist.index( name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]1612 elif name in dependentVars:1613 depDerivDict[ name][iBeg:iFin] += item[0]*dervDict[item[1]]1888 dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9] 1889 elif phfx+name in dependentVars: 1890 depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9] 1614 1891 if Ka2: 1615 depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]] 1616 for iPO in dIdPO: 1617 if iPO in varylist: 1618 dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int'] 1619 if Ka2: 1620 dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int'] 1621 elif iPO in dependentVars: 1622 depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int'] 1623 if Ka2: 1624 depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int'] 1625 for i,name in enumerate(['omega','chi','phi']): 1626 aname = pfx+'SH '+name 1627 if aname in varylist: 1628 dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int'] 1629 if Ka2: 1630 dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int'] 1631 elif aname in dependentVars: 1632 depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int'] 1633 if Ka2: 1634 depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int'] 1635 for iSH in dFdODF: 1636 if iSH in varylist: 1637 dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int'] 1638 if Ka2: 1639 dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int'] 1640 elif iSH in dependentVars: 1641 depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int'] 1642 if Ka2: 1643 depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int'] 1644 cellDervNames = cellVaryDerv(pfx,SGData,dpdA) 1645 for name,dpdA in cellDervNames: 1646 if name in varylist: 1647 dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos'] 1648 if Ka2: 1649 dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos'] 1650 elif name in dependentVars: 1651 depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos'] 1652 if Ka2: 1653 depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos'] 1654 dDijDict = GetHStrainShiftDerv(refl,SGData,phfx) 1655 for name in dDijDict: 1656 if name in varylist: 1657 dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos'] 1658 if Ka2: 1659 dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos'] 1660 elif name in dependentVars: 1661 depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos'] 1662 if Ka2: 1663 depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos'] 1664 sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict) 1665 for name in gamDict: 1666 if name in varylist: 1667 dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam'] 1668 if Ka2: 1669 dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam'] 1670 elif name in dependentVars: 1671 depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam'] 1672 if Ka2: 1673 depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam'] 1674 for name in sigDict: 1675 if name in varylist: 1676 dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig'] 1677 if Ka2: 1678 dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig'] 1679 elif name in dependentVars: 1680 depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig'] 1681 if Ka2: 1682 depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig'] 1683 for name in ['BabA','BabU']: 1684 if refl[9]: 1685 if phfx+name in varylist: 1686 dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9] 1687 if Ka2: 1688 dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9] 1689 elif phfx+name in dependentVars: 1690 depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9] 1691 if Ka2: 1692 depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9] 1693 elif 'T' in calcControls[hfx+'histType']: 1694 print 'TOF Undefined at present' 1695 raise Exception #no TOF yet 1892 depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9] 1696 1893 if not Phase['General'].get('doPawley'): 1697 1894 #do atom derivatives - for RB,F,X & U so far … … 1716 1913 if Ka2: 1717 1914 depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2 1718 # print 'profile derv time: %.3fs'%(time.time()-time0)1915 # print 'profile derv time: %.3fs'%(time.time()-time0) 1719 1916 # now process derivatives in constraints 1720 1917 G2mv.Dict2Deriv(varylist,depDerivDict,dMdv) … … 1722 1919 1723 1920 def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict): 1724 '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting1921 '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting 1725 1922 model (M) with respect to all parameters. Independent and dependant dM/dp arrays 1726 1923 are returned to either dervRefine or HessRefine. … … 1754 1951 for j,var in enumerate(varylist): 1755 1952 if var in dFdvDict: 1756 dMdvh[j][iref] = w*dFdvDict[var][iref]* dervCor*parmDict[phfx+'Scale']1953 dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor*ref[11] 1757 1954 for var in dependentVars: 1758 1955 if var in dFdvDict: 1759 depDerivDict[var][iref] = w*dFdvDict[var][iref]* dervCor*parmDict[phfx+'Scale']1956 depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor*ref[11] 1760 1957 if phfx+'Scale' in varylist: 1761 1958 dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor … … 1783 1980 for j,var in enumerate(varylist): 1784 1981 if var in dFdvDict: 1785 dMdvh[j][iref] = w*dFdvDict[var][iref]* dervCor*parmDict[phfx+'Scale']1982 dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor*ref[11] 1786 1983 for var in dependentVars: 1787 1984 if var in dFdvDict: 1788 depDerivDict[var][iref] = w*dFdvDict[var][iref]* dervCor*parmDict[phfx+'Scale']1985 depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor*ref[11] 1789 1986 if phfx+'Scale' in varylist: 1790 1987 dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
Note: See TracChangeset
for help on using the changeset viewer.