Changeset 391
 Timestamp:
 Oct 13, 2011 10:48:51 AM (11 years ago)
 Location:
 trunk
 Files:

 5 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIIO.py
r384 r391 1269 1269 if 'aniso_label' not in labels: 1270 1270 while S: 1271 atom = ['','','',0,0,0, 0,'','','I',0.01,0,0,0,0,0,0]1271 atom = ['','','',0,0,0,1.0,'','','I',0.01,0,0,0,0,0,0] 1272 1272 S.strip() 1273 1273 if len(S.split()) != len(labels): 
trunk/GSASIIlattice.py
r384 r391 85 85 return G,g 86 86 87 def A2Gmat(A ):87 def A2Gmat(A,inverse=True): 88 88 """Fill real & reciprocal metric tensor (G) from A 89 89 90 :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23] 90 :param 91 A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23] 92 inverse: if True return bot G and g; else just G 91 93 :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays) 92 94 … … 97 99 [A[3]/2.,A[1], A[5]/2.], 98 100 [A[4]/2.,A[5]/2., A[2]]] 99 g = nl.inv(G) 100 return G,g 101 if inverse: 102 g = nl.inv(G) 103 return G,g 104 else: 105 return G 101 106 102 107 def Gmat2A(G): 
trunk/GSASIIplot.py
r390 r391 1372 1372 varyList = Data['varyList'] 1373 1373 Xmax = len(varyList) 1374 covArray = Data['covariance'] 1374 covMatrix = Data['covMatrix'] 1375 sig = np.sqrt(np.diag(covMatrix)) 1376 xvar = np.outer(sig,np.ones_like(sig)) 1377 covArray = np.divide(np.divide(covMatrix,xvar),xvar.T) 1375 1378 1376 1379 def OnPlotKeyPress(event): 
trunk/GSASIIpwdGUI.py
r390 r391 1292 1292 self.dataFrame.SelectPhase.Enable(True) 1293 1293 rowLabels = [] 1294 refList = [] 1295 for h,k,l,m,d,pos,sig,gam,fo,fc,phi,x,x,x in data[self.RefList]: 1296 refList.append([h,k,l,m,d,pos,sig,gam,fo,fc,phi]) 1294 refList = [refl[:11] for refl in data[self.RefList]] 1297 1295 for i in range(len(refList)): rowLabels.append(str(i)) 1298 1296 colLabels = ['H','K','L','mul','d','pos','sig','gam','Fosq','Fcsq','phase',] 
trunk/GSASIIstruct.py
r390 r391 473 473 return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables 474 474 475 def SetPhaseData(parmDict,sigDict,Phases): 475 def getVCov(varyNames,varyList,covMatrix): 476 vcov = np.zeros((len(varyNames),len(varyNames))) 477 for i1,name1 in enumerate(varyNames): 478 for i2,name2 in enumerate(varyNames): 479 try: 480 vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)] 481 except ValueError: 482 vcov[i1][i2] = 0.0 483 return vcov 484 485 def getCellEsd(pfx,SGData,A,covData): 486 dpr = 180./np.pi 487 rVsq = G2lat.calc_rVsq(A) 488 G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors 489 cell = np.array(G2lat.Gmat2cell(g)) #real cell 490 cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell 491 scos = cosd(cellst[3:6]) 492 ssin = sind(cellst[3:6]) 493 scot = scos/ssin 494 rcos = cosd(cell[3:6]) 495 rsin = sind(cell[3:6]) 496 rcot = rcos/rsin 497 RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5'] 498 varyList = covData['varyList'] 499 covMatrix = covData['covMatrix'] 500 vcov = getVCov(RMnames,varyList,covMatrix) 501 Ax = np.array(A) 502 Ax[3:] /= 2. 503 drVdA = np.array([Ax[1]*Ax[2]Ax[5]**2,Ax[0]*Ax[2]Ax[4]**2,Ax[0]*Ax[1]Ax[3]**2, 504 Ax[4]*Ax[5]Ax[2]*Ax[3],Ax[3]*Ax[5]Ax[1]*Ax[4],Ax[3]*Ax[4]Ax[0]*Ax[5]]) 505 srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T)) 506 Vol = 1/np.sqrt(rVsq) 507 sigVol = Vol**3*np.sqrt(srcvlsq)/2. 508 R123 = Ax[0]*Ax[1]*Ax[2] 509 dsasdg = np.zeros((3,6)) 510 dadg = np.zeros((6,6)) 511 for i0 in range(3): #0 1 2 512 i1 = (i0+1)%3 #1 2 0 513 i2 = (i1+1)%3 #2 0 1 514 i3 = 5i2 #3 5 4 515 i4 = 5i1 #4 3 5 516 i5 = 5i0 #5 4 3 517 dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1] 518 dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2] 519 dsasdg[i0][i5] = scot[i0]/np.sqrt(Ax[i1]*Ax[i2]) 520 denmsq = Ax[i0]*(R123Ax[i1]*Ax[i4]**2Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2) 521 denom = np.sqrt(denmsq) 522 dadg[i5][i0] = Ax[i5]/denomrcos[i0]/denmsq*(R1230.5*Ax[i1]*Ax[i4]**20.5*Ax[i2]*Ax[i3]**2) 523 dadg[i5][i1] = 0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]Ax[i0]*Ax[i4]**2) 524 dadg[i5][i2] = 0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]Ax[i0]*Ax[i3]**2) 525 dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]Ax[i3]*Ax[i4]**2) 526 dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]Ax[i3]**2*Ax[i4]) 527 dadg[i5][i5] = Ax[i0]/denom 528 for i0 in range(3): 529 i1 = (i0+1)%3 530 i2 = (i1+1)%3 531 i3 = 5i2 532 for ij in range(6): 533 dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]dsasdg[i1][ij]/ssin[i1]) 534 if ij == i0: 535 dadg[i0][ij] = dadg[i0][ij]0.5*cell[i0]/Ax[i0] 536 dadg[i3][ij] = dadg[i3][ij]*rsin[2i0]*dpr 537 sigMat = np.inner(dadg,np.inner(vcov,dadg.T)) 538 var = np.diag(sigMat) 539 CS = np.where(var>0.,np.sqrt(var),0.) 540 cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order 541 return cellSig 542 543 def SetPhaseData(parmDict,sigDict,Phases,covData): 476 544 477 545 def cellFill(pfx,SGData,parmDict,sigDict): … … 572 640 if cell[0]: 573 641 A,sigA = cellFill(pfx,SGData,parmDict,sigDict) 642 cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol 574 643 print ' Reciprocal metric tensor: ' 575 644 ptfmt = "%15.9f" … … 590 659 cell[1:7] = G2lat.A2cell(A) 591 660 cell[7] = G2lat.calc_V(A) 592 print ' New unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]), \ 593 ' c = %.5f'%(cell[3]),' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]), \ 594 ' gamma = %.3f'%(cell[6]),' volume = %.3f'%(cell[7]) 661 print ' New unit cell:' 662 ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"] 663 names = ['a','b','c','alpha','beta','gamma','Volume'] 664 namstr = ' names :' 665 ptstr = ' values:' 666 sigstr = ' esds :' 667 for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig): 668 namstr += '%12s'%(name) 669 ptstr += fmt%(a) 670 if siga: 671 sigstr += fmt%(siga) 672 else: 673 sigstr += 12*' ' 674 print namstr 675 print ptstr 676 print sigstr 677 595 678 if 'Pawley' in Phase['General']['Type']: 596 679 pawleyRef = Phase['Pawley ref'] … … 1762 1845 ymb = np.array(yyb) 1763 1846 ycmb = np.array(ycyb) 1764 ratio = np.where(ycmb >0.,ymb/ycmb,1.0)1847 ratio = np.where(ycmb!=0.,ymb/ycmb,0.0) 1765 1848 xB = np.searchsorted(x,Limits[0]) 1766 1849 xF = np.searchsorted(x,Limits[1]) … … 1789 1872 iFin2 = np.searchsorted(x,pos2+fmax) 1790 1873 yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here 1791 refl[8] = np.sum( yp[iBeg:iFin2]*ratio[iBeg:iFin2])/(refl[13]*(1.+kRatio))1874 refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0)) 1792 1875 elif 'T' in calcControls[hfx+'histType']: 1793 1876 print 'TOF Undefined at present' … … 2083 2166 covMatrix = result[1]*GOF 2084 2167 sig = np.sqrt(np.diag(covMatrix)) 2085 xvar = np.outer(sig,np.ones_like(sig))2086 cov = np.divide(np.divide(covMatrix,xvar),xvar.T)2087 2168 if np.any(np.isnan(sig)): 2088 2169 print '*** Least squares aborted  some invalid esds possible ***' … … 2106 2187 GetFobsSq(Histograms,Phases,parmDict,calcControls) 2107 2188 sigDict = dict(zip(varyList,sig)) 2108 SetPhaseData(parmDict,sigDict,Phases) 2189 covData = {'variables':result[0],'varyList':varyList,'covMatrix':covMatrix} 2190 SetPhaseData(parmDict,sigDict,Phases,covData) 2109 2191 SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms) 2110 2192 SetHistogramData(parmDict,sigDict,Histograms) 2111 covData = {'variables':result[0],'varyList':varyList,'covariance':cov}2112 2193 SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData) 2113 2194 #for testing purposes!!!
Note: See TracChangeset
for help on using the changeset viewer.