Changeset 2546
- Timestamp:
- Nov 22, 2016 1:08:48 PM (7 years ago)
- Location:
- trunk
- Files:
-
- 26 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIElem.py
r2492 r2546 15 15 16 16 import math 17 import sys 17 18 import os.path 18 19 import GSASIIpath -
trunk/GSASIIElemGUI.py
r2486 r2546 59 59 def ElButton(self, name, pos, tip, color): 60 60 'Creates an element button widget' 61 Black = wx.Colour(0,0,0)62 61 if not self.ifNone and name[0] == 'None': 63 62 return -
trunk/GSASIIIO.py
r2544 r2546 602 602 fileSize = os.stat(filename).st_size 603 603 Npix = (fileSize-6000)/2 604 Head =File.read(6000)604 File.read(6000) 605 605 head = ['Rigaku R-Axis IV detector data',] 606 606 image = np.array(ar.array('H',File.read(fileSize-6000)),dtype=np.int32) … … 2498 2498 trans = False 2499 2499 Trans = [] 2500 instr = False2501 2500 for diff in lines: 2502 2501 diff = diff[:-1].lower() … … 2513 2512 continue 2514 2513 if diff.strip() == 'instrumental': 2515 instr = True2516 2514 continue 2517 2515 if diff.strip() == 'structural': 2518 instr = False2519 2516 struct = True 2520 2517 continue -
trunk/GSASIIconstrGUI.py
r2544 r2546 18 18 import wx 19 19 import wx.grid as wg 20 import wx.lib.scrolledpanel as wxscroll21 import time22 20 import random as ran 23 21 import numpy as np … … 1088 1086 else: # called directly, get current page 1089 1087 page = G2frame.dataDisplay.GetSelection() 1090 oldPage =G2frame.dataDisplay.ChangeSelection(page)1088 G2frame.dataDisplay.ChangeSelection(page) 1091 1089 text = G2frame.dataDisplay.GetPageText(page) 1092 1090 G2frame.dataFrame.ConstraintEdit.Enable(G2gd.wxID_EQUIVALANCEATOMS,False) … … 1210 1208 opId = oldPhase['pId'] 1211 1209 npId = newPhase['pId'] 1212 oph = '%d::'%(opId)1213 nph = '%d::'%(npId)1214 1210 cx,ct,cs,cia = newPhase['General']['AtomPtrs'] 1215 1211 nAtoms = newPhase['Atoms'] 1216 oAtoms = oldPhase['Atoms']1217 1212 oSGData = oldPhase['General']['SGData'] 1218 1213 nSGData = newPhase['General']['SGData'] 1219 1214 oAcof = G2lat.cell2A(oldPhase['General']['Cell'][1:7]) 1220 1215 nAcof = G2lat.cell2A(newPhase['General']['Cell'][1:7]) 1221 oGmat = G2lat.cell2Gmat(oldPhase['General']['Cell'][1:7])[1]1222 nGmat = G2lat.cell2Gmat(newPhase['General']['Cell'][1:7])[1]1223 1216 item = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,'Constraints') 1224 1217 constraints = G2frame.PatternTree.GetItemPyData(item) … … 1331 1324 else: 1332 1325 page = G2frame.dataDisplay.GetSelection() 1333 oldPage =G2frame.dataDisplay.ChangeSelection(page)1326 G2frame.dataDisplay.ChangeSelection(page) 1334 1327 text = G2frame.dataDisplay.GetPageText(page) 1335 1328 if text == 'Vector rigid bodies': … … 1367 1360 1368 1361 def getTextFile(): 1369 defDir = os.path.join(os.path.split(__file__)[0],'GSASIImacros')1370 1362 dlg = wx.FileDialog(G2frame,'Choose rigid body text file', '.', '', 1371 1363 "GSAS-II text file (*.txt)|*.txt|XYZ file (*.xyz)|*.xyz|" … … 1882 1874 try: 1883 1875 val = float(vecGrid.GetCellValue(r,c)) 1884 rbData['rb Vect'][imag][r][c] = val1876 rbData['rbXYZ'][r][c] = val 1885 1877 except ValueError: 1886 1878 pass -
trunk/GSASIIctrls.py
r2541 r2546 35 35 import GSASIIpath 36 36 GSASIIpath.SetVersionNumber("$Revision$") 37 import GSASIIgrid as G2gd 37 38 # import GSASIImath as G2mth 38 39 # import GSASIIIO as G2IO … … 681 682 try: 682 683 val = self.typ(tc.GetValue()) 683 except (ValueError, SyntaxError) as e:684 except (ValueError, SyntaxError): 684 685 if self.typ is float: # for float values, see if an expression can be evaluated 685 686 val = G2py3.FormulaEval(tc.GetValue()) … … 1621 1622 else: 1622 1623 refdictlst = None 1623 Id = G etPatternTreeItemId(G2frame,G2frame.root,hst)1624 hstData = G2frame.PatternTree.GetItemPyData(G etPatternTreeItemId(G2frame,Id,'Instrument Parameters'))[0]1624 Id = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,hst) 1625 hstData = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Instrument Parameters'))[0] 1625 1626 for h in copyList: 1626 Id = G etPatternTreeItemId(G2frame,G2frame.root,h)1627 instData = G2frame.PatternTree.GetItemPyData(G etPatternTreeItemId(G2frame,Id,'Instrument Parameters'))[0]1627 Id = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,h) 1628 instData = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Instrument Parameters'))[0] 1628 1629 if len(hstData) != len(instData) or hstData['Type'][0] != instData['Type'][0]: #don't mix data types or lam & lam1/lam2 parms! 1629 1630 print h+' not copied - instrument parameters not commensurate' 1630 1631 continue 1631 hData = G2frame.PatternTree.GetItemPyData(G etPatternTreeItemId(G2frame,Id,TreeItemType))1632 hData = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,TreeItemType)) 1632 1633 if TreeItemType == 'Instrument Parameters': 1633 1634 hData = hData[0] … … 2471 2472 dlg.SetSizer(sizer) 2472 2473 sizer.Fit(dlg) 2473 val =dlg.ShowModal()2474 dlg.ShowModal() 2474 2475 2475 2476 ################################################################################ … … 2953 2954 self._tc.SetInsertionPointEnd() 2954 2955 2955 def Clone(self ):2956 def Clone(self,grid): 2956 2957 return GridFractionEditor(grid) 2957 2958 … … 3819 3820 def ShowHelp(helpType,frame): 3820 3821 '''Called to bring up a web page for documentation.''' 3821 global htmlFirstUse 3822 global htmlFirstUse,htmlPanel,htmlFrame 3822 3823 # look up a definition for help info from dict 3823 3824 #helplink = helpLocDict.get(helpType) … … 3861 3862 '''Called to show a tutorial web page. 3862 3863 ''' 3863 global htmlFirstUse 3864 global htmlFirstUse,htmlPanel,htmlFrame 3864 3865 # determine if a web browser or the internal viewer should be used for help info 3865 3866 if GSASIIpath.GetConfigValue('Help_mode'): … … 4150 4151 GSASIIpath.svnUpdateDir(os.path.join(self.tutorialPath,i[0])) 4151 4152 else: 4152 fullpath = os.path.join(self.tutorialPath,i[0],i[1])4153 4153 fulldir = os.path.join(self.tutorialPath,i[0]) 4154 4154 URL = G2BaseURL+'/Tutorials/'+i[0]+'/' … … 4212 4212 print "Cancel" 4213 4213 dlg.Destroy() 4214 import sys;sys.exit()4214 sys.exit() 4215 4215 #====================================================================== 4216 4216 # test ScrolledMultiEditor -
trunk/GSASIIddataGUI.py
r2495 r2546 22 22 import GSASIIlattice as G2lat 23 23 import GSASIIspc as G2spc 24 import GSASIIElem as G2elem25 import GSASIIElemGUI as G2elemGUI26 24 import GSASIIplot as G2plt 27 25 import GSASIIgrid as G2gd 28 import GSASIIIO as G2IO29 import GSASIImath as G2mth30 26 import GSASIIpwd as G2pwd 31 27 import GSASIIphsGUI as G2phsGUI … … 293 289 def OnHstrainVal(event): 294 290 event.Skip() 295 Snames = G2spc.HStrainNames(SGData)296 291 Obj = event.GetEventObject() 297 292 hist,pid = Indx[Obj.GetId()] … … 987 982 h,k,l = hkl 988 983 Obj.SetValue('%3d %3d %3d'%(h,k,l)) 989 990 try: 991 histData = UseList[G2frame.hist] 992 except KeyError: 984 985 if G2frame.hist not in UseList: 993 986 G2frame.ErrorDialog('Missing data error', 994 987 G2frame.hist+' not in GSAS-II data tree') -
trunk/GSASIIexprGUI.py
r2326 r2546 24 24 ''' 25 25 import re 26 import sys27 26 import wx 28 import os.path29 27 import wx.lib.scrolledpanel as wxscroll 30 28 import numpy as np … … 400 398 try: 401 399 float(val) 402 except ValueError,TypeError:400 except (ValueError,TypeError): 403 401 invalid += 1 404 402 if msg: msg += "; " … … 953 951 newobj = dlg.Show(True) 954 952 print dlg.GetDepVar() 955 import sys956 #sys.exit()957 953 958 954 #app.MainLoop() -
trunk/GSASIIgrid.py
r2544 r2546 20 20 import time 21 21 import copy 22 import cPickle23 22 import sys 24 23 import os 25 24 import random as ran 26 25 import numpy as np 27 import numpy.ma as ma28 26 import scipy.optimize as so 29 27 import GSASIIpath … … 252 250 for line in text: 253 251 mainSizer.Add(wx.StaticText(self.panel,label=' %s '%(line)),0,WACV) 254 try:255 nops = self.names.index(' 1bar ')256 except ValueError:257 nops = len(self.names)258 252 ncol = self.table[0].count(',')+2 259 253 for ic,cent in enumerate(cents): … … 2702 2696 data['F**2'] = fsqRef.GetValue() 2703 2697 2704 def OnHatomFix(event):2705 data['HatomFix'] = Hfix.GetValue()2698 # def OnHatomFix(event): 2699 # data['HatomFix'] = Hfix.GetValue() 2706 2700 2707 2701 def OnUsrRej(event): … … 3379 3373 UseFlags = G2frame.SeqTable.GetColValues(0) 3380 3374 for obj in eqObjList: 3381 expr = obj.expression3382 3375 # assemble refined vars for this equation 3383 3376 varyValueDict.update({var:val for var,val in obj.GetVariedVarVal()}) … … 3490 3483 newobj = dlg.Show(True) 3491 3484 if newobj: 3492 calcobj = G2obj.ExpressionCalcObj(newobj)3493 3485 Controls['SeqParFitEqList'][selected] = newobj 3494 3486 EnableParFitEqMenus() … … 3545 3537 newobj = dlg.Show(True) 3546 3538 if newobj: 3547 calcobj = G2obj.ExpressionCalcObj(newobj)3548 3539 Controls['SeqParFitEqList'].append(newobj) 3549 3540 EnableParFitEqMenus() … … 3888 3879 G2mv.InitVars() 3889 3880 parmDict = data[name].get('parmDict') 3890 badVary = data[name].get('badVary',[])3891 3881 constraintInfo = data[name].get('constraintInfo',[[],[],{},[],seqnum]) 3892 3882 groups,parmlist,constrDict,fixedList,ihst = constraintInfo … … 3971 3961 G2frame.dataDisplay = G2G.GSGrid(parent=G2frame.dataFrame) 3972 3962 G2frame.SeqTable = G2G.Table( 3973 [list(c ) for cin zip(*colList)], # convert from columns to rows3963 [list(cl) for cl in zip(*colList)], # convert from columns to rows 3974 3964 colLabels=colLabels,rowLabels=histNames,types=Types) 3975 3965 G2frame.dataDisplay.SetTable(G2frame.SeqTable, True) … … 4020 4010 inp[2] = (inp[1] - inp[0])/(len(data[1][0])-1.) 4021 4011 names = ('start angle', 'end angle', 'step size') 4022 dictlst = [inp] * len(inp)4023 elemlst = range(len(inp))4024 4012 dlg = G2G.ScrolledMultiEditor( 4025 4013 G2frame,[inp] * len(inp), range(len(inp)), names, … … 4116 4104 Comments.append(" Transformation M*H = H' applied; M=") 4117 4105 Comments.append(str(Trans)) 4118 SG = 'P '+Laue4119 SGData = G2spc.SpcGroup(SG)[1]4120 4106 refList = G2lat.LaueUnique(Laue,refList) 4121 4107 dlg = wx.ProgressDialog('Build HKL dictonary','',len(refList)+1, … … 4197 4183 wtval.SetValue('%.3f'%(val)) 4198 4184 4199 def OnCompression(event):4200 data[0] = int(comp.GetValue())4185 # def OnCompression(event): 4186 # data[0] = int(comp.GetValue()) 4201 4187 4202 4188 def onCopyPlotCtrls(event): … … 4331 4317 tab = G2frame.G2plotNB.nb.GetPageText(page) 4332 4318 if '3D' in tab: 4333 Hmin = np.array([int(np.min(refList.T[0])),int(np.min(refList.T[1])),int(np.min(refList.T[2]))])4334 Hmax = np.array([int(np.max(refList.T[0])),int(np.max(refList.T[1])),int(np.max(refList.T[2]))])4335 Vpoint = [int(np.mean(refList.T[0])),int(np.mean(refList.T[1])),int(np.mean(refList.T[2]))]4336 4319 Page = G2frame.G2plotNB.nb.GetPage(page) 4337 4320 controls = Page.controls … … 4471 4454 data = G2frame.PatternTree.GetItemPyData(item) 4472 4455 Phases = G2frame.GetPhaseData() 4473 phase = ''4474 4456 phaseName = '' 4475 4457 if Phases: -
trunk/GSASIIimage.py
r2522 r2546 21 21 import numpy.linalg as nl 22 22 import numpy.ma as ma 23 import numpy.fft as fft24 import scipy.signal as signal25 23 import polymask as pm 26 24 from scipy.optimize import leastsq … … 32 30 import GSASIIpwd as G2pwd 33 31 import GSASIIspc as G2spc 34 import fellipse as fel35 32 36 33 # trig functions in degrees … … 91 88 ''' gives ellipse center coordinates 92 89 ''' 93 b,c,d,f, g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]90 b,c,d,f,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[0] 94 91 num = b*b-a*c 95 92 x0=(c*d-b*f)/num … … 101 98 range will be -90 to 90 deg 102 99 ''' 103 b,c, d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]100 b,c,a = p[1]/2, p[2], p[0] 104 101 return 0.5*npatand(2*b/(a-c)) 105 102 … … 170 167 phi0 = npatan2d(y-parms['det-Y'],x-parms['det-X']) 171 168 dxy = peneCorr(tth,parms['dep'],parms['tilt'],phi0) 172 ttth = nptand(tth)173 169 stth = npsind(tth) 174 170 cosb = npcosd(parms['tilt']) … … 176 172 tbm = nptand((tth-parms['tilt'])/2.) 177 173 tbp = nptand((tth+parms['tilt'])/2.) 178 sinb = npsind(parms['tilt'])179 174 d = parms['dist']+dxy 180 175 fplus = d*tanb*stth/(cosb+stth) … … 277 272 ''' 278 273 radii = [0,0] 279 ttth = tand(tth)280 274 stth = sind(tth) 281 ctth = cosd(tth)282 275 cosb = cosd(tilt) 283 276 tanb = tand(tilt) … … 332 325 cent = data['center'] 333 326 tth = 2.0*asind(data['wavelength']/(2.*dsp)) 334 ttth = tand(tth)335 327 stth = sind(tth) 336 ctth = cosd(tth)337 328 cosb = cosd(tilt) 338 329 if radii[0] > 0.: … … 680 671 Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,G2frame.ImageZ)[0] 681 672 ttth = nptand(tth) 682 stth = npsind(tth)683 673 ctth = npcosd(tth) 684 674 #1st estimate of tilt; assume ellipse - don't know sign though … … 879 869 numAzms = data['outAzimuths'] 880 870 numChans = data['outChannels'] 881 azmOff = data['azmthOff']882 871 Dazm = (LRazm[1]-LRazm[0])/numAzms 883 872 if '2-theta' in data.get('binType','2-theta'): … … 893 882 NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32) 894 883 H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32) 895 imageN = len(image)896 884 Nx,Ny = data['size'] 897 885 nXBlks = (Nx-1)/blkSize+1 -
trunk/GSASIIindex.py
r2500 r2546 20 20 import time 21 21 import numpy as np 22 import numpy.linalg as nl23 22 import GSASIIpath 24 23 GSASIIpath.SetVersionNumber("$Revision$") … … 191 190 if hkl[3] < d20: 192 191 break 193 eta = diff/Nobs20194 192 Q20 = 1.0/d20**2 195 193 if diff: … … 221 219 if hkl[4] < d20: 222 220 break 223 eta = diff/Nobs20224 221 Q20 = 1.0/d20**2 225 222 if diff: … … 1031 1028 def TestData(): 1032 1029 'needs a doc string' 1033 array = np.array1034 1030 global NeedTestData 1035 1031 NeedTestData = False … … 1099 1095 def test0(): 1100 1096 if NeedTestData: TestData() 1101 msg = 'test FitHKL'1102 1097 ibrav,cell,bestcell,Pwr,peaks = TestData 1103 1098 print 'best cell:',bestcell … … 1123 1118 def test1(): 1124 1119 if NeedTestData: TestData() 1125 msg = 'test FitHKL'1126 1120 ibrav,A,Pwr,peaks = TestData2 1127 1121 print 'bad cell:',G2lat.A2cell(A) -
trunk/GSASIIlattice.py
r2512 r2546 317 317 if Phase['General']['Type'] == 'magnetic': 318 318 cm = cx+4 319 unit = np.zeros(3)320 319 for iat,atom in enumerate(Atoms): 321 320 XYZ = np.array(atom[cx:cx+3]) … … 931 930 932 931 """ 933 import math934 932 if Bravais in [9,11]: 935 933 Cent = 'C' … … 1228 1226 matd3q = np.array([[0,0,-1],[-1,0,0],[0,1,0]]) #hkl -> -l,-h,k 1229 1227 matd3t = np.array([[0,0,-1],[1,0,0],[0,-1,0]]) #hkl -> -l,h,-k 1230 matd3p = np.array([[0,1,0],[0,0,-1],[-1,0,0]]) #hkl -> k,-l,-h1231 1228 mat6 = np.array([[1,1,0],[-1,0,0],[0,0,1]]) #hkl -> h+k,-h,l really 65 1232 1229 matdm = np.array([[0,1,0],[1,0,0],[0,0,1]]) #hkl -> k,h,l 1233 matdmt = np.array([[0,-1,0],[-1,0,0],[0,0,1]]) #hkl -> -k,-h,l1234 1230 matdmp = np.array([[-1,-1,0],[0,1,0],[0,0,1]]) #hkl -> -h-k,k,l 1235 matdmq = np.array([[-1,0,0],[1,1,0],[0,0,1]]) #hkl -> -h,h+k,l1236 1231 matkm = np.array([[-1,0,0],[1,1,0],[0,0,1]]) #hkl -> -h,h+k,l 1237 matkmp = np.array([[1,0,0],[-1,-1,0],[0,0,1]]) #hkl -> h,-h-k,l1238 1232 matd2 = np.array([[0,1,0],[1,0,0],[0,0,-1]]) #hkl -> k,h,-l 1239 matd2p = np.array([[-1,-1,0],[0,1,0],[0,0,-1]]) #hkl -> -h-k,k,-l1240 1233 matdm3 = np.array([[1,0,0],[0,0,1],[0,1,0]]) #hkl -> h,l,k 1241 1234 mat2d43 = np.array([[0,1,0],[1,0,0],[0,0,1]]) #hkl -> k,-h,l 1242 math2 = np.array([[0,-1,0],[-1,0,0],[0,0,-1]]) #hkl -> -k,-h,-l1243 1235 matk2 = np.array([[-1,0,0],[1,1,0],[0,0,-1]]) #hkl -> -h,-i,-l 1244 1236 #triclinic … … 2231 2223 spdict = spc.SpcGroup(key)[1] 2232 2224 cell = sgtbxlattinp.sgtbx8[key][0] 2233 center = spdict['SGLatt']2234 Laue = spdict['SGLaue']2235 2225 Axis = spdict['SGUniq'] 2236 2226 system = spdict['SGSys'] -
trunk/GSASIImapvars.py
r2519 r2546 210 210 211 211 import numpy as np 212 import sys 212 213 import GSASIIpath 213 214 GSASIIpath.SetVersionNumber("$Revision$") … … 238 239 fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations 239 240 consNum = 0 # number of the next constraint to be created 240 fixedVarList = []241 241 242 242 def VarKeys(constr): … … 816 816 if len(printlist) == 0: return 817 817 s1 = '' 818 s2 = '' 819 s3 = '' 818 820 print >>pFile,130*'-' 819 821 print >>pFile,"Variables generated by constraints" … … 845 847 for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList): 846 848 #if invmultarr is None: continue # probably not needed 847 try:848 valuelist = [parmDict[var] for var in mapvars]849 except KeyError:850 continue849 # try: 850 # valuelist = [parmDict[var] for var in mapvars] 851 # except KeyError: 852 # continue 851 853 # get the v-covar matrix for independent parameters 852 854 vcov = np.zeros((len(mapvars),len(mapvars))) … … 1114 1116 the first m rows are not linearly independent 1115 1117 ''' 1116 n = len(collist)1117 1118 for i in range(m): 1118 1119 if np.allclose(arr[i,i],0): -
trunk/GSASIImath.py
r2544 r2546 376 376 return XYZ 377 377 378 def TransformAtoms(Atoms,cx,cia,Trans,Vec):379 for Atom in Atoms:380 XYZ = Atom[cx:cx+3]381 if 'A' in Atom[cia]:382 U6 = Atom[cia+2:cia+8]378 #def TransformAtoms(Atoms,cx,cia,Trans,Vec): 379 # for Atom in Atoms: 380 # XYZ = Atom[cx:cx+3] 381 # if 'A' in Atom[cia]: 382 # U6 = Atom[cia+2:cia+8] 383 383 384 384 … … 564 564 return [],[] 565 565 566 def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj): #unfinished & not used567 '''default doc string568 569 :param type name: description570 571 :returns: type name: description572 573 '''574 for atom in atomData:575 XYZ = np.inner(Amat,atom[cx:cx+3])576 if atom[cia] == 'A':577 UIJ = atom[cia+2:cia+8]566 #def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj): #unfinished & not used 567 # '''default doc string 568 # 569 # :param type name: description 570 # 571 # :returns: type name: description 572 # 573 # ''' 574 # for atom in atomData: 575 # XYZ = np.inner(Amat,atom[cx:cx+3]) 576 # if atom[cia] == 'A': 577 # UIJ = atom[cia+2:cia+8] 578 578 579 579 def TLS2Uij(xyz,g,Amat,rbObj): #not used anywhere, but could be? … … 1159 1159 # GSASIIpath.IPyBreak() 1160 1160 if nWaves[0]: 1161 tauF = np.arange(1.,nWaves[0]+1 -nf)[:,nxs]*glTau #Fwaves x ngl1161 tauF = np.arange(1.,nWaves[0]+1)[:,nxs]*glTau #Fwaves x ngl 1162 1162 StauF = np.ones_like(Af)[:,:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf 1163 1163 CtauF = np.ones_like(Bf)[:,:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf … … 2898 2898 phasep = complex(a,b) 2899 2899 phasem = complex(a,-b) 2900 h,k,l = hkl+Hmax 2901 Ehkl[h,k,l] = E*phasep 2902 h,k,l = -hkl+Hmax #Friedel pair refl. 2903 Ehkl[h,k,l] = E*phasem 2900 Ehkl[hkl+Hmax] = E*phasep 2901 Ehkl[-hkl+Hmax] = E*phasem 2904 2902 # Ehkl[Hmax] = 0.00001 #this to preserve F[0,0,0] 2905 2903 testHKL = np.array(flipData['testHKL'])+Hmax … … 3515 3513 ins[x] = Parms[x][ind] 3516 3514 if ifQ: #qplot - convert back to 2-theta 3517 pos = 2.0*asind(pos* wave/(4*math.pi))3515 pos = 2.0*asind(pos*getWave(Parms)/(4*math.pi)) 3518 3516 sig = getCWsig(ins,pos) 3519 3517 gam = getCWgam(ins,pos) … … 4344 4342 lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False), 4345 4343 ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar) 4346 M =mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)4344 mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict) 4347 4345 # for ref in refs.T: 4348 4346 # print ' %4d %4d %4d %10.3f %10.3f %10.3f'%(int(ref[0]),int(ref[1]),int(ref[2]),ref[4],ref[5],ref[6]) -
trunk/GSASIIobj.py
r2504 r2546 1877 1877 try: 1878 1878 exprast = ast.parse(expr) 1879 except SyntaxError as err:1879 except SyntaxError: 1880 1880 s = '' 1881 1881 import traceback -
trunk/GSASIIphsGUI.py
r2544 r2546 4521 4521 SGData = generalData['SGData'] 4522 4522 SpnFlp = SGData.get('SpnFlp',[]) 4523 MagMom = SGData.get('MagMom',[])4523 # MagMom = SGData.get('MagMom',[]) 4524 4524 wx.BeginBusyCursor() 4525 4525 try: -
trunk/GSASIIplot.py
r2544 r2546 1952 1952 Parms = ParmList[N] 1953 1953 Sample = SampleList[N] 1954 if 'C' in Parms['Type'][0]:1955 wave = G2mth.getWave(Parms)1956 else:1957 difC = Parms['difC'][1]1958 1954 ifpicked = False 1959 1955 LimitId = 0 … … 2391 2387 PlotISFG(G2frame,newPlot=newPlot,type=type) 2392 2388 2393 def OnKeyBox(event):2394 if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index(type):2395 event.key = cb.GetValue()[0]2396 cb.SetValue(' key press')2397 wx.CallAfter(OnPlotKeyPress,event)2398 Page.canvas.SetFocus() # redirect the Focus from the button back to the plot2399 2400 2389 def OnMotion(event): 2401 2390 xpos = event.xdata … … 2772 2761 wx.CallAfter(PlotXYZ,G2frame,XY,Z,labelX,labelY,False,Title) 2773 2762 2774 def OnKeyBox(event):2775 if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index(type):2776 event.key = cb.GetValue()[0]2777 cb.SetValue(' key press')2778 wx.CallAfter(OnKeyPress,event)2779 Page.canvas.SetFocus() # redirect the Focus from the button back to the plot2780 2781 2763 def OnMotion(event): 2782 2764 xpos = event.xdata … … 3039 3021 Plot.set_xlabel(r'$Q, \AA^{-1}$',fontsize=14) 3040 3022 Plot.set_ylabel(r'$\Delta Q/Q, \Delta d/d$',fontsize=14) 3041 try: 3042 Xmin,Xmax = limits[1] 3043 X = np.linspace(Xmin,Xmax,num=101,endpoint=True) 3044 Q = 4.*np.pi*npsind(X/2.)/lam 3045 Z = np.ones_like(X) 3046 data = G2mth.setPeakparms(Parms,Parms2,X,Z) 3047 s = np.sqrt(data[4])*np.pi/18000. #var -> sig(radians) 3048 g = data[6]*np.pi/18000. #centideg -> radians 3049 G = G2pwd.getgamFW(g,s)/2. #delt-theta 3050 Y = s/nptand(X/2.) 3051 Z = g/nptand(X/2.) 3052 W = G/nptand(X/2.) 3053 Plot.plot(Q,Y,color='r',label='Gaussian') 3054 Plot.plot(Q,Z,color='g',label='Lorentzian') 3055 Plot.plot(Q,W,color='b',label='G+L') 3056 3057 fit = G2mth.setPeakparms(Parms,Parms2,X,Z,useFit=True) 3058 sf = np.sqrt(fit[4])*np.pi/18000. 3059 gf = fit[6]*np.pi/18000. 3060 Gf = G2pwd.getgamFW(gf,sf)/2. 3061 Yf = sf/nptand(X/2.) 3062 Zf = gf/nptand(X/2.) 3063 Wf = Gf/nptand(X/2.) 3064 Plot.plot(Q,Yf,color='r',dashes=(5,5),label='Gaussian fit') 3065 Plot.plot(Q,Zf,color='g',dashes=(5,5),label='Lorentzian fit') 3066 Plot.plot(Q,Wf,color='b',dashes=(5,5),label='G+L fit') 3067 3068 X = [] 3069 Y = [] 3070 Z = [] 3071 W = [] 3072 for peak in peaks: 3073 X.append(4.0*math.pi*sind(peak[0]/2.0)/lam) 3074 try: 3075 s = math.sqrt(peak[4])*math.pi/18000. 3076 except ValueError: 3077 s = 0.01 3078 g = peak[6]*math.pi/18000. 3079 G = G2pwd.getgamFW(g,s)/2. 3080 Y.append(s/tand(peak[0]/2.)) 3081 Z.append(g/tand(peak[0]/2.)) 3082 W.append(G/tand(peak[0]/2.)) 3083 if len(peaks): 3084 Plot.plot(X,Y,'+',color='r',label='G peak') 3085 Plot.plot(X,Z,'+',color='g',label='L peak') 3086 Plot.plot(X,W,'+',color='b',label='G+L peak') 3087 Plot.legend(loc='best') 3088 Page.canvas.draw() 3089 except ValueError: 3090 print '**** ERROR - default U,V,W profile coefficients yield sqrt of negative value at 2theta =', \ 3091 '%.3f'%(2*theta) 3092 G2frame.G2plotNB.Delete('Peak Widths') 3023 Xmin,Xmax = limits[1] 3024 X = np.linspace(Xmin,Xmax,num=101,endpoint=True) 3025 Q = 4.*np.pi*npsind(X/2.)/lam 3026 Z = np.ones_like(X) 3027 data = G2mth.setPeakparms(Parms,Parms2,X,Z) 3028 s = np.sqrt(data[4])*np.pi/18000. #var -> sig(radians) 3029 g = data[6]*np.pi/18000. #centideg -> radians 3030 G = G2pwd.getgamFW(g,s)/2. #delt-theta 3031 Y = s/nptand(X/2.) 3032 Z = g/nptand(X/2.) 3033 W = G/nptand(X/2.) 3034 Plot.plot(Q,Y,color='r',label='Gaussian') 3035 Plot.plot(Q,Z,color='g',label='Lorentzian') 3036 Plot.plot(Q,W,color='b',label='G+L') 3037 3038 fit = G2mth.setPeakparms(Parms,Parms2,X,Z,useFit=True) 3039 sf = np.sqrt(fit[4])*np.pi/18000. 3040 gf = fit[6]*np.pi/18000. 3041 Gf = G2pwd.getgamFW(gf,sf)/2. 3042 Yf = sf/nptand(X/2.) 3043 Zf = gf/nptand(X/2.) 3044 Wf = Gf/nptand(X/2.) 3045 Plot.plot(Q,Yf,color='r',dashes=(5,5),label='Gaussian fit') 3046 Plot.plot(Q,Zf,color='g',dashes=(5,5),label='Lorentzian fit') 3047 Plot.plot(Q,Wf,color='b',dashes=(5,5),label='G+L fit') 3048 3049 X = [] 3050 Y = [] 3051 Z = [] 3052 W = [] 3053 for peak in peaks: 3054 X.append(4.0*math.pi*sind(peak[0]/2.0)/lam) 3055 try: 3056 s = math.sqrt(peak[4])*math.pi/18000. 3057 except ValueError: 3058 s = 0.01 3059 g = peak[6]*math.pi/18000. 3060 G = G2pwd.getgamFW(g,s)/2. 3061 Y.append(s/tand(peak[0]/2.)) 3062 Z.append(g/tand(peak[0]/2.)) 3063 W.append(G/tand(peak[0]/2.)) 3064 if len(peaks): 3065 Plot.plot(X,Y,'+',color='r',label='G peak') 3066 Plot.plot(X,Z,'+',color='g',label='L peak') 3067 Plot.plot(X,W,'+',color='b',label='G+L peak') 3068 Plot.legend(loc='best') 3069 Page.canvas.draw() 3093 3070 else: #'T'OF 3094 3071 Plot.set_title('Instrument and sample peak coefficients') … … 3397 3374 except ValueError: 3398 3375 sfac *= 1.05 3399 Z = [lut(r *np.pi/180.,p*np.pi/180.) for r,p in zip(list(R),list(P))]3376 Z = [lut(ri*np.pi/180.,p*np.pi/180.) for ri,p in zip(list(R),list(P))] 3400 3377 print 'IVP for histogramn: %s: interpolate sfactor: %.2f'%(hist,sfac) 3401 3378 except AttributeError: … … 4241 4218 wx.CallAfter(PlotImage,G2frame,newPlot=True) 4242 4219 4243 def OnKeyBox(event):4244 if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index('2D Powder Image'):4245 event.key = cb.GetValue()[0]4246 cb.SetValue(' key press')4247 if event.key in ['l','s','a','r','p','x','y']:4248 wx.CallAfter(OnImPlotKeyPress,event)4249 Page.canvas.SetFocus() # redirect the Focus from the button back to the plot4250 4251 4220 def OnImPick(event): 4252 4221 if G2frame.PatternTree.GetItemText(G2frame.PickId) not in ['Image Controls','Masks']: … … 4862 4831 drawAtoms = drawingData.get('Atoms',[]) 4863 4832 mapData = {} 4864 flipData = {}4865 4833 showBonds = False 4866 4834 if 'Map' in generalData: 4867 4835 mapData = generalData['Map'] 4868 4836 showBonds = mapData.get('Show bonds',False) 4869 if 'Flip' in generalData:4870 flipData = generalData['Flip']4871 4837 Wt = np.array([255,255,255]) 4872 4838 Rd = np.array([255,0,0]) -
trunk/GSASIIpwd.py
r2522 r2546 20 20 21 21 import numpy as np 22 import scipy as sp23 22 import numpy.linalg as nl 24 23 import numpy.ma as ma … … 34 33 import GSASIIspc as G2spc 35 34 import GSASIIElem as G2elem 36 import GSASIIgrid as G2gd37 import GSASIIIO as G2IO38 35 import GSASIImath as G2mth 39 36 import pypowder as pyd … … 180 177 181 178 Sth2 = npsind(Tth/2.0)**2 182 Cth2 = 1.-Sth2183 179 if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample 184 180 if 'array' in str(type(MuR)): … … 308 304 #Apply angle dependent corrections 309 305 Tth = xydata['IofQ'][1][0] 310 dt = (Tth[1]-Tth[0])311 306 MuR = Abs*data['Diam']/20.0 312 307 xydata['IofQ'][1][1] /= Absorb(data['Geometry'],MuR,Tth) … … 318 313 #convert to Q 319 314 if 'C' in inst['Type'][0]: 320 hc = 12.397639321 315 wave = G2mth.getWave(inst) 322 keV = hc/wave323 316 minQ = npT2q(Tth[0],wave) 324 317 maxQ = npT2q(Tth[-1],wave) … … 559 552 extended by exponential coeff. 560 553 ''' 561 lnf = 3.3 # =log(0.001/2)562 554 widths = [np.sqrt(sig),gam] 563 555 fwhm = 2.355*widths[0]+2.*widths[1] … … 887 879 Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl) 888 880 # Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl) 889 sumDf = np.sum(Df)890 881 return Df,dFdp,dFds,dFdg,dFdsh 891 882 … … 901 892 902 893 Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam) 903 sumDf = np.sum(Df)904 894 return Df,dFdp,dFds,dFdg 905 895 … … 913 903 'needs a doc string' 914 904 Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam) 915 sumDf = np.sum(Df)916 905 return Df,dFdp,dFda,dFdb,dFds,dFdg 917 906 … … 1688 1677 if controls: 1689 1678 Ftol = controls['min dM/M'] 1690 derivType = controls['deriv type']1691 1679 else: 1692 1680 Ftol = 0.0001 1693 derivType = 'analytic'1694 1681 if oneCycle: 1695 1682 Ftol = 1.0 … … 1825 1812 Parms = [] 1826 1813 #cell parms 1827 cell = Layers['Cell']1828 1814 if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']: 1829 1815 Parms.append('cellA') … … 1836 1822 Parms.append('cellG') 1837 1823 #Transition parms 1838 Trans = Layers['Transitions']1839 1824 for iY in range(len(Layers['Layers'])): 1840 1825 for iX in range(len(Layers['Layers'])): … … 2044 2029 SFdat = [] 2045 2030 for atType in atTypes: 2046 if atType == 'H':2047 blen = -.37412048 else:2049 blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]2050 2031 Adat = atmdata.XrayFF[atType] 2051 2032 SF = np.zeros(9) … … 2206 2187 SetStackingTrans(Layers,Nlayers) 2207 2188 # result as Sadp 2208 mirror = laueId in [-1,2,3,7,8,9,10]2209 2189 Nspec = 20001 2210 2190 spec = np.zeros(Nspec,dtype='double') … … 2237 2217 'needs a doc string' 2238 2218 # global NeedTestData 2239 NeedTestData = False2240 2219 global bakType 2241 2220 bakType = 'chebyschev' … … 2275 2254 def test0(): 2276 2255 if NeedTestData: TestData() 2277 msg = 'test '2278 2256 gplot = plotter.add('FCJ-Voigt, 11BM').gca() 2279 2257 gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0]) … … 2287 2265 time0 = time.time() 2288 2266 for i in range(100): 2289 y =getPeakProfile(parmDict1,xdata,varyList,bakType)2267 getPeakProfile(parmDict1,xdata,varyList,bakType) 2290 2268 print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0 2291 2269 -
trunk/GSASIIpwdGUI.py
r2544 r2546 24 24 import numpy.ma as ma 25 25 import math 26 import time27 26 import copy 28 27 import random as ran … … 557 556 prevVaryList = [] 558 557 Names = [] 558 peaks = None 559 varyList = None 559 560 if Reverse: 560 561 names.reverse() … … 1238 1239 G2frame.dataFrame.SetLabel('Background') 1239 1240 if not G2frame.dataFrame.GetStatusBar(): 1240 Status =G2frame.dataFrame.CreateStatusBar()1241 G2frame.dataFrame.CreateStatusBar() 1241 1242 G2frame.Bind(wx.EVT_MENU,OnBackCopy,id=G2gd.wxID_BACKCOPY) 1242 1243 G2frame.Bind(wx.EVT_MENU,OnBackFlagCopy,id=G2gd.wxID_BACKFLAGCOPY) … … 1434 1435 G2frame.dataFrame.SetLabel('Limits') 1435 1436 if not G2frame.dataFrame.GetStatusBar(): 1436 Status =G2frame.dataFrame.CreateStatusBar()1437 G2frame.dataFrame.CreateStatusBar() 1437 1438 G2frame.Bind(wx.EVT_MENU,OnLimitCopy,id=G2gd.wxID_LIMITCOPY) 1438 1439 G2frame.Bind(wx.EVT_MENU,OnAddExcl,id=G2gd.wxID_ADDEXCLREGION) … … 1924 1925 continue 1925 1926 nDig = (10,3) 1926 fmt = '%10.3f'1927 1927 if 'beta' in item: 1928 fmt = '%12.6g'1929 1928 nDig = (12,6) 1930 1929 instSizer.Add( … … 2044 2043 Inst2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame, 2045 2044 G2frame.PatternId,'Instrument Parameters'))[1] 2046 try:2047 histoName = G2frame.PatternTree.GetItemPyData(G2frame.PatternId)[-1]2048 ifHisto = IsHistogramInAnyPhase(G2frame,histoName)2049 except TypeError: #PKS data never used in a phase as data2050 ifhisto = False2051 2045 G2gd.SetDataMenuBar(G2frame) 2052 2046 #patch … … 2419 2413 G2frame.dataFrame.SetScale.Enable(True) 2420 2414 if not G2frame.dataFrame.GetStatusBar(): 2421 Status =G2frame.dataFrame.CreateStatusBar()2415 G2frame.dataFrame.CreateStatusBar() 2422 2416 G2frame.dataDisplay = wx.Panel(G2frame.dataFrame) 2423 2417 Controls = G2frame.PatternTree.GetItemPyData( … … 2643 2637 G2frame.dataFrame.Clear() 2644 2638 if not G2frame.dataFrame.GetStatusBar(): 2645 Status =G2frame.dataFrame.CreateStatusBar()2639 G2frame.dataFrame.CreateStatusBar() 2646 2640 if 'PWD' in G2frame.PatternTree.GetItemText(G2frame.PatternId): 2647 2641 G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.IndPeaksMenu) … … 3534 3528 General = G2frame.PatternTree.GetItemPyData(phaseId)['General'] 3535 3529 Super = General.get('Super',0) 3536 SuperVec = General.get('SuperVec',[])3537 3530 else: 3538 3531 Super = 0 … … 3988 3981 G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.SubstanceMenu) 3989 3982 if not G2frame.dataFrame.GetStatusBar(): 3990 Status =G2frame.dataFrame.CreateStatusBar()3983 G2frame.dataFrame.CreateStatusBar() 3991 3984 G2frame.dataDisplay = wxscroll.ScrolledPanel(G2frame.dataFrame) 3992 3985 G2frame.dataFrame.SetLabel('Substances') … … 4162 4155 if Reverse: 4163 4156 names.reverse() 4157 JModel = None 4164 4158 try: 4165 4159 for i,name in enumerate(names): … … 4798 4792 Obj.SetValue(fmt%(value)) 4799 4793 data[fileKey][itemKey] = value 4800 auxPlot = ComputePDF(data) 4801 G2plt.PlotISFG(G2frame,newPlot=True) 4794 OnComputePDF(None) 4802 4795 4803 4796 item = data[key] … … 4829 4822 sumVol += Avol*ElList[El]['FormulaNo'] 4830 4823 return sumVol 4831 auxPlot = ComputePDF(data) 4832 G2plt.PlotISFG(G2frame,newPlot=True) 4824 OnComputePDF(None) 4833 4825 4834 4826 def FillElemSizer(elemSizer,ElData): … … 4845 4837 formVol.SetValue('%.2f'%(data['Form Vol'])) 4846 4838 wx.CallAfter(UpdatePDFGrid,G2frame,data) 4847 auxPlot = ComputePDF(data) 4848 G2plt.PlotISFG(G2frame,newPlot=True) 4839 OnComputePDF(None) 4849 4840 4850 4841 elemSizer.Add(wx.StaticText(parent=G2frame.dataDisplay, … … 4861 4852 data['Geometry'] = geometry.GetValue() 4862 4853 UpdatePDFGrid(G2frame,data) 4863 auxPlot = ComputePDF(data) 4864 G2plt.PlotISFG(G2frame,newPlot=True) 4854 OnComputePDF(None) 4865 4855 4866 4856 def OnDetType(event): 4867 4857 data['DetType'] = detType.GetValue() 4868 4858 UpdatePDFGrid(G2frame,data) 4869 auxPlot = ComputePDF(data) 4870 G2plt.PlotISFG(G2frame,newPlot=True) 4859 OnComputePDF(None) 4871 4860 4872 4861 def OnFormVol(event): … … 4880 4869 data['Form Vol'] = value 4881 4870 UpdatePDFGrid(G2frame,data) 4882 auxPlot = ComputePDF(data) 4883 G2plt.PlotISFG(G2frame,newPlot=False) 4871 OnComputePDF(None) 4884 4872 4885 4873 def OnDiameter(event): … … 4893 4881 data['Diam'] = value 4894 4882 UpdatePDFGrid(G2frame,data) 4895 auxPlot = ComputePDF(data) 4896 G2plt.PlotISFG(G2frame,newPlot=False) 4883 OnComputePDF(None) 4897 4884 4898 4885 # def OnPolaVal(event): … … 4936 4923 data['ObliqCoeff'] = value 4937 4924 obliqCoeff.SetValue('%.3f'%(value)) 4938 auxPlot = ComputePDF(data) 4939 G2plt.PlotISFG(G2frame,newPlot=False) 4925 OnComputePDF(None) 4940 4926 4941 4927 def OnBackVal(event): … … 4948 4934 data['BackRatio'] = value 4949 4935 backVal.SetValue('%.3f'%(value)) 4950 auxPlot = ComputePDF(data) 4951 G2plt.PlotISFG(G2frame,newPlot=False) 4936 OnComputePDF(None) 4952 4937 4953 4938 def OnBackSlider(event): … … 4955 4940 data['BackRatio'] = value 4956 4941 backVal.SetValue('%.3f'%(data['BackRatio'])) 4957 auxPlot = ComputePDF(data) 4958 G2plt.PlotISFG(G2frame,newPlot=False) 4942 OnComputePDF(None) 4959 4943 4960 4944 def OnRulandWdt(event): … … 4970 4954 data['Ruland'] = value 4971 4955 rulandWdt.SetValue('%.3f'%(value)) 4972 auxPlot = ComputePDF(data) 4973 G2plt.PlotISFG(G2frame,newPlot=False) 4956 OnComputePDF(None) 4974 4957 4975 4958 def OnRulSlider(event): … … 4977 4960 data['Ruland'] = max(0.001,value) 4978 4961 rulandWdt.SetValue('%.3f'%(data['Ruland'])) 4979 auxPlot = ComputePDF(data) 4980 G2plt.PlotISFG(G2frame,newPlot=False) 4962 OnComputePDF(None) 4981 4963 4982 4964 def OnLorch(event): 4983 4965 data['Lorch'] = lorch.GetValue() 4984 auxPlot = ComputePDF(data) 4985 G2plt.PlotISFG(G2frame,newPlot=False) 4966 OnComputePDF(None) 4986 4967 4987 4968 def OnPacking(event): … … 4995 4976 data['Pack'] = value 4996 4977 UpdatePDFGrid(G2frame,data) 4997 auxPlot = ComputePDF(data) 4998 G2plt.PlotISFG(G2frame,newPlot=False) 4978 OnComputePDF(None) 4999 4979 5000 4980 def OnSQmin(event): … … 5008 4988 data['QScaleLim'][0] = value 5009 4989 SQmin.SetValue('%.1f'%(value)) 5010 auxPlot = ComputePDF(data) 5011 G2plt.PlotISFG(G2frame,newPlot=True) 4990 OnComputePDF(None) 5012 4991 5013 4992 def OnSQmax(event): … … 5024 5003 SQmin.SetValue('%.1f'%(data['QScaleLim'][0])) 5025 5004 SQmax.SetValue('%.1f'%(value)) 5026 auxPlot = ComputePDF(data) 5027 G2plt.PlotISFG(G2frame,newPlot=True) 5005 OnComputePDF(None) 5028 5006 5029 5007 def OnResetQ(event): … … 5033 5011 data['QScaleLim'][0] = 0.9*qLimits[1] 5034 5012 SQmin.SetValue('%.1f'%(data['QScaleLim'][0])) 5035 auxPlot = ComputePDF(data) 5036 G2plt.PlotISFG(G2frame,newPlot=True) 5013 OnComputePDF(None) 5037 5014 5038 5015 def OnNoRing(event): 5039 5016 data['noRing'] = not data['noRing'] 5040 auxPlot = ComputePDF(data) 5041 G2plt.PlotISFG(G2frame,newPlot=True) 5042 5017 OnComputePDF(None) 5043 5018 5044 5019 def GetFileList(fileType): … … 5132 5107 5133 5108 def OnComputePDF(event): 5134 print 'Calculating PDF:'5109 # print 'Calculating PDF:' 5135 5110 auxPlot = ComputePDF(data) 5136 print 'Done calculating PDF:'5111 # print 'Done calculating PDF:' 5137 5112 Status.SetStatusText('PDF computed') 5138 5113 for plot in auxPlot: 5139 5114 XY = np.array(plot[:2]) 5140 5115 G2plt.PlotXY(G2frame,[XY,],Title=plot[2]) 5141 5142 G2plt.PlotISFG(G2frame,newPlot=True,type='I(Q)')5143 G2plt.PlotISFG(G2frame,newPlot=True,type='S(Q)')5144 G2plt.PlotISFG(G2frame,newPlot=True,type='F(Q)')5116 if event is not None: 5117 G2plt.PlotISFG(G2frame,newPlot=True,type='I(Q)') 5118 G2plt.PlotISFG(G2frame,newPlot=True,type='S(Q)') 5119 G2plt.PlotISFG(G2frame,newPlot=True,type='F(Q)') 5145 5120 G2plt.PlotISFG(G2frame,newPlot=True,type='G(R)') 5146 5121 … … 5153 5128 if 'PDF' in Name.split()[0]: 5154 5129 Data = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,id,'PDF Controls')) 5155 auxPlot =ComputePDF(Data)5130 ComputePDF(Data) 5156 5131 id, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie) 5157 5132 Status.SetStatusText('All PDFs computed') -
trunk/GSASIIpy3.py
r2065 r2546 111 111 try: 112 112 return fmt.format(float(val)).strip() 113 except ValueError as err:113 except ValueError: 114 114 print 'FormatValue Error with val,maxdigits,fmt=',val,maxdigits,fmt 115 115 return str(val) … … 162 162 try: 163 163 return fmt.format(float(val)).strip() 164 except ValueError as err:164 except ValueError: 165 165 print 'FormatValue Error with val,maxdigits, sigfigs, fmt=',val, maxdigits,sigfigs, fmt 166 166 return str(val) -
trunk/GSASIIrestrGUI.py
r2530 r2546 17 17 import wx 18 18 import wx.grid as wg 19 import time20 19 import numpy as np 21 20 import numpy.ma as ma … … 129 128 print '**** ERROR - wrong restraint macro file selected, try again ****' 130 129 macro = [] 131 else: # cancel was pressed132 macxro = []133 130 finally: 134 131 dlg.Destroy() … … 280 277 281 278 def AddAABondRestraint(bondRestData): 282 Radii = dict(zip(General['AtomTypes'],General['BondRadii']))283 279 macro = getMacroFile('bond') 284 280 if not macro: … … 389 385 ids = [vecta[4],vecta[0],vectb[4]] 390 386 ops = [vecta[5],vecta[1],vectb[5]] 391 XYZ = np.array([vecta[6],vecta[2],vectb[6]])392 387 angle = [ids,ops,value,1.0] 393 388 if angle not in angleRestData['Angles']: … … 551 546 def makeChains(Names,Ids): 552 547 Chains = {} 553 chain = ''554 548 atoms = zip(Names,Ids) 555 549 for name,id in atoms: … … 837 831 def OnCellChange(event): 838 832 r,c = event.GetRow(),event.GetCol() 839 val = bondRestData['Bonds'][r][c]840 833 try: 841 834 new = float(bondTable.GetValue(r,c)) … … 889 882 if BondRestr.GetSizer(): 890 883 BondRestr.GetSizer().Clear(True) 891 dataDisplay = wx.Panel(BondRestr)892 884 mainSizer = wx.BoxSizer(wx.VERTICAL) 893 885 mainSizer.Add((5,5),0) … … 969 961 def OnCellChange(event): 970 962 r,c = event.GetRow(),event.GetCol() 971 val = angleRestData['Angles'][r][c]972 963 try: 973 964 new = float(angleTable.GetValue(r,c)) … … 1020 1011 if AngleRestr.GetSizer(): 1021 1012 AngleRestr.GetSizer().Clear(True) 1022 dataDisplay = wx.Panel(AngleRestr)1023 1013 mainSizer = wx.BoxSizer(wx.VERTICAL) 1024 1014 mainSizer.Add((5,5),0) … … 1105 1095 def OnCellChange(event): 1106 1096 r,c = event.GetRow(),event.GetCol() 1107 val = planeRestData['Planes'][r][c]1108 1097 try: 1109 1098 new = float(planeTable.GetValue(r,c)) … … 1142 1131 if PlaneRestr.GetSizer(): 1143 1132 PlaneRestr.GetSizer().Clear(True) 1144 dataDisplay = wx.Panel(PlaneRestr)1145 1133 mainSizer = wx.BoxSizer(wx.VERTICAL) 1146 1134 mainSizer.Add((5,5),0) … … 1227 1215 def OnCellChange(event): 1228 1216 r,c = event.GetRow(),event.GetCol() 1229 val = chiralRestData['Volumes'][r][c]1230 1217 try: 1231 1218 new = float(volumeTable.GetValue(r,c)) … … 1278 1265 if ChiralRestr.GetSizer(): 1279 1266 ChiralRestr.GetSizer().Clear(True) 1280 dataDisplay = wx.Panel(ChiralRestr)1281 1267 mainSizer = wx.BoxSizer(wx.VERTICAL) 1282 1268 mainSizer.Add((5,5),0) … … 1358 1344 def OnCellChange(event): 1359 1345 r,c = event.GetRow(),event.GetCol() 1360 val = torsionRestData['Torsions'][r][c]1361 1346 try: 1362 1347 new = float(torsionTable.GetValue(r,c)) … … 1369 1354 1370 1355 def OnDeleteRestraint(event): 1371 rows = GetSelectedRows(Torsion s)1356 rows = GetSelectedRows(TorsionRestr.Torsions) 1372 1357 if not rows: 1373 1358 return … … 1379 1364 1380 1365 def OnChangeEsd(event): 1381 rows = GetSelectedRows(Torsion s)1366 rows = GetSelectedRows(TorsionRestr.Torsions) 1382 1367 if not rows: 1383 1368 return 1384 Torsion s.ClearSelection()1369 TorsionRestr.Torsions.ClearSelection() 1385 1370 val = torsionList[rows[0]][4] 1386 1371 dlg = G2G.SingleFloatDialog(G2frame,'New value','Enter new esd for torsion restraints',val,[0.,5.],'%.2f') … … 1395 1380 if TorsionRestr.GetSizer(): 1396 1381 TorsionRestr.GetSizer().Clear(True) 1397 dataDisplay = wx.Panel(TorsionRestr)1398 1382 mainSizer = wx.BoxSizer(wx.VERTICAL) 1399 1383 mainSizer.Add((5,5),0) … … 1479 1463 def OnCellChange(event): 1480 1464 r,c = event.GetRow(),event.GetCol() 1481 val = ramaRestData['Ramas'][r][c]1482 1465 try: 1483 1466 new = float(ramaTable.GetValue(r,c)) … … 1490 1473 1491 1474 def OnDeleteRestraint(event): 1492 rows = GetSelectedRows(Rama s)1475 rows = GetSelectedRows(RamaRestr.Ramas) 1493 1476 if not rows: 1494 1477 return … … 1500 1483 1501 1484 def OnChangeEsd(event): 1502 rows = GetSelectedRows(Rama s)1485 rows = GetSelectedRows(RamaRestr.Ramas) 1503 1486 if not rows: 1504 1487 return 1505 Rama s.ClearSelection()1488 RamaRestr.Ramas.ClearSelection() 1506 1489 val = ramaList[rows[0]][4] 1507 1490 dlg = G2G.SingleFloatDialog(G2frame,'New value','Enter new esd for energy',val,[0.,5.],'%.2f') … … 1516 1499 if RamaRestr.GetSizer(): 1517 1500 RamaRestr.GetSizer().Clear(True) 1518 dataDisplay = wx.Panel(RamaRestr)1519 1501 mainSizer = wx.BoxSizer(wx.VERTICAL) 1520 1502 mainSizer.Add((5,5),0) … … 1654 1636 if ChemCompRestr.GetSizer(): 1655 1637 ChemCompRestr.GetSizer().Clear(True) 1656 dataDisplay = wx.Panel(ChemCompRestr)1657 1638 mainSizer = wx.BoxSizer(wx.VERTICAL) 1658 1639 mainSizer.Add((5,5),0) … … 1745 1726 def OnCellChange(event): 1746 1727 r,c = event.GetRow(),event.GetCol() 1747 val = textureRestData['HKLs'][r][c]1748 1728 try: 1749 1729 if c == 1: #grid size 1750 1730 new = int(textureTable.GetValue(r,c)) 1751 1731 if new < 6 or new > 24: 1752 raise valueError1732 raise ValueError 1753 1733 elif c in [2,4]: #esds 1754 1734 new = float(textureTable.GetValue(r,c)) … … 1765 1745 if TextureRestr.GetSizer(): 1766 1746 TextureRestr.GetSizer().Clear(True) 1767 dataDisplay = wx.Panel(TextureRestr)1768 1747 mainSizer = wx.BoxSizer(wx.VERTICAL) 1769 1748 mainSizer.Add((5,5),0) … … 1775 1754 1776 1755 textureList = textureRestData['HKLs'] 1777 textureData = General['SH Texture']1778 1756 if len(textureList): 1779 1757 table = [] -
trunk/GSASIIsasd.py
r2308 r2546 564 564 lam2 = (-1.0*qb+np.sqrt(radic))/(2.0*qa) 565 565 lam = min(lam1,lam2) 566 test = 1.0 + 2.0*eta567 566 mu = lam*eta*etam1 568 567 alp = (1.0 + 2.0*eta - mu)/(etam1*etam1) … … 939 938 print " MaxEnt trial/max: %3d/%3d" % ((iter+1), IterMax) 940 939 print " Residual: %5.2lf%% Entropy: %8lg" % (100*test, S) 941 if iter > 0:942 value = 100*( math.sqrt(chisq/chtarg)-1)943 else:944 value = 0945 # print " %12.5lg %10.4lf" % ( math.sqrt(chtarg/npt), value )946 940 print " Function sum: %.6lg Change from last: %.2lf%%\n" % (fSum,100*fChange/fSum) 947 941 … … 1057 1051 'Spherical shell': [SphericalShellFF, SphericalShellVol]} 1058 1052 Shape = data['Size']['Shape'][0] 1059 SlitLen = Sample.get('SlitLen',0.0)1060 1053 Parms = data['Size']['Shape'][1:] 1061 1054 if data['Size']['logBins']: … … 1092 1085 Ic[Ibeg:Ifin] += Back[0] 1093 1086 print ' Final chi^2: %.3f'%(chisq) 1094 Vols = shapes[Shape][1](Bins,Parms)1095 1087 data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)] 1096 1088 … … 1269 1261 result = so.leastsq(calcSASD,values,full_output=True,epsfcn=1.e-8, #ftol=Ftol, 1270 1262 args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],levelTypes,parmDict,varyList)) 1271 ncyc = int(result[2]['nfev']/2)1272 1263 parmDict.update(zip(varyList,result[0])) 1273 1264 chisq = np.sum(result[2]['fvec']**2) … … 1323 1314 # pdb.set_trace() 1324 1315 partData = sasdData['Particle'] 1325 matFrac = partData['Matrix']['VolFrac'] #[value,flag]1326 Scale = Sample['Scale'] #[value,flag]1327 SlitLen = Sample.get('SlitLen',0.0)1328 1316 Back = sasdData['Back'] 1329 1317 Q,Io,wt,Ic,Ib,Ifb = Profile[:6] 1330 1318 Qmin = Limits[1][0] 1331 1319 Qmax = Limits[1][1] 1332 wtFactor = ProfDict['wtFactor']1333 1320 Ibeg = np.searchsorted(Q,Qmin) 1334 1321 Ifin = np.searchsorted(Q,Qmax)+1 #include last point … … 1416 1403 distFxn = 'LogNormalDist' 1417 1404 cumeFxn = 'LogNormalCume' 1418 pos = distDict['MinSize']1419 args = [distDict['Mean'],distDict['StdDev']]1420 step = 4.*np.sqrt(np.exp(distDict['StdDev']**2)*(np.exp(distDict['StdDev']**2)-1.))1421 mode = distDict['MinSize']+distDict['Mean']/np.exp(distDict['StdDev']**2)1422 minX = 1. #pos1423 maxX = 1.e4 #np.min([mode+nPoints*step*40,1.e5])1424 1405 elif 'Gauss' in DistName: 1425 1406 distFxn = 'GaussDist' 1426 1407 cumeFxn = 'GaussCume' 1427 pos = distDict['Mean']1428 args = [distDict['StdDev']]1429 step = 0.02*distDict['StdDev']1430 mode = distDict['Mean']1431 minX = np.max([mode-4.*distDict['StdDev'],1.])1432 maxX = np.min([mode+4.*distDict['StdDev'],1.e5])1433 1408 elif 'LSW' in DistName: 1434 1409 distFxn = 'LSWDist' 1435 1410 cumeFxn = 'LSWCume' 1436 pos = distDict['Mean']1437 args = []1438 minX,maxX = [0.,2.*pos]1439 1411 elif 'Schulz' in DistName: 1440 1412 distFxn = 'SchulzZimmDist' 1441 1413 cumeFxn = 'SchulzZimmCume' 1442 pos = distDict['Mean']1443 args = [distDict['StdDev']]1444 minX = np.max([1.,pos-4.*distDict['StdDev']])1445 maxX = np.min([pos+4.*distDict['StdDev'],1.e5])1446 1414 nP = 500 1447 1415 Diam = np.logspace(0.,5.,nP,True) … … 1483 1451 raise Exception("file not found: " + filename) 1484 1452 buf = [line.split() for line in open(filename, 'r').readlines()] 1485 M = len(buf)1486 1453 buf = zip(*buf) # transpose rows and columns 1487 1454 q = np.array(buf[0], dtype=np.float64) -
trunk/GSASIIspc.py
r2544 r2546 515 515 UsymOp = [] 516 516 OprFlg = [] 517 if Nsyms in [1,3]: NunqOp = 0 #Triclinic acentric OR trigonal 3 518 elif Nsyms == 2: #Centric triclinic or acentric monoclinic 519 NunqOp = 1 517 if Nsyms == 2: #Centric triclinic or acentric monoclinic 520 518 UsymOp.append(OprNames[1]) 521 519 OprFlg.append(SGData['SGGen'][1]) 522 520 elif Nsyms == 4: #Point symmetry 2/m, 222, 22m, or 4 523 521 if '4z' in OprNames[1]: #Point symmetry 4 or -4 524 NunqOp = 1525 522 UsymOp.append(OprNames[1]) 526 523 OprFlg.append(SGData['SGGen'][1]) 527 524 elif not SGData['SGInv']: #Acentric Orthorhombic 528 525 if 'm' in OprNames[1:4]: #22m, 2m2 or m22 529 NunqOp = 2530 526 if '2' in OprNames[1]: #Acentric orthorhombic, 2mm 531 527 UsymOp.append(OprNames[2]) … … 544 540 OprFlg.append(SGData['SGGen'][2]) 545 541 else: #Acentric orthorhombic, 222 546 NunqOp = -3547 542 SGData['SGGen'][1:] = [4,2,1] 548 543 UsymOp.append(OprNames[1]) … … 553 548 OprFlg.append(SGData['SGGen'][3]) 554 549 else: #Centric Monoclinic 555 NunqOp = 2556 550 UsymOp.append(OprNames[1]) 557 551 OprFlg.append(SGData['SGGen'][1]) … … 559 553 OprFlg.append(SGData['SGGen'][3]) 560 554 elif Nsyms == 6: #Point symmetry 32, 3m or 6 561 NunqOp = 1562 555 if '6' in OprNames[1]: #Hexagonal 6/m Laue symmetry 563 556 UsymOp.append(OprNames[1]) … … 569 562 elif Nsyms == 8: #Point symmetry mmm, 4/m, or 422, etc 570 563 if '4' in OprNames[1]: #Tetragonal 571 NunqOp = 2572 564 if SGData['SGInv']: #4/m 573 565 UsymOp.append(OprNames[1]) … … 587 579 OprFlg.append(19) 588 580 else: #Orthorhombic, mmm 589 NunqOp = 3590 581 UsymOp.append(OprNames[1]) 591 582 OprFlg.append(SGData['SGGen'][1]) … … 595 586 OprFlg.append(SGData['SGGen'][7]) 596 587 elif Nsyms == 12 and '3' in OprNames[1] and SGData['SGSys'] != 'cubic': #Trigonal 597 NunqOp = 2598 588 UsymOp.append(OprNames[3]) 599 589 OprFlg.append(SGData['SGGen'][3]) … … 601 591 OprFlg.append(SGData['SGGen'][9]) 602 592 elif Nsyms == 12 and '6' in OprNames[1]: #Hexagonal 603 NunqOp = 2604 593 if 'mz' in OprNames[9]: #6/m 605 594 UsymOp.append(OprNames[1]) … … 615 604 elif Nsyms in [16,24]: 616 605 if '3' in OprNames[1]: 617 NunqOp = 1618 606 UsymOp.append('') 619 607 OprFlg.append(SGData['SGGen'][3]) … … 627 615 OprFlg[-1] = 24 628 616 else: #4/mmm or 6/mmm 629 NunqOp = 3630 617 UsymOp.append(' mz ') 631 618 OprFlg.append(1) … … 642 629 else: #System is cubic 643 630 if Nsyms == 48: 644 NunqOp = 2645 631 UsymOp.append(' mx ') 646 632 OprFlg.append(4) 647 633 UsymOp.append(' m110 ') 648 634 OprFlg.append(24) 649 else:650 NunqOp = 0651 635 ncv = len(SGData['SGCen']) 652 636 if ncv > 1: -
trunk/GSASIIstrIO.py
r2520 r2546 12 12 13 13 ''' 14 import sys15 14 import os 16 15 import os.path as ospath … … 18 17 import math 19 18 import copy 20 import random21 19 import cPickle 22 20 import numpy as np 23 21 import numpy.ma as ma 24 import numpy.linalg as nl25 import scipy.optimize as so26 22 import GSASIIpath 27 23 GSASIIpath.SetVersionNumber("$Revision$") 28 24 import GSASIIElem as G2el 29 import GSASIIgrid as G2gd30 25 import GSASIIlattice as G2lat 31 26 import GSASIIspc as G2spc … … 312 307 ''' 313 308 fl = open(GPXfile,'rb') 314 General = {}315 Atoms = []316 309 while True: 317 310 try: … … 787 780 for i,mag in enumerate(RBmags): 788 781 name = '::RBV;'+str(i)+':'+str(irb) 789 mag = parmDict[name]790 782 if name in sigDict: 791 783 VectSig.append(sigDict[name]) … … 1078 1070 phaseVary = [] 1079 1071 phaseDict = {} 1080 phaseConstr = {}1081 1072 pawleyLookup = {} 1082 1073 FFtables = {} #scattering factors - xrays … … 1084 1075 BLtables = {} # neutrons 1085 1076 Natoms = {} 1086 AtMults = {}1087 AtIA = {}1088 1077 maxSSwave = {} 1089 1078 shModels = ['cylindrical','none','shear - 2/m','rolling - mmm'] … … 1516 1505 print ' HKL grid neg esd sum neg num neg use unit? unit esd ' 1517 1506 for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]: 1518 PH = np.array(hkl)1519 1507 phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData) 1520 1508 ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData) … … 1686 1674 print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:' \ 1687 1675 %(at[ct-1],at[cs],Stype,waveType) 1688 line = ''1689 1676 for iw,wave in enumerate(Waves): 1690 1677 stiw = ':'+str(i)+':'+str(iw) … … 2176 2163 hapVary = [] 2177 2164 controlDict = {} 2178 poType = {}2179 poAxes = {}2180 spAxes = {}2181 spType = {}2182 2165 2183 2166 for phase in Phases: … … 2203 2186 limits = Histogram['Limits'][1] 2204 2187 inst = Histogram['Instrument Parameters'][0] 2205 Zero = inst['Zero'][0]2206 2188 if 'C' in inst['Type'][1]: 2207 2189 try: … … 2226 2208 if hapData['HStrain'][1][i]: 2227 2209 hapVary.append(pfx+name) 2228 DIJS = G2spc.HStrainVals(HSvals,SGData)2229 2210 controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0] 2230 2211 if hapData['Pref.Ori.'][0] == 'MD': … … 3082 3063 Back = Background[0] 3083 3064 DebyePeaks = Background[1] 3084 lenBack = len(Back[3:])3085 3065 valstr = ' value : ' 3086 3066 sigstr = ' sig : ' -
trunk/GSASIIstrMain.py
r2500 r2546 13 13 ########### SVN repository information ################### 14 14 import sys 15 import os16 15 import os.path as ospath 17 16 import time 18 17 import math 19 18 import copy 20 import random21 import cPickle22 19 import numpy as np 23 import numpy.ma as ma24 20 import numpy.linalg as nl 25 21 import scipy.optimize as so … … 108 104 # for item in table: print item,table[item] #useful debug - are things shifting? 109 105 break #refinement succeeded - finish up! 110 except TypeError ,FloatingPointError: #result[1] is None on singular matrix106 except TypeError: #result[1] is None on singular matrix 111 107 IfOK = False 112 108 if not len(varyList): … … 316 312 saveVaryList[i] = item 317 313 SeqResult['varyList'] = saveVaryList 318 origvaryList = varyList[:]319 314 parmDict = {} 320 315 parmDict.update(phaseDict) … … 493 488 varyList = covData['varyList'] 494 489 pfx = str(DisAglData['pId'])+'::' 495 A = G2lat.cell2A(Cell[:6])496 cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)497 names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']498 valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]499 490 500 491 Factor = DisAglCtls['Factors'] … … 579 570 :param file out: file object for output. Defaults to sys.stdout. 580 571 ''' 581 import numpy.ma as ma582 572 def MyPrint(s): 583 573 out.write(s+'\n') … … 608 598 if 'covData' in DisAglData: 609 599 covData = DisAglData['covData'] 610 covMatrix = covData['covMatrix']611 varyList = covData['varyList']612 600 pfx = str(DisAglData['pId'])+'::' 613 601 A = G2lat.cell2A(Cell[:6]) … … 626 614 AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData) 627 615 origAtoms = DisAglData['OrigAtoms'] 628 targAtoms = DisAglData['TargAtoms']629 616 for Oatom in origAtoms: 630 617 i = Oatom[0] … … 673 660 if 'covData' in DATData: 674 661 covData = DATData['covData'] 675 covMatrix = covData['covMatrix']676 varyList = covData['varyList']677 662 pfx = str(DATData['pId'])+'::' 678 663 Datoms = [] … … 752 737 print 'ERROR - ',GPXfile," doesn't exist!" 753 738 exit() 754 GPXpath = ospath.dirname(arg[1])755 739 Refine(GPXfile,None) 756 740 else: -
trunk/GSASIIstrMath.py
r2539 r2546 16 16 import numpy.ma as ma 17 17 import numpy.linalg as nl 18 import scipy.optimize as so19 18 import scipy.stats as st 20 19 import GSASIIpath … … 148 147 atxIds = ['dAx:','dAy:','dAz:'] 149 148 atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:'] 150 TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']151 LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']152 SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']153 PIds = ['Px:','Py:','Pz:']154 149 OIds = ['Oa:','Oi:','Oj:','Ok:'] 155 150 RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds … … 177 172 VModel = RBData['Vector'][RBObj['RBId']] 178 173 Q = RBObj['Orient'][0] 179 Pos = RBObj['Orig'][0]180 174 jrb = VRBIds.index(RBObj['RBId']) 181 175 rbsx = str(irb)+':'+str(jrb) … … 238 232 for irb,RBObj in enumerate(RBModels.get('Residue',[])): 239 233 Q = RBObj['Orient'][0] 240 Pos = RBObj['Orig'][0]241 234 jrb = RRBIds.index(RBObj['RBId']) 242 235 torData = RBData['Residue'][RBObj['RBId']]['rbSeq'] … … 521 514 hkl = np.array(hkl) 522 515 if np.any(lasthkl-hkl): 523 PH = np.array(hkl)524 516 phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData) 525 517 ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData) … … 556 548 hkl = np.array(HKLs[int(pnames[3])]) 557 549 if np.any(lasthkl-hkl): 558 PH = np.array(hkl)559 550 phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData) 560 551 SH3Coef = {} … … 589 580 IAdata = Natoms*[' ',] 590 581 Fdata = np.zeros(Natoms) 591 FFdata = []592 BLdata = []593 582 Xdata = np.zeros((3,Natoms)) 594 583 dXdata = np.zeros((3,Natoms)) … … 718 707 #reflection processing begins here - big arrays! 719 708 iBeg = 0 720 time0 = time.time()721 709 while iBeg < nRef: 722 710 iFin = min(iBeg+blkSize,nRef) … … 833 821 SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 834 822 SGT = np.array([ops[1] for ops in SGData['SGOps']]) 835 Ncen = len(SGData['SGCen'])836 Nops = len(SGMT)837 823 FFtables = calcControls['FFtables'] 838 824 BLtables = calcControls['BLtables'] 839 MFtables = calcControls['MFtables']840 825 Amat,Bmat = G2lat.Gmat2AB(G) 841 826 nRef = len(refDict['RefList']) … … 862 847 Flack = 1.-2.*parmDict[phfx+'Flack'] 863 848 time0 = time.time() 864 nref = len(refDict['RefList'])/100865 849 #reflection processing begins here - big arrays! 866 850 iBeg = 0 … … 896 880 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT) 897 881 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 898 Hij = np.reshape(np.array([G2lat.UijtoU6( Uij) for Uij in Hij]),(-1,len(SGT),6))882 Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6)) 899 883 fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr 900 884 if len(FPP.shape) > 1: … … 989 973 :returns: dict dFdvDict: dictionary of derivatives 990 974 ''' 991 phfx = pfx.split(':')[0]+hfx992 975 ast = np.sqrt(np.diag(G)) 993 976 Mast = twopisq*np.multiply.outer(ast,ast) … … 996 979 Ncen = len(SGData['SGCen']) 997 980 Nops = len(SGMT)*Ncen*(1+SGData['SGInv']) 998 MFtables = calcControls['MFtables']999 981 Amat,Bmat = G2lat.Gmat2AB(G) 1000 982 nRef = len(refDict['RefList']) … … 1029 1011 dFdua = np.zeros((nRef,mSize,6)) 1030 1012 time0 = time.time() 1031 nref = len(refDict['RefList'])/1001032 1013 #reflection processing begins here - big arrays! 1033 1014 iBeg = 0 … … 1048 1029 Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T 1049 1030 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 1050 Hij = np.reshape(np.array([G2lat.UijtoU6( Uij) for Uij in Hij]),(-1,len(SGT),6))1031 Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6)) 1051 1032 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1052 1033 MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T #Nref,Natm … … 1212 1193 HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1) 1213 1194 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 1214 Hij = np.reshape(np.array([G2lat.UijtoU6( Uij) for Uij in Hij]),(-1,nTwin,len(SGT),6))1195 Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,nTwin,len(SGT),6)) 1215 1196 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 1216 1197 Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*Mdata*Fdata/len(SGMT) … … 1239 1220 dfbdfr = np.sum(np.sum(fb/occ,axis=-2),axis=0) #Fdata != 0 avoids /0. problem 1240 1221 dfadba /= 2. 1241 dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.1222 # dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2. 1242 1223 dfbdui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),axis=0) 1243 1224 dfbdx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fbx[:,:,:,:,:,nxs],axis=-3),axis=0) … … 1248 1229 dfbdui = np.zeros_like(dfadui) 1249 1230 dfbdua = np.zeros_like(dfadua) 1250 dfbdba = np.zeros_like(dfadba)1231 # dfbdba = np.zeros_like(dfadba) 1251 1232 SA = fas[0]+fas[1] 1252 1233 SB = fbs[0]+fbs[1] … … 1306 1287 SGInv = SGData['SGInv'] 1307 1288 SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 1308 SGT = np.array([ops[1] for ops in SGData['SGOps']])1309 1289 SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']]) 1310 1290 SSGT = np.array([ops[1] for ops in SSGData['SSGOps']]) 1311 1291 FFtables = calcControls['FFtables'] 1312 1292 BLtables = calcControls['BLtables'] 1293 MFtables = calcControls['MFtables'] 1313 1294 Flack = 1.0 1314 1295 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: … … 1425 1406 SGInv = SGData['SGInv'] 1426 1407 SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 1427 SGT = np.array([ops[1] for ops in SGData['SGOps']])1428 1408 SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']]) 1429 1409 SSGT = np.array([ops[1] for ops in SSGData['SSGOps']]) 1430 1410 FFtables = calcControls['FFtables'] 1431 1411 BLtables = calcControls['BLtables'] 1412 MFtables = calcControls['MFtables'] 1432 1413 Flack = 1.0 1433 1414 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: … … 1565 1546 SGInv = SGData['SGInv'] 1566 1547 SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 1567 SGT = np.array([ops[1] for ops in SGData['SGOps']])1568 1548 SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']]) 1569 1549 SSGT = np.array([ops[1] for ops in SSGData['SSGOps']]) … … 1600 1580 dFdbab = np.zeros((nRef,2)) 1601 1581 dFdfl = np.zeros((nRef)) 1602 dFdtw = np.zeros((nRef))1603 1582 dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2)) 1604 1583 dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6)) … … 1636 1615 HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1) #ops x atoms 1637 1616 Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3 1638 Hij = np.array([G2lat.UijtoU6( Uij) for Uij in Hij]) #atoms x 61617 Hij = np.array([G2lat.UijtoU6(uij) for uij in Hij]) #atoms x 6 1639 1618 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) #ops x atoms 1640 1619 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0] #ops x atoms … … 1813 1792 dFdbab = np.zeros((nRef,2)) 1814 1793 dFdfl = np.zeros((nRef)) 1815 dFdtw = np.zeros((nRef))1816 1794 dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2)) 1817 1795 dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6)) … … 1838 1816 FP = np.repeat(FP.T,len(SGT),axis=0) 1839 1817 FPP = np.repeat(FPP.T,len(SGT),axis=0) 1840 dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)1818 # dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor) 1841 1819 Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)) 1842 1820 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) … … 1859 1837 HbH = -np.sum(UniqP[:,:,nxs,:3]*np.inner(UniqP[:,:,:3],bij),axis=-1) #ops x atoms 1860 1838 Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in np.reshape(UniqP,(-1,3))]) #atoms x 3x3 1861 Hij = np.reshape(np.array([G2lat.UijtoU6( Uij) for Uij in Hij]),(iFin-iBeg,-1,6)) #atoms x 61839 Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(iFin-iBeg,-1,6)) #atoms x 6 1862 1840 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) #ops x atoms 1863 1841 # GSASIIpath.IPyBreak() … … 1883 1861 dfadfr = np.sum(fag/occ,axis=-2) #Fdata != 0 ever avoids /0. problem 1884 1862 dfbdfr = np.sum(fbg/occ,axis=-2) #Fdata != 0 avoids /0. problem 1885 dfadba = np.sum(-cosp*Tcorr,axis=-2)1886 dfbdba = np.sum(-sinp*Tcorr,axis=-2)1863 # dfadba = np.sum(-cosp*Tcorr,axis=-2) 1864 # dfbdba = np.sum(-sinp*Tcorr,axis=-2) 1887 1865 dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fag,axis=-2) 1888 1866 dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fbg,axis=-2) … … 2000 1978 SGInv = SGData['SGInv'] 2001 1979 SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 2002 SGT = np.array([ops[1] for ops in SGData['SGOps']])2003 1980 SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']]) 2004 1981 SSGT = np.array([ops[1] for ops in SSGData['SSGOps']]) … … 2011 1988 NM = calcControls[phfx+'TwinNMN']+1 2012 1989 TwinLaw = calcControls[phfx+'TwinLaw'] 2013 TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])2014 1990 TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1)) 2015 1991 nTwin = len(TwinLaw) … … 2043 2019 dFdua = np.zeros((nRef,nTwin,mSize,6)) 2044 2020 dFdbab = np.zeros((nRef,nTwin,2)) 2045 dFdfl = np.zeros((nRef,nTwin))2046 2021 dFdtw = np.zeros((nRef,nTwin)) 2047 2022 dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1])) … … 2088 2063 HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1) #ops x atoms 2089 2064 Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3 2090 Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6( Uij) for Uij in Hij]),(nTwin,-1,6)))2065 Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(nTwin,-1,6))) 2091 2066 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) #ops x atoms 2092 2067 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0] #ops x atoms … … 2129 2104 dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1) 2130 2105 # GSASIIpath.IPyBreak() 2131 if not SGData['SGInv'] and len(TwinLaw) == 1: #Flack derivative2132 dfadfl = np.sum(-FPP*Tcorr*sinp)2133 dfbdfl = np.sum(FPP*Tcorr*cosp)2134 else:2135 dfadfl = 1.02136 dfbdfl = 1.02137 2106 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3! 2138 2107 SA = fas[0]+fas[1] #float = A+A' (might be array[nTwin]) … … 2780 2749 sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2 2781 2750 else: #generalized - P.W. Stephens model OK 2782 pwrs = calcControls[phfx+'MuPwrs']2783 2751 Strms = G2spc.MustrainCoeff(refl[:3],SGData) 2784 2752 const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3 … … 3102 3070 if Phase['General'].get('Modulated',False): 3103 3071 SSGData = Phase['General']['SSGData'] 3104 SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])3105 3072 im = 1 #offset in SS reflection list 3106 3073 #?? … … 3113 3080 Vst = np.sqrt(nl.det(G)) #V* 3114 3081 if not Phase['General'].get('doPawley'): 3115 time0 = time.time()3116 3082 if im: 3117 3083 SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) … … 3119 3085 StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict) 3120 3086 # print 'sf calc time: %.3fs'%(time.time()-time0) 3121 time0 = time.time()3122 3087 badPeak = False 3123 3088 for iref,refl in enumerate(refDict['RefList']): … … 3234 3199 depDerivDict[j] = np.zeros(shape=(len(x))) 3235 3200 #print 'dependent vars',dependentVars 3236 lenX = len(x)3237 3201 hId = Histogram['hId'] 3238 3202 hfx = ':%d:'%(hId) … … 3281 3245 if Phase['General'].get('Modulated',False): 3282 3246 SSGData = Phase['General']['SSGData'] 3283 SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])3284 3247 im = 1 #offset in SS reflection list 3285 3248 #?? … … 3292 3255 GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies 3293 3256 if not Phase['General'].get('doPawley'): 3294 time0 = time.time()3295 3257 if im: 3296 3258 dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) … … 3302 3264 # print 'sf-derv time %.3fs'%(time.time()-time0) 3303 3265 ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase) 3304 time0 = time.time()3305 3266 for iref,refl in enumerate(refDict['RefList']): 3306 3267 if im: … … 3537 3498 :returns: 3538 3499 ''' 3539 nobs = Histogram['Residuals']['Nobs']3540 3500 hId = Histogram['hId'] 3541 3501 hfx = ':%d:'%(hId) … … 3546 3506 if Phase['General'].get('Modulated',False): 3547 3507 SSGData = Phase['General']['SSGData'] 3548 SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])3549 3508 im = 1 #offset in SS reflection list 3550 3509 A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] … … 3647 3606 G2mv.Dict2Map(parmDict,varylist) 3648 3607 Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases 3649 nvar = len(varylist)3650 3608 dMdv = np.empty(0) 3651 3609 histoList = Histograms.keys() … … 3702 3660 #fixup H atom positions here? 3703 3661 ApplyRBModels(parmDict,Phases,rigidbodyDict) #,Update=True?? 3704 nvar = len(varylist)3705 3662 Hess = np.empty(0) 3663 Vec = np.empty(0) 3706 3664 histoList = Histograms.keys() 3707 3665 histoList.sort() … … 3840 3798 if Phase['General'].get('Modulated',False): 3841 3799 SSGData = Phase['General']['SSGData'] 3842 SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])3843 3800 im = 1 #offset in SS reflection list 3844 3801 A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] 3845 3802 G,g = G2lat.A2Gmat(A) #recip & real metric tensors 3846 3803 refDict = Histogram['Data'] 3847 time0 = time.time()3848 3804 if im: 3849 3805 if len(TwinLaw) > 1: -
trunk/imports/G2img_1TIF.py
r2539 r2546 134 134 elif Type == 3: 135 135 Value = st.unpack(byteOrd+nVal*'h',File.read(nVal*2)) 136 st.unpack(byteOrd+nVal*'h',File.read(nVal*2)) 136 137 elif Type == 4: 137 138 if Tag in [273,279]:
Note: See TracChangeset
for help on using the changeset viewer.