Changeset 342

Timestamp:

Aug 5, 2011 2:35:43 PM (12 years ago)

Author:

vondreele

Message:

major modifications to get 1st "working" version of Refine
Add GSASIImapvars.py
fix cell refinement in indexing window
single cycle for peak refinement

Location:

trunk

Files:

• GSASII.py (modified) (6 diffs)
• GSASIIIO.py (modified) (4 diffs)
• GSASIIgrid.py (modified) (9 diffs)
• GSASIIimgGUI.py (modified) (1 diff)
• GSASIIindex.py (modified) (4 diffs)
• GSASIIlattice.py (modified) (5 diffs)
• GSASIImapvars.py (added)
• GSASIIphsGUI.py (modified) (25 diffs)
• GSASIIplot.py (modified) (7 diffs)
• GSASIIpwd.py (modified) (13 diffs)
• GSASIIpwdGUI.py (modified) (17 diffs)
• GSASIIspc.py (modified) (5 diffs)
• GSASIIstruct.py (modified) (10 diffs)

trunk/GSASII.py

@@ -147 +147 @@
         self.Refine = parent.Append(help='', id=wxID_REFINE, kind=wx.ITEM_NORMAL,
             text='Refine')
-        self.Refine.Enable(False)
+        self.Refine.Enable(True)
         self.Bind(wx.EVT_MENU, self.OnRefine, id=wxID_REFINE)
         self.Solve = parent.Append(help='', id=wxID_SOLVE, kind=wx.ITEM_NORMAL,
@@ -391 +391 @@
                             self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Peak List'),[])
                             self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
-                            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])             
+                            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])
+                            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Reflection Lists'),{})             
                             self.PatternId = G2gd.GetPatternTreeItemId(self,Id,'Limits')
                     finally:
@@ -746 +747 @@
                     Id = self.PatternTree.AppendItem(parent=self.root,text=outname)
                     if Id:
-                        Sample = {'Scale':[1.0,True],'Type':'Debye-Scherrer','Absorption':[0.0,False],'DisplaceX':[0.0,False],
-                            'DisplaceY':[0.0,False],'Diffuse':[],'Temperature':300.,'Pressure':1.0,'Humidity':0.0,'Voltage':0.0,'Force':0.0}
+                        Sample = {'Scale':[1.0,True],'Type':'Debye-Scherrer','Absorption':[0.0,False],
+                            'DisplaceX':[0.0,False],'DisplaceY':[0.0,False],'Diffuse':[],
+                            'Temperature':300.,'Pressure':1.0,'Humidity':0.0,
+                            'Voltage':0.0,'Force':0.0,'Gonio. radius':200.0}
                         self.PatternTree.SetItemPyData(Id,[[''],[Xsum,Ysum,Wsum,YCsum,YBsum,YDsum]])
                         self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Comments'),Comments)                    
@@ -757 +760 @@
                         self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
                         self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])             
+                        self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Reflection Lists'),{})             
                         self.PatternTree.SelectItem(Id)
                         self.PatternTree.Expand(Id)
@@ -878 +882 @@
         sub = self.PatternTree.AppendItem(parent=sub,text=PhaseName)
         E,SGData = G2spc.SpcGroup('P 1')
-        self.PatternTree.SetItemPyData(sub, \
-            {'General':{'Name':PhaseName,'Type':'nuclear','SGData':SGData,
-            'Cell':[False,10.,10.,10.,90.,90.,90,1000.],
-            'Pawley dmin':1.0},'Atoms':[],'Drawing':{},'Histograms':{},'Pawley ref':[],
-            'Models':{},'SH Texture':{'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
-            'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
-            'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}})
+        self.PatternTree.SetItemPyData(sub, {
+            'General':{
+                'Name':PhaseName,
+                'Type':'nuclear',
+                'SGData':SGData,
+                'Cell':[False,10.,10.,10.,90.,90.,90,1000.],
+                'Pawley dmin':1.0},
+            'Atoms':[],
+            'Drawing':{},
+            'Histograms':{},
+            'Pawley ref':[],
+            'Models':{},
+            'SH Texture':{
+                'Order':0,
+                'Model':'cylindrical',
+                'Sample omega':[False,0.0],
+                'Sample chi':[False,0.0],
+                'Sample phi':[False,0.0],
+                'SH Coeff':[False,{}],
+                'SHShow':False,
+                'PFhkl':[0,0,1],
+                'PFxyz':[0,0,1],
+                'PlotType':'Pole figure'}
+            })
         
     def OnDeletePhase(self,event):
+        #Hmm, also need to delete this phase from Reflection Lists for each PWDR histogram
         if self.dataFrame:
             self.dataFrame.Clear() 
@@ -1323 +1345 @@
        
     def OnRefine(self,event):
+        self.OnFileSave(event)
         #works - but it'd be better if it could restore plots
         G2str.Refine(self.GSASprojectfile)

trunk/GSASIIIO.py

@@ -218 +218 @@
         sig = float(item[6])/rt2ln2x2
         Iobs = float(item[7])*mult
-        PawleyPeaks.append([h,k,l,mult,tth,sig,False,Iobs,0.0,[]])
+        PawleyPeaks.append([h,k,l,mult,tth,False,Iobs,0.0])
         S = File.readline()
         item = S.split()
@@ -801 +801 @@
             for datus in data[1:]:
                 print '    load: ',datus[0]
-                                
                 sub = self.PatternTree.AppendItem(Id,datus[0])
                 self.PatternTree.SetItemPyData(sub,datus[1])
-                                                
             if 'IMG' in datum[0]:                   #retreive image default flag & data if set
                 Data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Image Controls'))
@@ -1017 +1015 @@
     textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
         'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
-        'SHShow':False,'PFhkl':[0,0,1]}
+        'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}
     shNcof = 0
     file = open(filename, 'Ur')
-    Phase = {}
     S = 1
     Expr = [{},{},{},{},{},{},{},{},{}]
@@ -1121 +1118 @@
                 shCoef[key] = float(val)
         textureData['SH Coeff'] = [False,shCoef]
+        
+    Phase = {
+            'General':{
+                'Name':PhaseName,
+                'Type':Ptype,
+                'SGData':SGData,
+                'Cell':[False,]+abc+angles+[Volume,],
+                'Pawley dmin':1.0},
+            'Atoms':Atoms,
+            'Drawing':{},
+            'Histograms':{},
+            'Pawley ref':[],
+            'Models':{},
+            'SH Texture':textureData
+            }
             
-    Phase['General'] = {'Name':PhaseName,'Type':Ptype,'SGData':SGData,
-        'Cell':[False,]+abc+angles+[Volume,],'Pawley ref':[],'Models':{},'SH Texture':textureData}
-    Phase['Atoms'] = Atoms
-    Phase['Drawing'] = {}
-    Phase['Histograms'] = {}
     return Phase
        

trunk/GSASIIgrid.py

@@ -46 +46 @@
 ] = [wx.NewId() for _init_coll_IndPeaks_Items in range(1)]
 
-[ wxID_UNDO,wxID_LSQPEAKFIT,wxID_BFGSPEAKFIT,wxID_RESETSIGGAM,
-] = [wx.NewId() for _init_coll_PEAK_Items in range(4)]
+[ wxID_UNDO,wxID_LSQPEAKFIT,wxID_LSQONECYCLE,wxID_BFGSPEAKFIT,wxID_RESETSIGGAM,
+] = [wx.NewId() for _init_coll_PEAK_Items in range(5)]
 
 [  wxID_INDEXPEAKS, wxID_REFINECELL, wxID_COPYCELL, wxID_MAKENEWPHASE,
 ] = [wx.NewId() for _init_coll_INDEX_Items in range(4)]
+
+[ wxID_SELECTPHASE,
+] = [wx.NewId() for _init_coll_Refl_Items in range(1)]
 
 [ wxID_PDFCOPYCONTROLS, wxID_PDFSAVECONTROLS, wxID_PDFLOADCONTROLS, 
@@ -92 +95 @@
         parent.Append(menu=self.IndexEdit, title='Cell Index/Refine')
         
+    def _init_coll_ReflMenu(self,parent):
+        parent.Append(menu=self.ReflEdit, title='Reflection List')
+
     def _init_coll_PDFMenu(self,parent):
         parent.Append(menu=self.PDFEdit, title='PDF Controls')
@@ -148 +154 @@
         parent.Append(id=wxID_PAWLEYLOAD, kind=wx.ITEM_NORMAL,text='Pawley create',
             help='Initialize Pawley reflection list')
-        parent.Append(id=wxID_PAWLEYIMPORT, kind=wx.ITEM_NORMAL,text='Pawley import',
-            help='Import Pawley reflection list')
+#        parent.Append(id=wxID_PAWLEYIMPORT, kind=wx.ITEM_NORMAL,text='Pawley import',
+#            help='Import Pawley reflection list')
         parent.Append(id=wxID_PAWLEYDELETE, kind=wx.ITEM_NORMAL,text='Pawley delete',
             help='Delete Pawley reflection list')
@@ -156 +162 @@
         parent.Append(help='Load/Reload index peaks from peak list',id=wxID_INDXRELOAD, 
             kind=wx.ITEM_NORMAL,text='Load/Reload')
+            
+    def _init_coll_Refl_Items(self,parent):
+        self.SelectPhase = parent.Append(help='Select phase for reflection list',id=wxID_SELECTPHASE, 
+            kind=wx.ITEM_NORMAL,text='Select phase')
             
     def _init_coll_Image_Items(self,parent):
@@ -191 +201 @@
         self.PeakFit = parent.Append(id=wxID_LSQPEAKFIT, kind=wx.ITEM_NORMAL,text='LSQ PeakFit', 
             help='Peak fitting via least-squares' )
+        self.PFOneCycle = parent.Append(id=wxID_LSQONECYCLE, kind=wx.ITEM_NORMAL,text='LSQ one cycle', 
+            help='One cycle of Peak fitting via least-squares' )
 #        self.PeakFit = parent.Append(id=wxID_BFGSPEAKFIT, kind=wx.ITEM_NORMAL,text='BFGS PeakFit', 
 #            help='Peak fitting via BFGS algorithm' )
@@ -236 +248 @@
         self.IndPeaksMenu = wx.MenuBar()
         self.IndexMenu = wx.MenuBar()
+        self.ReflMenu = wx.MenuBar()
         self.PDFMenu = wx.MenuBar()
         self.AtomEdit = wx.Menu(title='')
@@ -247 +260 @@
         self.IndPeaksEdit = wx.Menu(title='')
         self.IndexEdit = wx.Menu(title='')
+        self.ReflEdit = wx.Menu(title='')
         self.PDFEdit = wx.Menu(title='')
         self._init_coll_AtomsMenu(self.AtomsMenu)
@@ -268 +282 @@
         self._init_coll_IndexMenu(self.IndexMenu)
         self._init_coll_Index_Items(self.IndexEdit)
+        self._init_coll_ReflMenu(self.ReflMenu)
+        self._init_coll_Refl_Items(self.ReflEdit)
         self._init_coll_PDFMenu(self.PDFMenu)
         self._init_coll_PDF_Items(self.PDFEdit)
         self.UnDo.Enable(False)
         self.PeakFit.Enable(False)
+        self.PFOneCycle.Enable(False)
         self.IndexPeaks.Enable(False)
         self.CopyCell.Enable(False)
@@ -879 +896 @@
         else:
             G2plt.PlotPatterns(self)
+    elif self.PatternTree.GetItemText(item) == 'Reflection Lists':
+        self.PatternId = self.PatternTree.GetItemParent(item)
+        self.PickId = item
+        data = self.PatternTree.GetItemPyData(item)
+        self.RefList = ''
+        if len(data):
+            self.RefList = data.keys()[0]
+        G2pdG.UpdateReflectionGrid(self,data)
+        G2plt.PlotPatterns(self)
      

trunk/GSASIIimgGUI.py

@@ -8 +8 @@
 ########### SVN repository information ###################
 import wx
-import wx.grid as wg
 import matplotlib as mpl
 import math

trunk/GSASIIindex.py

@@ -268 +268 @@
     
     values = A2values(ibrav,A)
-    result = so.leastsq(errFit,values,args=(ibrav,Peaks[7],Peaks[4:7],Pwr),
-        full_output=True,factor=0.1)
+    result = so.leastsq(errFit,values,full_output=True,ftol=0.0001,
+        args=(ibrav,Peaks[7],Peaks[4:7],Pwr))
     A = Values2A(ibrav,result[0])
     return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Pwr)**2),A,result
@@ -286 +286 @@
     values = A2values(ibrav,A)
     values.append(Z)
-    result = so.leastsq(errFit,values,args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Pwr),
-        full_output=True)
+    result = so.leastsq(errFit,values,full_output=True,ftol=0.001,
+        args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Pwr))
     A = Values2A(ibrav,result[0][:-1])
     Z = result[0][-1]
@@ -348 +348 @@
     dmin = getDmin(peaks)
     smin = 1.0e10
-    pwr = 4
+    pwr = 8
     maxTries = 10
     OK = False
@@ -380 +380 @@
         Peaks = np.array(peaks).T
         H = Peaks[4:7]
-        Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A))
+        try:
+            Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A))
+        except FloatingPointError:
+            print G2lat.calc_rDsq(H,A)
+            Peaks[8] = 1.0
         peaks = Peaks.T
         

trunk/GSASIIlattice.py

@@ -162 +162 @@
     B = nl.inv(A)
     return A,B
+    
+def U6toUij(U6):
+    '''Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23]
+    returns 
+    '''
+    U = np.zeros(shape=(3,3))
+    U = [
+        [U6[0],  U6[3],  U6[4]], 
+        [U6[3],  U6[1],  U6[5]], 
+        [U6[4],  U6[5],  U6[2]]]
+    return U
         
 def Uij2betaij(Uij,G):
@@ -173 +184 @@
     '''
     pass
+    
+def CosSinAngle(U,V,G):
+    ''' calculate sin & cos of angle betwee U & V in generalized coordinates 
+    defined by metric tensor G
+    input:
+        U & V - 3-vectors assume numpy arrays
+        G - metric tensor for U & V defined space assume numpy array
+    return:
+        cos(phi) & sin(phi)
+    '''
+    u = U/nl.norm(U)
+    v = V/nl.norm(V)
+    cosP = np.inner(u,np.inner(G,v))
+    sinP = np.sqrt(1.0-cosP**2)
+    return cosP,sinP
     
 def criticalEllipse(prob):
@@ -485 +511 @@
     return sortHKLd(HKL,True,False)
     
-def GenHLaue(dmin,SGLaue,SGLatt,SGUniq,A):
+def GenHLaue(dmin,SGData,A):
     '''Generate the crystallographically unique powder diffraction reflections
     for a lattice and Bravais type
-    '''
-# dmin - minimum d-spacing
-# SGLaue - Laue group symbol = '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm',
-#                            '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
-# SGLatt - lattice centering = 'P','A','B','C','I','F'
-# SGUniq - code for unique monoclinic axis = 'a','b','c'
-# A - 6 terms as defined in calc_rDsq
-# returns - HKL = list of [h,k,l,d] sorted with largest d first and is unique 
-# part of reciprocal space ignoring anomalous dispersion
+    Input:
+        dmin - minimum d-spacing
+        SGData - space group dictionary with at least:
+            SGLaue - Laue group symbol = '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm',
+                     '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
+            SGLatt - lattice centering = 'P','A','B','C','I','F'
+            SGUniq - code for unique monoclinic axis = 'a','b','c'
+        A - 6 terms as defined in calc_rDsq
+    Return;
+        HKL = list of [h,k,l,d] sorted with largest d first and is unique 
+            part of reciprocal space ignoring anomalous dispersion
+    '''
     import math
+    SGLaue = SGData['SGLaue']
+    SGLatt = SGData['SGLatt']
+    SGUniq = SGData['SGUniq']
     #finds maximum allowed hkl for given A within dmin
     if SGLaue in ['3R','3mR']:        #Rhombohedral axes
@@ -645 +677 @@
                     if OdfChk(SGLaue,iord,n):
                         coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
+            else:
+                for n in [i-iord for i in range(2*iord+1)]:
+                    if OdfChk(SGLaue,iord,n):
+                        coeffNames.append('C(%d,%d)'%(iord,n))
     return coeffNames
 
@@ -1130 +1166 @@
         system = spdict['SGSys']
 
-        g2list = GenHLaue(dmin,Laue,center,Axis,cell2A(cell))
+        g2list = GenHLaue(dmin,spdict,cell2A(cell))
         #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
         #    print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])

trunk/GSASIIphsGUI.py

@@ -197 +197 @@
             self.dataFrame.SetLabel('Phase Data for '+generalData['Name'])
             self.PatternTree.SetItemText(Item,generalData['Name'])
+            #Hmm, need to change phase name key in Reflection Lists for each histogram
                         
         def OnPhaseType(event):
@@ -208 +209 @@
                         self.dataDisplay.DeletePage(self.dataDisplay.FindPage('Draw Options'))
                         self.dataDisplay.DeletePage(self.dataDisplay.FindPage('Draw Atoms'))
-                        self.dataDisplay.AdvanceSelection()
-                    if not self.dataDisplay.FindPage('Pawley reflections'):      
-                        self.dataDisplay.AddPage(G2gd.GSGrid(self.dataDisplay),'Pawley reflections')
+                    if not self.dataDisplay.FindPage('Pawley reflections'):
+                        self.PawleyRefl = G2gd.GSGrid(self.dataDisplay)      
+                        self.dataDisplay.AddPage(self.PawleyRefl,'Pawley reflections')
             else:
-                TypeTxt.SetValue(generalData['Type'])
-                
+                TypeTxt.SetValue(generalData['Type'])                
                 
         def OnSpaceGroup(event):
@@ -482 +482 @@
             
         AAchoice = ": ,ALA,ARG,ASN,ASP,CYS,GLN,GLU,GLY,HIS,ILE,LEU,LYS,MET,PHE,PRO,SER,THR,TRP,TYR,VAL,MSE,HOH,UNK"
-        Types = [wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,wg.GRID_VALUE_CHOICE+": ,X,XU,U,F,FX,FXU,FU",
-            wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,4', #x,y,z,frac
+        Types = [wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,wg.GRID_VALUE_CHOICE+": ,X,XU,U,F,FX,FXU,FU",]+ \
+            3*[wg.GRID_VALUE_FLOAT+':10,5',]+[wg.GRID_VALUE_FLOAT+':10,4', #x,y,z,frac
             wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,wg.GRID_VALUE_CHOICE+":I,A",]
         Types += 7*[wg.GRID_VALUE_FLOAT+':10,4',]
@@ -490 +490 @@
             colLabels += ['Mx','My','Mz']
             Types[2] = wg.GRID_VALUE_CHOICE+": ,X,XU,U,M,MX,MXU,MU,F,FX,FXU,FU,FM,FMX,FMU,"
-            Types += [
-                wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_FLOAT+':10,4']
+            Types += 3*[wg.GRID_VALUE_FLOAT+':10,4',]
         elif generalData['Type'] == 'macromolecular':
             colLabels = ['res no','residue','chain'] + colLabels
@@ -622 +621 @@
         def ChangeAtomCell(event):
             
-            def chkUij(Uij,CSI):
+            def chkUij(Uij,CSI): #needs to do something!!!
                 return Uij
 
@@ -1031 +1030 @@
         cx,ct,cs = drawingData['atomPtrs']
         atomData = drawingData['Atoms']
-        Types = [wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,
-            wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,5',    #x,y,z
-            wg.GRID_VALUE_STRING,wg.GRID_VALUE_CHOICE+": ,lines,vdW balls,sticks,balls & sticks,ellipsoids,polyhedra",
+        Types = [wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,5',]+ \
+            [wg.GRID_VALUE_STRING,wg.GRID_VALUE_CHOICE+": ,lines,vdW balls,sticks,balls & sticks,ellipsoids,polyhedra",
             wg.GRID_VALUE_CHOICE+": ,type,name,number",wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,]
         styleChoice = [' ','lines','vdW balls','sticks','balls & sticks','ellipsoids','polyhedra']
@@ -1040 +1038 @@
         if generalData['Type'] == 'macromolecular':
             colLabels = ['Residue','1-letter','Chain'] + colLabels
-            Types = [wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING]+Types
+            Types = 3*[wg.GRID_VALUE_STRING,]+Types
             Types[8] = wg.GRID_VALUE_CHOICE+": ,lines,vdW balls,sticks,balls & sticks,ellipsoids,backbone,ribbons,schematic"
             styleChoice = [' ','lines','vdW balls','sticks','balls & sticks','ellipsoids','backbone','ribbons','schematic']
@@ -1990 +1988 @@
             Obj = event.GetEventObject()
             hist,pid = Indx[Obj.GetId()]
-            UseList[hist]['Size'][2][pid] = Obj.GetValue()
+            if UseList[hist]['Size'][0] == 'ellipsoidal':
+                UseList[hist]['Size'][5][pid] = Obj.GetValue()                
+            else:
+                UseList[hist]['Size'][2][pid] = Obj.GetValue()
             
         def OnSizeVal(event):
             Obj = event.GetEventObject()
             hist,pid = Indx[Obj.GetId()]
-            try:
-                size = float(Obj.GetValue())
-                if pid == 0 and size <= 0:
-                    raise ValueError
-                elif pid == 1 and size <= -UseList[hist]['Size'][1][0]:
-                    raise ValueError
-                UseList[hist]['Size'][1][pid] = size
-            except ValueError:
-                pass
+            if UseList[hist]['Size'][0] == 'ellipsoidal':
+                try:
+                    size = float(Obj.GetValue())
+                    if pid < 3 and size <= 0.01:            #10A lower limit!
+                        raise ValueError                    
+                except ValueError:
+                    pass
+            else:
+                try:
+                    size = float(Obj.GetValue())
+                    if size <= 0.01:            #10A lower limit!
+                        raise ValueError
+                    UseList[hist]['Size'][1][pid] = size
+                except ValueError:
+                    pass
             Obj.SetValue("%.1f"%(UseList[hist]['Size'][1][pid]))          #reset in case of error
             
@@ -2046 +2053 @@
                     UseList[hist]['Mustrain'][4][pid] = strain
                 else:
-                    if pid == 0 and strain < 0:
-                        raise ValueError
-                    elif pid == 1 and strain < -UseList[hist]['Mustrain'][1][0]:
+                    if strain <= 0:
                         raise ValueError
                     UseList[hist]['Mustrain'][1][pid] = strain
@@ -2054 +2059 @@
                 pass
             if UseList[hist]['Mustrain'][0] == 'generalized':
-                Obj.SetValue("%.5f"%(UseList[hist]['Mustrain'][4][pid]))          #reset in case of error
+                Obj.SetValue("%.3f"%(UseList[hist]['Mustrain'][4][pid]))          #reset in case of error
             else:
-                Obj.SetValue("%.5f"%(UseList[hist]['Mustrain'][1][pid]))          #reset in case of error
+                Obj.SetValue("%.3f"%(UseList[hist]['Mustrain'][1][pid]))          #reset in case of error
             G2plt.PlotStrain(self,data)
             
@@ -2073 +2078 @@
             Obj.SetValue('%3d %3d %3d'%(h,k,l)) 
             G2plt.PlotStrain(self,data)
-
-        def OnMDRef(event):
+            
+        def OnPOType(event):
             Obj = event.GetEventObject()
             hist = Indx[Obj.GetId()]
-            UseList[hist]['MDtexture'][1] = Obj.GetValue()
-            
-        def OnMDVal(event):
+            if 'March' in POType.GetValue():
+                UseList[hist]['Pref.Ori.'][0] = 'MD'
+            else:
+                UseList[hist]['Pref.Ori.'][0] = 'SH'
+            UpdateDData()            
+
+        def OnPORef(event):
+            Obj = event.GetEventObject()
+            hist = Indx[Obj.GetId()]
+            UseList[hist]['Pref.Ori.'][2] = Obj.GetValue()
+            
+        def OnPOVal(event):
             Obj = event.GetEventObject()
             hist = Indx[Obj.GetId()]
@@ -2085 +2099 @@
                 mdVal = float(Obj.GetValue())
                 if mdVal > 0:
-                    UseList[hist]['MDtexture'][0] = mdVal
+                    UseList[hist]['Pref.Ori.'][1] = mdVal
             except ValueError:
                 pass
-            Obj.SetValue("%.3f"%(UseList[hist]['MDtexture'][0]))          #reset in case of error
-            
-        def OnMDAxis(event):
+            Obj.SetValue("%.3f"%(UseList[hist]['Pref.Ori.'][1]))          #reset in case of error
+            
+        def OnPOAxis(event):
             Obj = event.GetEventObject()
             hist = Indx[Obj.GetId()]
@@ -2097 +2111 @@
                 hkl = [int(Saxis[i]) for i in range(3)]
             except (ValueError,IndexError):
-                hkl = UseList[hist]['MDtexture'][2]
+                hkl = UseList[hist]['Pref.Ori.'][3]
             if not np.any(np.array(hkl)):
-                hkl = UseList[hist]['MDtexture'][2]
-            UseList[hist]['MDtexture'][2] = hkl
+                hkl = UseList[hist]['Pref.Ori.'][3]
+            UseList[hist]['Pref.Ori.'][3] = hkl
             h,k,l = hkl
             Obj.SetValue('%3d %3d %3d'%(h,k,l)) 
             
+        def OnPOOrder(event):
+            Obj = event.GetEventObject()
+            hist = Indx[Obj.GetId()]
+            Order = int(Obj.GetValue())
+            UseList[hist]['Pref.Ori.'][4] = Order
+            UseList[hist]['Pref.Ori.'][5] = SetPOCoef(Order,hist)
+            UpdateDData()
+
+        def SetPOCoef(Order,hist):
+            cofNames = G2lat.GenSHCoeff(SGData['SGLaue'],None,Order)     #cylindrical sample symmetry
+            newPOCoef = dict(zip(cofNames,np.zeros(len(cofNames))))
+            POCoeff = UseList[hist]['Pref.Ori.'][5]
+            for cofName in POCoeff:
+                if cofName in  cofNames:
+                    newPOCoef[cofName] = POCoeff[cofName]
+            return newPOCoef
+        
         def OnExtRef(event):
             Obj = event.GetEventObject()
@@ -2144 +2175 @@
         shOrder.Bind(wx.EVT_COMBOBOX,OnShOrder)
         shSizer.Add(shOrder,0,wx.ALIGN_CENTER_VERTICAL)
-        shSizer.Add((5,0),0)
-        shRef = wx.CheckBox(dataDisplay,label=' Refine texture?')
-        shRef.SetValue(textureData['SH Coeff'][0])
-        shRef.Bind(wx.EVT_CHECKBOX, OnSHRefine)
-        shSizer.Add(shRef,0,wx.ALIGN_CENTER_VERTICAL)
-        shShow = wx.CheckBox(dataDisplay,label=' Show coeff.?')
-        shShow.SetValue(textureData['SHShow'])
-        shShow.Bind(wx.EVT_CHECKBOX, OnSHShow)
-        shSizer.Add(shShow,0,wx.ALIGN_CENTER_VERTICAL)
-        mainSizer.Add(shSizer,0,0)
-        mainSizer.Add((0,5),0)
-        PTSizer = wx.FlexGridSizer(2,4,5,5)
-        PTSizer.Add(wx.StaticText(dataDisplay,-1,' Texture plot type: '),0,wx.ALIGN_CENTER_VERTICAL)
-        choices = ['Axial pole distribution','Pole figure','Inverse pole figure']            
-        pfType = wx.ComboBox(dataDisplay,-1,value=str(textureData['PlotType']),choices=choices,
-            style=wx.CB_READONLY|wx.CB_DROPDOWN)
-        pfType.Bind(wx.EVT_COMBOBOX,OnPfType)
-        PTSizer.Add(pfType,0,wx.ALIGN_CENTER_VERTICAL)
-        PTSizer.Add(wx.StaticText(dataDisplay,-1,' Projection type: '),0,wx.ALIGN_CENTER_VERTICAL)
-        projSel = wx.ComboBox(dataDisplay,-1,value=self.Projection,choices=['equal area','stereographic'],
-            style=wx.CB_READONLY|wx.CB_DROPDOWN)
-        projSel.Bind(wx.EVT_COMBOBOX,OnProjSel)
-        PTSizer.Add(projSel,0,wx.ALIGN_CENTER_VERTICAL)
-        if textureData['PlotType'] in ['Pole figure','Axial pole distribution']:
-            PTSizer.Add(wx.StaticText(dataDisplay,-1,' Pole figure HKL: '),0,wx.ALIGN_CENTER_VERTICAL)
-            PH = textureData['PFhkl']
-            pfVal = wx.TextCtrl(dataDisplay,-1,'%d,%d,%d'%(PH[0],PH[1],PH[2]),style=wx.TE_PROCESS_ENTER)
-        else:
-            PTSizer.Add(wx.StaticText(dataDisplay,-1,' Inverse pole figure XYZ: '),0,wx.ALIGN_CENTER_VERTICAL)
-            PX = textureData['PFxyz']
-            pfVal = wx.TextCtrl(dataDisplay,-1,'%3.1f,%3.1f,%3.1f'%(PX[0],PX[1],PX[2]),style=wx.TE_PROCESS_ENTER)
-        pfVal.Bind(wx.EVT_TEXT_ENTER,OnPFValue)
-        pfVal.Bind(wx.EVT_KILL_FOCUS,OnPFValue)
-        PTSizer.Add(pfVal,0,wx.ALIGN_CENTER_VERTICAL)
-        PTSizer.Add(wx.StaticText(dataDisplay,-1,' Color scheme'),0,wx.ALIGN_CENTER_VERTICAL)
-        choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
-        choice.sort()
-        colorSel = wx.ComboBox(dataDisplay,-1,value=self.ContourColor,choices=choice,
-            style=wx.CB_READONLY|wx.CB_DROPDOWN)
-        colorSel.Bind(wx.EVT_COMBOBOX,OnColorSel)
-        PTSizer.Add(colorSel,0,wx.ALIGN_CENTER_VERTICAL)        
-        mainSizer.Add(PTSizer,0,wx.ALIGN_CENTER_VERTICAL)
-        
-        mainSizer.Add((0,5),0)
-        if textureData['SHShow']:
-            mainSizer.Add(wx.StaticText(dataDisplay,-1,'Spherical harmonic coefficients: '),0,wx.ALIGN_CENTER_VERTICAL)
+        if textureData['Order']:
+            shSizer.Add((5,0),0)
+            shRef = wx.CheckBox(dataDisplay,label=' Refine texture?')
+            shRef.SetValue(textureData['SH Coeff'][0])
+            shRef.Bind(wx.EVT_CHECKBOX, OnSHRefine)
+            shSizer.Add(shRef,0,wx.ALIGN_CENTER_VERTICAL)
+            shShow = wx.CheckBox(dataDisplay,label=' Show coeff.?')
+            shShow.SetValue(textureData['SHShow'])
+            shShow.Bind(wx.EVT_CHECKBOX, OnSHShow)
+            shSizer.Add(shShow,0,wx.ALIGN_CENTER_VERTICAL)
+            mainSizer.Add(shSizer,0,0)
             mainSizer.Add((0,5),0)
-            ODFSizer = wx.FlexGridSizer(2,8,2,2)
-            ODFIndx = {}
-            ODFkeys = textureData['SH Coeff'][1].keys()
-            ODFkeys.sort()
-            for item in ODFkeys:
-                ODFSizer.Add(wx.StaticText(dataDisplay,-1,item),0,wx.ALIGN_CENTER_VERTICAL)
-                ODFval = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%8.3f'%(textureData['SH Coeff'][1][item]),style=wx.TE_PROCESS_ENTER)
-                ODFIndx[ODFval.GetId()] = item
-                ODFval.Bind(wx.EVT_TEXT_ENTER,OnODFValue)
-                ODFval.Bind(wx.EVT_KILL_FOCUS,OnODFValue)
-                ODFSizer.Add(ODFval,0,wx.ALIGN_CENTER_VERTICAL)
-            mainSizer.Add(ODFSizer,0,wx.ALIGN_CENTER_VERTICAL)
+            PTSizer = wx.FlexGridSizer(2,4,5,5)
+            PTSizer.Add(wx.StaticText(dataDisplay,-1,' Texture plot type: '),0,wx.ALIGN_CENTER_VERTICAL)
+            choices = ['Axial pole distribution','Pole figure','Inverse pole figure']            
+            pfType = wx.ComboBox(dataDisplay,-1,value=str(textureData['PlotType']),choices=choices,
+                style=wx.CB_READONLY|wx.CB_DROPDOWN)
+            pfType.Bind(wx.EVT_COMBOBOX,OnPfType)
+            PTSizer.Add(pfType,0,wx.ALIGN_CENTER_VERTICAL)
+            PTSizer.Add(wx.StaticText(dataDisplay,-1,' Projection type: '),0,wx.ALIGN_CENTER_VERTICAL)
+            projSel = wx.ComboBox(dataDisplay,-1,value=self.Projection,choices=['equal area','stereographic'],
+                style=wx.CB_READONLY|wx.CB_DROPDOWN)
+            projSel.Bind(wx.EVT_COMBOBOX,OnProjSel)
+            PTSizer.Add(projSel,0,wx.ALIGN_CENTER_VERTICAL)
+            if textureData['PlotType'] in ['Pole figure','Axial pole distribution']:
+                PTSizer.Add(wx.StaticText(dataDisplay,-1,' Pole figure HKL: '),0,wx.ALIGN_CENTER_VERTICAL)
+                PH = textureData['PFhkl']
+                pfVal = wx.TextCtrl(dataDisplay,-1,'%d,%d,%d'%(PH[0],PH[1],PH[2]),style=wx.TE_PROCESS_ENTER)
+            else:
+                PTSizer.Add(wx.StaticText(dataDisplay,-1,' Inverse pole figure XYZ: '),0,wx.ALIGN_CENTER_VERTICAL)
+                PX = textureData['PFxyz']
+                pfVal = wx.TextCtrl(dataDisplay,-1,'%3.1f,%3.1f,%3.1f'%(PX[0],PX[1],PX[2]),style=wx.TE_PROCESS_ENTER)
+            pfVal.Bind(wx.EVT_TEXT_ENTER,OnPFValue)
+            pfVal.Bind(wx.EVT_KILL_FOCUS,OnPFValue)
+            PTSizer.Add(pfVal,0,wx.ALIGN_CENTER_VERTICAL)
+            PTSizer.Add(wx.StaticText(dataDisplay,-1,' Color scheme'),0,wx.ALIGN_CENTER_VERTICAL)
+            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
+            choice.sort()
+            colorSel = wx.ComboBox(dataDisplay,-1,value=self.ContourColor,choices=choice,
+                style=wx.CB_READONLY|wx.CB_DROPDOWN)
+            colorSel.Bind(wx.EVT_COMBOBOX,OnColorSel)
+            PTSizer.Add(colorSel,0,wx.ALIGN_CENTER_VERTICAL)        
+            mainSizer.Add(PTSizer,0,wx.ALIGN_CENTER_VERTICAL)
+            
             mainSizer.Add((0,5),0)
-        mainSizer.Add((0,5),0)
-        mainSizer.Add(wx.StaticText(dataDisplay,-1,'Sample orientation angles: '),0,wx.ALIGN_CENTER_VERTICAL)
-        mainSizer.Add((0,5),0)
-        angSizer = wx.BoxSizer(wx.HORIZONTAL)
-        angIndx = {}
-        valIndx = {}
-        for item in ['Sample omega','Sample chi','Sample phi']:
-            angRef = wx.CheckBox(dataDisplay,label=item+': ')
-            angRef.SetValue(textureData[item][0])
-            angIndx[angRef.GetId()] = item
-            angRef.Bind(wx.EVT_CHECKBOX, OnAngRef)
-            angSizer.Add(angRef,0,wx.ALIGN_CENTER_VERTICAL)
-            angVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%8.2f'%(textureData[item][1]),style=wx.TE_PROCESS_ENTER)
-            valIndx[angVal.GetId()] = item
-            angVal.Bind(wx.EVT_TEXT_ENTER,OnAngValue)
-            angVal.Bind(wx.EVT_KILL_FOCUS,OnAngValue)
-            angSizer.Add(angVal,0,wx.ALIGN_CENTER_VERTICAL)
-            angSizer.Add((5,0),0)
-        mainSizer.Add(angSizer,0,wx.ALIGN_CENTER_VERTICAL)
+            if textureData['SHShow']:
+                mainSizer.Add(wx.StaticText(dataDisplay,-1,'Spherical harmonic coefficients: '),0,wx.ALIGN_CENTER_VERTICAL)
+                mainSizer.Add((0,5),0)
+                ODFSizer = wx.FlexGridSizer(2,8,2,2)
+                ODFIndx = {}
+                ODFkeys = textureData['SH Coeff'][1].keys()
+                ODFkeys.sort()
+                for item in ODFkeys:
+                    ODFSizer.Add(wx.StaticText(dataDisplay,-1,item),0,wx.ALIGN_CENTER_VERTICAL)
+                    ODFval = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%8.3f'%(textureData['SH Coeff'][1][item]),style=wx.TE_PROCESS_ENTER)
+                    ODFIndx[ODFval.GetId()] = item
+                    ODFval.Bind(wx.EVT_TEXT_ENTER,OnODFValue)
+                    ODFval.Bind(wx.EVT_KILL_FOCUS,OnODFValue)
+                    ODFSizer.Add(ODFval,0,wx.ALIGN_CENTER_VERTICAL)
+                mainSizer.Add(ODFSizer,0,wx.ALIGN_CENTER_VERTICAL)
+                mainSizer.Add((0,5),0)
+            mainSizer.Add((0,5),0)
+            mainSizer.Add(wx.StaticText(dataDisplay,-1,'Sample orientation angles: '),0,wx.ALIGN_CENTER_VERTICAL)
+            mainSizer.Add((0,5),0)
+            angSizer = wx.BoxSizer(wx.HORIZONTAL)
+            angIndx = {}
+            valIndx = {}
+            for item in ['Sample omega','Sample chi','Sample phi']:
+                angRef = wx.CheckBox(dataDisplay,label=item+': ')
+                angRef.SetValue(textureData[item][0])
+                angIndx[angRef.GetId()] = item
+                angRef.Bind(wx.EVT_CHECKBOX, OnAngRef)
+                angSizer.Add(angRef,0,wx.ALIGN_CENTER_VERTICAL)
+                angVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%8.2f'%(textureData[item][1]),style=wx.TE_PROCESS_ENTER)
+                valIndx[angVal.GetId()] = item
+                angVal.Bind(wx.EVT_TEXT_ENTER,OnAngValue)
+                angVal.Bind(wx.EVT_KILL_FOCUS,OnAngValue)
+                angSizer.Add(angVal,0,wx.ALIGN_CENTER_VERTICAL)
+                angSizer.Add((5,0),0)
+            mainSizer.Add(angSizer,0,wx.ALIGN_CENTER_VERTICAL)
+        else:  #finish the texture output when order = 0
+            mainSizer.Add(shSizer,0,0)
+            mainSizer.Add((0,5),0)            
         mainSizer.Add(wx.StaticText(dataDisplay,-1,'Histogram data for '+PhaseName+':'),0,wx.ALIGN_CENTER_VERTICAL)
         for item in keyList:
@@ -2235 +2270 @@
             if UseList[item]['Show']:
                 scaleSizer = wx.BoxSizer(wx.HORIZONTAL)
-                scaleRef = wx.CheckBox(dataDisplay,label=' Scale factor: ')
+                scaleRef = wx.CheckBox(dataDisplay,label=' Phase fraction: ')
                 scaleRef.SetValue(UseList[item]['Scale'][1])
                 Indx[scaleRef.GetId()] = item
@@ -2252 +2287 @@
                 mainSizer.Add((0,5),0)
                 sizeSizer = wx.BoxSizer(wx.HORIZONTAL)
-                choices = ['isotropic','uniaxial',]
+                sizeSizer.Add(wx.StaticText(dataDisplay,-1,' Size model: '),0,wx.ALIGN_CENTER_VERTICAL)
+                choices = ['isotropic','uniaxial','ellipsoidal']
                 sizeType = wx.ComboBox(dataDisplay,wx.ID_ANY,value=UseList[item]['Size'][0],choices=choices,
                     style=wx.CB_READONLY|wx.CB_DROPDOWN)
@@ -2300 +2336 @@
                     sizeSizer.Add((5,0),0)                    
                     mainSizer.Add(sizeSizer)
+                elif UseList[item]['Size'][0] == 'ellipsoidal':
+                    mainSizer.Add(sizeSizer)
+                    mainSizer.Add((0,5),0)
+                    sizeSizer = wx.BoxSizer(wx.HORIZONTAL)
+                    parms = zip(['S11','S22','S33','S12','S13','S23'],UseList[item]['Size'][4],
+                        UseList[item]['Size'][5],range(6))
+                    for Pa,val,ref,id in parms:
+                        sizeRef = wx.CheckBox(dataDisplay,label=Pa)
+                        sizeRef.SetValue(ref)
+                        Indx[sizeRef.GetId()] = [item,id]
+                        sizeRef.Bind(wx.EVT_CHECKBOX, OnSizeRef)
+                        sizeSizer.Add(sizeRef,0,wx.ALIGN_CENTER_VERTICAL)
+                        sizeVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%.1f'%(val),style=wx.TE_PROCESS_ENTER)
+                        Indx[sizeVal.GetId()] = [item,id]
+                        sizeVal.Bind(wx.EVT_TEXT_ENTER,OnSizeVal)
+                        sizeVal.Bind(wx.EVT_KILL_FOCUS,OnSizeVal)
+                        sizeSizer.Add(sizeVal,0,wx.ALIGN_CENTER_VERTICAL)
+                        sizeSizer.Add((5,0),0)
+                    sizeSizer.Add((5,0),0)                    
+                    mainSizer.Add(sizeSizer)
                 
                 strainSizer = wx.BoxSizer(wx.HORIZONTAL)
+                strainSizer.Add(wx.StaticText(dataDisplay,-1,' Mustrain model: '),0,wx.ALIGN_CENTER_VERTICAL)
                 choices = ['isotropic','uniaxial','generalized',]
                 strainType = wx.ComboBox(dataDisplay,wx.ID_ANY,value=UseList[item]['Mustrain'][0],choices=choices,
@@ -2373 +2430 @@
                         strainSizer.Add(strainVal,0,wx.ALIGN_CENTER_VERTICAL)
                     mainSizer.Add(strainSizer)
-                #MD texture  'MDtexture':[1.0,False,[0,0,1]]
-                mdSizer = wx.BoxSizer(wx.HORIZONTAL)
-                mdRef = wx.CheckBox(dataDisplay,label=' March-Dollase texture ratio: ')
-                mdRef.SetValue(UseList[item]['MDtexture'][1])
-                Indx[mdRef.GetId()] = item
-                mdRef.Bind(wx.EVT_CHECKBOX, OnMDRef)
-                mdSizer.Add(mdRef,0,wx.ALIGN_CENTER_VERTICAL)
-                mdVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,
-                    '%.3f'%(UseList[item]['MDtexture'][0]),style=wx.TE_PROCESS_ENTER)
-                Indx[mdVal.GetId()] = item
-                mdVal.Bind(wx.EVT_TEXT_ENTER,OnMDVal)
-                mdVal.Bind(wx.EVT_KILL_FOCUS,OnMDVal)
-                mdSizer.Add(mdVal,0,wx.ALIGN_CENTER_VERTICAL)
-                mdSizer.Add(wx.StaticText(dataDisplay,-1,' Unique axis, H K L: '),0,wx.ALIGN_CENTER_VERTICAL)
-                h,k,l = UseList[item]['MDtexture'][2]
-                mdAxis = wx.TextCtrl(dataDisplay,-1,'%3d %3d %3d'%(h,k,l),style=wx.TE_PROCESS_ENTER)
-                Indx[mdAxis.GetId()] = item
-                mdAxis.Bind(wx.EVT_TEXT_ENTER,OnMDAxis)
-                mdAxis.Bind(wx.EVT_KILL_FOCUS,OnMDAxis)
-                mdSizer.Add(mdAxis,0,wx.ALIGN_CENTER_VERTICAL)
-                mainSizer.Add(mdSizer)
-                mainSizer.Add((0,5),0)
-                
+                #texture  'Pref. Ori.':['MD',1.0,False,[0,0,1],0,[]] last two for 'SH' are SHorder & coeff
+#                poSizer = wx.BoxSizer(wx.HORIZONTAL)
+                poSizer = wx.FlexGridSizer(1,6,5,5)
+                POData = UseList[item]['Pref.Ori.']
+                choice = ['March-Dollase','Spherical harmonics']
+                POtype = choice[['MD','SH'].index(POData[0])]
+                poSizer.Add(wx.StaticText(dataDisplay,-1,' Preferred orientation model '),0,wx.ALIGN_CENTER_VERTICAL)
+                POType = wx.ComboBox(dataDisplay,wx.ID_ANY,value=POtype,choices=choice,
+                    style=wx.CB_READONLY|wx.CB_DROPDOWN)
+                Indx[POType.GetId()] = item
+                POType.Bind(wx.EVT_COMBOBOX, OnPOType)
+                poSizer.Add(POType)
+                if POData[0] == 'MD':
+                    mainSizer.Add(poSizer)
+                    poSizer = wx.BoxSizer(wx.HORIZONTAL)
+                    poRef = wx.CheckBox(dataDisplay,label=' March-Dollase ratio: ')
+                    poRef.SetValue(POData[2])
+                    Indx[poRef.GetId()] = item
+                    poRef.Bind(wx.EVT_CHECKBOX,OnPORef)
+                    poSizer.Add(poRef,0,wx.ALIGN_CENTER_VERTICAL)
+                    poVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,
+                        '%.3f'%(POData[1]),style=wx.TE_PROCESS_ENTER)
+                    Indx[poVal.GetId()] = item
+                    poVal.Bind(wx.EVT_TEXT_ENTER,OnPOVal)
+                    poVal.Bind(wx.EVT_KILL_FOCUS,OnPOVal)
+                    poSizer.Add(poVal,0,wx.ALIGN_CENTER_VERTICAL)
+                    poSizer.Add(wx.StaticText(dataDisplay,-1,' Unique axis, H K L: '),0,wx.ALIGN_CENTER_VERTICAL)
+                    h,k,l =POData[3]
+                    poAxis = wx.TextCtrl(dataDisplay,-1,'%3d %3d %3d'%(h,k,l),style=wx.TE_PROCESS_ENTER)
+                    Indx[poAxis.GetId()] = item
+                    poAxis.Bind(wx.EVT_TEXT_ENTER,OnPOAxis)
+                    poAxis.Bind(wx.EVT_KILL_FOCUS,OnPOAxis)
+                    poSizer.Add(poAxis,0,wx.ALIGN_CENTER_VERTICAL)
+                    mainSizer.Add(poSizer)
+                else:           #'SH'
+                    poSizer.Add(wx.StaticText(dataDisplay,-1,' Harmonic order: '),0,wx.ALIGN_CENTER_VERTICAL)
+                    poOrder = wx.ComboBox(dataDisplay,wx.ID_ANY,value=str(POData[4]),choices=[str(2*i) for i in range(18)],
+                        style=wx.CB_READONLY|wx.CB_DROPDOWN)
+                    Indx[poOrder.GetId()] = item
+                    poOrder.Bind(wx.EVT_COMBOBOX,OnPOOrder)
+                    poSizer.Add(poOrder,0,wx.ALIGN_CENTER_VERTICAL)
+                    poRef = wx.CheckBox(dataDisplay,label=' Refine? ')
+                    poRef.SetValue(POData[2])
+                    Indx[poRef.GetId()] = item
+                    poRef.Bind(wx.EVT_CHECKBOX,OnPORef)
+                    poSizer.Add(poRef,0,wx.ALIGN_CENTER_VERTICAL)
+                    mainSizer.Add(poSizer)
+                    if POData[4]:
+                        mainSizer.Add(wx.StaticText(dataDisplay,-1,'Spherical harmonic coefficients: '),0,wx.ALIGN_CENTER_VERTICAL)
+                        mainSizer.Add((0,5),0)
+                        ODFSizer = wx.FlexGridSizer(2,8,2,2)
+                        ODFIndx = {}
+                        ODFkeys = POData[5].keys()
+                        ODFkeys.sort()
+                        for odf in ODFkeys:
+                            ODFSizer.Add(wx.StaticText(dataDisplay,-1,odf),0,wx.ALIGN_CENTER_VERTICAL)
+                            ODFval = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%8.3f'%(POData[5][odf]),style=wx.TE_PROCESS_ENTER)
+                            ODFIndx[ODFval.GetId()] = odf
+    #                        ODFval.Bind(wx.EVT_TEXT_ENTER,OnODFValue)
+    #                        ODFval.Bind(wx.EVT_KILL_FOCUS,OnODFValue)
+                            ODFSizer.Add(ODFval,0,wx.ALIGN_CENTER_VERTICAL)
+                        mainSizer.Add(ODFSizer,0,wx.ALIGN_CENTER_VERTICAL)
+                        mainSizer.Add((0,5),0)
+                mainSizer.Add((0,5),0)                
                 #Extinction  'Extinction':[0.0,False]
                 extSizer = wx.BoxSizer(wx.HORIZONTAL)
@@ -2469 +2567 @@
                     for i in result: 
                         histoName = TextList[i]
+                        pId = G2gd.GetPatternTreeItemId(self,self.root,histoName)
                         UseList[histoName] = {'Histogram':histoName,'Show':False,
-                            'Scale':[1.0,False],'MDtexture':[1.0,False,[0,0,1]],
-                            'Size':['isotropic',[10000.,0,],[False,False],[0,0,1]],
-                            'Mustrain':['isotropic',[1.0,0.0],[False,False],[0,0,1],
+                            'Scale':[1.0,False],'Pref.Ori.':['MD',1.0,False,[0,0,1],0,{}],
+                            'Size':['isotropic',[10.,10.,],[False,False],[0,0,1],6*[0.0,],6*[False,]],
+                            'Mustrain':['isotropic',[1.0,1.0],[False,False],[0,0,1],
                                 NShkl*[0.01,],NShkl*[False,]],                            
                             'Extinction':[0.0,False]}
+                        refList = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,pId,'Reflection Lists'))
+                        refList[generalData['Name']] = []                       
                     data['Histograms'] = UseList
                     UpdateDData()
@@ -2501 +2602 @@
 
     def FillPawleyReflectionsGrid():
-        if data['Histograms']:
-            self.dataFrame.PawleyMenu.FindItemById(G2gd.wxID_PAWLEYLOAD).Enable(True)
-            self.dataFrame.PawleyMenu.FindItemById(G2gd.wxID_PAWLEYIMPORT).Enable(True)
-        else:
-            self.dataFrame.PawleyMenu.FindItemById(G2gd.wxID_PAWLEYLOAD).Enable(False)
-            self.dataFrame.PawleyMenu.FindItemById(G2gd.wxID_PAWLEYIMPORT).Enable(False)
                         
         def KeyEditPawleyGrid(event):
-            colList = PawleyRefl.GetSelectedCols()
+            colList = self.PawleyRefl.GetSelectedCols()
             PawleyPeaks = data['Pawley ref']
             if event.GetKeyCode() == wx.WXK_RETURN:
@@ -2518 +2613 @@
                 event.Skip(True)
             elif colList:
-                PawleyRefl.ClearSelection()
+                self.PawleyRefl.ClearSelection()
                 key = event.GetKeyCode()
                 for col in colList:
@@ -2532 +2627 @@
             rowLabels = []
             for i in range(len(PawleyPeaks)): rowLabels.append(str(i+1))
-            colLabels = ['h','k','l','mul','2-theta','sigma','refine','Iobs','Icalc']
-            Types = [wg.GRID_VALUE_LONG,wg.GRID_VALUE_LONG,wg.GRID_VALUE_LONG,wg.GRID_VALUE_LONG,
-                wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
-                wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,5']
+            colLabels = ['h','k','l','mul','d','refine','Fsq(hkl)','sig(Fsq)']
+            Types = 4*[wg.GRID_VALUE_LONG,]+[wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,]+ \
+                2*[wg.GRID_VALUE_FLOAT+':10,2',]
             PawleyTable = G2gd.Table(PawleyPeaks,rowLabels=rowLabels,colLabels=colLabels,types=Types)
-            PawleyRefl.SetTable(PawleyTable, True)
-            PawleyRefl.Bind(wx.EVT_KEY_DOWN, KeyEditPawleyGrid)                 
-            PawleyRefl.SetMargins(0,0)
-            PawleyRefl.AutoSizeColumns(False)
+            self.PawleyRefl.SetTable(PawleyTable, True)
+            self.PawleyRefl.Bind(wx.EVT_KEY_DOWN, KeyEditPawleyGrid)                 
+            self.PawleyRefl.SetMargins(0,0)
+            self.PawleyRefl.AutoSizeColumns(False)
             self.dataFrame.setSizePosLeft([500,300])
                     
     def OnPawleyLoad(event):
-        sig = lambda Th,U,V,W: 1.17741*math.sqrt(U*tand(Th)**2+V*tand(Th)+W)/100.
-        gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))/100.
-        gamFW = lambda s,g: math.sqrt(s**2+(0.4654996*g)**2)+.5345004*g  #Ubaldo Bafile - private communication
-        choice = data['Histograms'].keys()
-        dlg = wx.SingleChoiceDialog(self,'Select','Powder histogram',choice)
-        if dlg.ShowModal() == wx.ID_OK:
-            histogram = choice[dlg.GetSelection()]
-            Id  = G2gd.GetPatternTreeItemId(self, self.root, histogram)
-            Iparms = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Instrument Parameters'))
-            try:                                            #try single wavelength
-                lam = Iparms[1][Iparms[3].index('Lam')]
-            except ValueError:                              #Ka1 & Ka2 present
-                lam = Iparms[1][Iparms[3].index('Lam1')]
-            GU = Iparms[1][Iparms[3].index('U')]               
-            GV = Iparms[1][Iparms[3].index('V')]               
-            GW = Iparms[1][Iparms[3].index('W')]               
-            LX = Iparms[1][Iparms[3].index('X')]               
-            LY = Iparms[1][Iparms[3].index('Y')]
-        dlg.Destroy()
         generalData = data['General']
         dmin = generalData['Pawley dmin']
@@ -2568 +2643 @@
         A = G2lat.cell2A(cell)
         SGData = generalData['SGData']
-        Laue = SGData['SGLaue']
-        SGLatt = SGData['SGLatt']
-        SGUniq = SGData['SGUniq']        
-        HKLd = G2lat.GenHLaue(dmin,Laue,SGLatt,SGUniq,A)
+        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
         PawleyPeaks = []
         wx.BeginBusyCursor()
         try:
             for h,k,l,d in HKLd:
-                ext,mul = G2spc.GenHKL([h,k,l],SGData)[:2]
+                ext,mul = G2spc.GenHKLf([h,k,l],SGData)[:2]
                 if not ext:
-                    th = asind(lam/(2.0*d))
-                    H = gamFW(sig(th,GU,GV,GW),gam(th,LX,LY))/2.35482
-                    PawleyPeaks.append([h,k,l,mul,2*th,H,False,0,0])
+                    PawleyPeaks.append([h,k,l,mul,d,False,100.0,1.0])
         finally:
             wx.EndBusyCursor()
         data['Pawley ref'] = PawleyPeaks
         FillPawleyReflectionsGrid()
-                
-            
-    def OnPawleyImport(event):
-        dlg = wx.FileDialog(self, 'Choose file with Pawley reflections', '.', '', 
-            'GSAS Pawley files (*.RFL)|*.RFL',wx.OPEN)
-        if self.dirname:
-            dlg.SetDirectory(self.dirname)
-        try:
-            if dlg.ShowModal() == wx.ID_OK:
-                PawleyFile = dlg.GetPath()
-                self.dirname = dlg.GetDirectory()
-                choice = data['Histograms'].keys()
-                dlg2 = wx.SingleChoiceDialog(self,'Select','Powder histogram',choice)
-                if dlg2.ShowModal() == wx.ID_OK:
-                    histogram = choice[dlg2.GetSelection()]
-                    Id  = G2gd.GetPatternTreeItemId(self, self.root, histogram)
-                    Iparms = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Instrument Parameters'))
-                dlg2.Destroy()
-                
-                PawleyPeaks = G2IO.GetPawleyPeaks(PawleyFile)
-                data['Pawley ref'] = PawleyPeaks
-                FillPawleyReflectionsGrid()
-        finally:
-            dlg.Destroy()
-            
+                            
     def OnPawleyDelete(event):
         dlg = wx.MessageDialog(self,'Do you really want to delete Pawley reflections?','Delete', 
@@ -2668 +2714 @@
             self.dataFrame.SetMenuBar(self.dataFrame.PawleyMenu)
             self.dataFrame.Bind(wx.EVT_MENU, OnPawleyLoad, id=G2gd.wxID_PAWLEYLOAD)
-            self.dataFrame.Bind(wx.EVT_MENU, OnPawleyImport, id=G2gd.wxID_PAWLEYIMPORT)
             self.dataFrame.Bind(wx.EVT_MENU, OnPawleyDelete, id=G2gd.wxID_PAWLEYDELETE)            
             FillPawleyReflectionsGrid()            
@@ -2684 +2729 @@
         DData = wx.ScrolledWindow(self.dataDisplay)
         self.dataDisplay.AddPage(DData,'Data')
-        PawleyRefl = G2gd.GSGrid(self.dataDisplay)
-        self.dataDisplay.AddPage(PawleyRefl,'Pawley reflections')
+        self.PawleyRefl = G2gd.GSGrid(self.dataDisplay)
+        self.dataDisplay.AddPage(self.PawleyRefl,'Pawley reflections')
     else:
         DData = wx.ScrolledWindow(self.dataDisplay)

trunk/GSASIIplot.py

@@ -456 +456 @@
                 if self.itemPicked:
                     Page.canvas.SetToolTipString('%9.3f'%(xpos))
-                if self.PickId and self.PatternTree.GetItemText(self.PickId) in ['Index Peak List','Unit Cells List']:
+                if self.PickId:
                     found = []
-                    if len(HKL):
-                        view = Page.toolbar._views.forward()[0][:2]
-                        wid = view[1]-view[0]
-                        found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
-                    if len(found):
-                        h,k,l = found[0][:3] 
-                        Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
-                    else:
-                        Page.canvas.SetToolTipString('')
+                    if self.PatternTree.GetItemText(self.PickId) in ['Index Peak List','Unit Cells List','Reflection Lists']:
+                        if len(HKL):
+                            view = Page.toolbar._views.forward()[0][:2]
+                            wid = view[1]-view[0]
+                            found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
+                        if len(found):
+                            h,k,l = found[0][:3] 
+                            Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
+                        else:
+                            Page.canvas.SetToolTipString('')
 
             except TypeError:
@@ -573 +574 @@
     Ymax = 1.0
     lenX = 0
-    HKL = np.array(self.HKL)
+    if self.PatternTree.GetItemText(PickId) in ['Reflection Lists']:
+        Phases = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId,'Reflection Lists'))
+        HKL = []
+        for peak in Phases[self.RefList]:
+            HKL.append(peak[:6])
+        HKL = np.array(HKL)
+    else:
+        HKL = np.array(self.HKL)
     for Pattern in PlotList:
         xye = Pattern[1]
@@ -659 +667 @@
                 else:
                     Plot.plot(X,Y,colors[N%6],picker=False)
-    if PickId and self.PatternTree.GetItemText(PickId) in ['Index Peak List','Unit Cells List']:
+    if PickId:
         Values,Names = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Instrument Parameters'))[1::2]
         Parms = dict(zip(Names,Values))
@@ -666 +674 @@
         except KeyError:
             wave = Parms['Lam1']
-        peaks = np.array((self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Index Peak List'))))
-        for peak in peaks:
-            if self.qPlot:
-                Plot.axvline(4*np.pi*sind(peak[0]/2.0)/wave,color='b')
-            else:
-                Plot.axvline(peak[0],color='b')
-        for hkl in self.HKL:
-            if self.qPlot:
-                Plot.axvline(4*np.pi*sind(hkl[5]/2.0)/wave,color='r',dashes=(5,5))
-            else:
-                Plot.axvline(hkl[5],color='r',dashes=(5,5))
+        if self.PatternTree.GetItemText(PickId) in ['Index Peak List','Unit Cells List']:
+            peaks = np.array((self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Index Peak List'))))
+            for peak in peaks:
+                if self.qPlot:
+                    Plot.axvline(4*np.pi*sind(peak[0]/2.0)/wave,color='b')
+                else:
+                    Plot.axvline(peak[0],color='b')
+            for hkl in self.HKL:
+                if self.qPlot:
+                    Plot.axvline(4*np.pi*sind(hkl[5]/2.0)/wave,color='r',dashes=(5,5))
+                else:
+                    Plot.axvline(hkl[5],color='r',dashes=(5,5))
+        elif self.PatternTree.GetItemText(PickId) in ['Reflection Lists']:
+            Phases = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId,'Reflection Lists'))
+            for pId,phase in enumerate(Phases):
+                pId += 1
+                peaks = Phases[phase]
+                for peak in peaks:
+                    if self.qPlot:
+                        Plot.plot(2*np.pi/peak[4],offset*N-pId*Ymax*.005,colors[pId%6]+'|',mew=1,ms=8,picker=2.)
+                    else:
+                        Plot.plot(peak[5],offset*N-pId*Ymax*.005,colors[pId%6]+'|',mew=1,ms=8,picker=2.)
+            
     if self.Contour:
         acolor = mpl.cm.get_cmap(self.ContourColor)
@@ -1157 +1177 @@
             elif muStrain[0] == 'uniaxial':
                 def uniaxMustrain(xyz,muiso,muaniso,axes):
-                    cp = abs(np.dot(xyz,axes))
-                    S = muiso+muaniso*cp
-                    return S*xyz*math.pi/0.018      #centidegrees to radians!
+                    Z = np.array(axes)
+                    cp = abs(np.dot(xyz,Z))
+                    sp = np.sqrt(1.-cp**2)
+                    R = muiso*muaniso/np.sqrt((muiso*cp)**2+(muaniso*sp)**2)
+#                    S = muiso+muaniso*cp           #old GSAS - wrong math!!
+                    return R*xyz*math.pi/0.018      #centidegrees to radians!
                 muiso,muaniso = muStrain[1][:2]
                 axes = np.inner(A,np.array(muStrain[3]))
@@ -1176 +1199 @@
                     uvw = np.inner(A.T,xyz)
                     Strm = np.array(G2spc.MustrainCoeff(uvw,SGData))
-                    sum = np.sqrt(np.sum(np.multiply(Shkl,Strm)))*math.pi/0.018      #centidegrees to radians!
+                    sum = np.sum(np.multiply(Shkl,Strm))
+                    sum = np.where(sum > 0.01,sum,0.01)
+                    sum = np.sqrt(sum)*math.pi/0.018      #centidegrees to radians!
                     return sum*xyz
                 Shkl = np.array(muStrain[4])
@@ -1189 +1214 @@
                     
             if np.any(X) and np.any(Y) and np.any(Z):
-                Plot.plot_surface(X,Y,Z,rstride=1,cstride=1,color='g')
+                Plot.plot_surface(X,Y,Z,rstride=1,cstride=1,color='g',linewidth=1)
                 
             Plot.set_xlabel('X')

trunk/GSASIIpwd.py

@@ -92 +92 @@
     return plist
 
-def ValuesOut(parmdict, varylist):
-    '''Use before call to leastsq to setup list of values for the parameters 
-    in parmdict, as selected by key in varylist'''
-    return [parmdict[key] for key in varylist] 
-    
-def ValuesIn(parmdict, varylist, values):
-    ''' Use after call to leastsq to update the parameter dictionary with 
-    values corresponding to keys in varylist'''
-    parmdict.update(zip(varylist,values))
-    
 def Transmission(Geometry,Abs,Diam):
 #Calculate sample transmission
@@ -482 +472 @@
 fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
                 
-def getBackground(parmDict,bakType,xdata):
+def getBackground(pfx,parmDict,bakType,xdata):
+    yb = np.zeros_like(xdata)
     if bakType == 'chebyschev':
-        yb = np.zeros_like(xdata)
         iBak = 0
         while True:
-            key = 'Back'+str(iBak)
+            key = pfx+'Back:'+str(iBak)
             try:
                 yb += parmDict[key]*(xdata-xdata[0])**iBak
                 iBak += 1
             except KeyError:
-                return yb
-                
+                break
+    return yb
+
+def calcPeakFFT(x,fxns,widths,pos,args):
+    dx = x[1]-x[0]
+    z = np.empty([len(fxns),len(x)])
+    for i,(fxn,width,arg) in enumerate(zip(fxns,widths,args)):
+        z[i] = fxn.pdf(x,loc=pos,scale=width,*arg)*dx
+    Z = fft(z)
+    if len(fxns)%2:
+        return ifft(Z.prod(axis=0)).real
+    else:
+        return fftshift(ifft(Z.prod(axis=0))).real
+                    
+def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
+    
+    DX = xdata[1]-xdata[0]
+    widths,fmin,fmax = getWidths(pos,sig,gam,shl)
+    x = np.linspace(pos-fmin,pos+fmin,1024)
+    if pos > 90:
+        fmin,fmax = [fmax,fmin]          
+    dx = x[1]-x[0]
+    arg = [pos,shl/10.,dx,]
+    Df = fcjde.pdf(x,*arg,loc=pos,scale=1.0)
+    if len(np.nonzero(Df)[0])>5:
+        fxns = [st.norm,st.cauchy,fcjde]
+        args = [[],[],arg]
+    else:
+        fxns = [st.norm,st.cauchy]
+        args = [[],[]]
+    Df = calcPeakFFT(x,fxns,widths,pos,args)
+    Df /= np.sum(Df)
+    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
+    return intens*Df(xdata)*DX/dx       #*10 to get close to old fxn???
+                                    
 def getPeakProfile(parmDict,xdata,varyList,bakType):
     
-    def calcPeakFFT(x,fxns,widths,pos,args):
-        dx = x[1]-x[0]
-        z = np.empty([len(fxns),len(x)])
-        for i,(fxn,width,arg) in enumerate(zip(fxns,widths,args)):
-            z[i] = fxn.pdf(x,loc=pos,scale=width,*arg)*dx
-        Z = fft(z)
-        if len(fxns)%2:
-            return ifft(Z.prod(axis=0)).real
-        else:
-            return fftshift(ifft(Z.prod(axis=0))).real
-                        
-    def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
-        
-        DX = xdata[1]-xdata[0]
-        widths,fmin,fmax = getWidths(pos,sig,gam,shl)
-        x = np.linspace(pos-fmin,pos+fmin,1024)
-        if pos > 90:
-            fmin,fmax = [fmax,fmin]          
-        dx = x[1]-x[0]
-        arg = [pos,shl/10.,dx,]
-        Df = fcjde.pdf(x,*arg,loc=pos,scale=1.0)
-        if len(np.nonzero(Df)[0])>5:
-            fxns = [st.norm,st.cauchy,fcjde]
-            args = [[],[],arg]
-        else:
-            fxns = [st.norm,st.cauchy]
-            args = [[],[]]
-        Df = calcPeakFFT(x,fxns,widths,pos,args)
-        Df /= np.sum(Df)
-        Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
-        return intens*Df(xdata)*DX/dx       #*10 to get close to old fxn???
-                    
-    yb = getBackground(parmDict,bakType,xdata)
+    yb = getBackground('',parmDict,bakType,xdata)
     yc = np.zeros_like(yb)
     U = parmDict['U']
@@ -591 +582 @@
     return widths,fmin,fmax
                 
-def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data):
+def Dict2Values(parmdict, varylist):
+    '''Use before call to leastsq to setup list of values for the parameters 
+    in parmdict, as selected by key in varylist'''
+    return [parmdict[key] for key in varylist] 
+    
+def Values2Dict(parmdict, varylist, values):
+    ''' Use after call to leastsq to update the parameter dictionary with 
+    values corresponding to keys in varylist'''
+    parmdict.update(zip(varylist,values))
+    
+def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False):
     
     def SetBackgroundParms(Background):
         bakType,bakFlag = Background[:2]
         backVals = Background[3:]
-        backNames = ['Back'+str(i) for i in range(len(backVals))]
+        backNames = ['Back:'+str(i) for i in range(len(backVals))]
         if bakFlag: #returns backNames as varyList = backNames
             return bakType,dict(zip(backNames,backVals)),backNames
@@ -606 +607 @@
         while True:
             try:
-                bakName = 'Back'+str(iBak)
+                bakName = 'Back:'+str(iBak)
                 Background[iBak+3] = parmList[bakName]
                 iBak += 1
@@ -620 +621 @@
             for i,back in enumerate(Background[3:]):
                 ptstr += ptfmt % (back)
-                sigstr += ptfmt % (sigDict['Back'+str(i)])
+                sigstr += ptfmt % (sigDict['Back:'+str(i)])
             print ptstr
             print sigstr
@@ -675 +676 @@
         print sigstr
 
-    def setPeaksParms(Peaks):
+    def SetPeaksParms(Peaks):
         peakNames = []
         peakVary = []
@@ -734 +735 @@
         return M
         
+    Ftol = 0.01
+    if oneCycle:
+        Ftol = 1.0
     x,y,w,yc,yb,yd = data               #these are numpy arrays!
     xBeg = np.searchsorted(x,Limits[0])
@@ -739 +743 @@
     bakType,bakDict,bakVary = SetBackgroundParms(Background)
     dataType,insDict,insVary = SetInstParms(Inst)
-    peakDict,peakVary = setPeaksParms(Peaks)
+    peakDict,peakVary = SetPeaksParms(Peaks)
     parmDict = {}
     parmDict.update(bakDict)
@@ -747 +751 @@
     while True:
         begin = time.time()
-        values =  np.array(ValuesOut(parmDict, varyList))
+        values =  np.array(Dict2Values(parmDict, varyList))
         if FitPgm == 'LSQ':
             dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0, 
@@ -755 +759 @@
             dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
             try:
-                result = so.leastsq(errPeakProfile,values,
-                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg),full_output=True)
+                result = so.leastsq(errPeakProfile,values,full_output=True,ftol=Ftol,
+                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
             finally:
                 dlg.Destroy()
@@ -762 +766 @@
             chisq = np.sum(result[2]['fvec']**2)
             ncyc = int(result[2]['nfev']/len(varyList))
-            ValuesIn(parmDict, varyList, result[0])
+            Values2Dict(parmDict, varyList, result[0])
             Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
             GOF = chisq/(xFin-xBeg-len(varyList))
@@ -786 +790 @@
         
     sigDict = dict(zip(varyList,sig))
-    yb[xBeg:xFin] = getBackground(parmDict,bakType,x[xBeg:xFin])
+    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
     yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType)
     yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
@@ -913 +917 @@
     msg = 'test '
     gplot = plotter.add('FCJ-Voigt, 11BM').gca()
-    gplot.plot(xdata,getBackground(parmDict0,bakType,xdata))   
+    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
     gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
     fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
-    fplot.plot(xdata,getBackground(parmDict1,bakType,xdata))   
+    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
     fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
+    
+def test1():
+    if NeedTestData: TestData()
     time0 = time.time()
     for i in range(100):
         y = getPeakProfile(parmDict1,xdata,varyList,bakType)
-    print '100+6*Ka1-2 peaks=3600 peaks',time.time()-time0
+    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
+    
     
 

trunk/GSASIIpwdGUI.py

@@ -81 +81 @@
         OnPeakFit('LSQ')
         
+    def OnOneCycle(event):
+        OnPeakFit('LSQ',oneCycle=True)
+        
     def OnBGFSPeakFit(event):
         OnPeakFit('BFGS')
                 
-    def OnPeakFit(FitPgm):
+    def OnPeakFit(FitPgm,oneCycle=False):
         SaveState()
         print 'Peak Fitting:'
@@ -99 +102 @@
         wx.BeginBusyCursor()
         try:
-            G2pwd.DoPeakFit(FitPgm,peaks,background,limits,inst,data)
+            G2pwd.DoPeakFit(FitPgm,peaks,background,limits,inst,data,oneCycle)
         finally:
             wx.EndBusyCursor()    
@@ -237 +240 @@
     self.Bind(wx.EVT_MENU, OnUnDo, id=G2gd.wxID_UNDO)
     self.Bind(wx.EVT_MENU, OnLSQPeakFit, id=G2gd.wxID_LSQPEAKFIT)
+    self.Bind(wx.EVT_MENU, OnOneCycle, id=G2gd.wxID_LSQONECYCLE)
 #    self.Bind(wx.EVT_MENU, OnBGFSPeakFit, id=G2gd.wxID_BFGSPEAKFIT)
     self.Bind(wx.EVT_MENU, OnResetSigGam, id=G2gd.wxID_RESETSIGGAM)
@@ -242 +246 @@
     if data:
         self.dataFrame.PeakFit.Enable(True)
+        self.dataFrame.PFOneCycle.Enable(True)
     self.PickTable = []
     rowLabels = []
@@ -346 +351 @@
     colLabels = ['Tmin','Tmax']
     rowLabels = ['original','changed']
-    Types = [wg.GRID_VALUE_FLOAT+':10,3',wg.GRID_VALUE_FLOAT+':10,3']
+    Types = 2*[wg.GRID_VALUE_FLOAT+':10,3',]
     self.LimitsTable = G2gd.Table(data,rowLabels=rowLabels,colLabels=colLabels,types=Types)
     self.dataFrame.SetLabel('Limits')
@@ -598 +603 @@
         Status = self.dataFrame.CreateStatusBar()    
     self.dataDisplay = wx.Panel(self.dataFrame)
-    
+
+#patch
+    if not 'Gonio. radius' in data:
+        data['Gonio. radius'] = 200.0
+#patch end
+    
+    parms = [['Gonio. radius',' Goniometer radius(mm): ','%.2f',]]
     if data['Type'] == 'Debye-Scherrer':
-        parms = [['DisplaceX',' Sample X displacement: ','%.4f',],
-            ['DisplaceY',' Sample Y displacement: ','%.4f',],
-            ['Absorption',' Sample absorption: ','%.4f',],]
+        parms += [['DisplaceX',' Sample X displacement(\xb5m): ','%.2f',],
+            ['DisplaceY',' Sample Y displacement(\xb5m): ','%.2f',],
+            ['Absorption',' Sample absorption(\xb5r): ','%.4f',],]
     elif data['Type'] == 'Bragg-Brentano':
-        parms = [['Shift',' Sample displacement: ','%.4f',],
-            ['Transparency',' Sample transparency: ','%.4f'],]
-    parms.append(['Temperature',' Sample temperature: ','%.2f'])
-    parms.append(['Pressure',' Sample pressure: ','%.3f'])
-    parms.append(['Humidity',' Sample humidity: ','%.1f'])
-    parms.append(['Voltage',' Sample voltage: ','%.3f'])
-    parms.append(['Force',' Applied load: ','%.3f'])
+        parms += [['Shift',' Sample displacement(\xb5m): ','%.2f',],
+            ['Transparency',' Sample transparency(1/\xb5eff,\xc5): ','%.4f'],]
+    parms.append(['Temperature',' Sample temperature(K): ','%.2f'])
+    parms.append(['Pressure',' Sample pressure(MPa): ','%.3f'])
+    parms.append(['Humidity',' Sample humidity(%): ','%.1f'])
+    parms.append(['Voltage',' Sample voltage(V): ','%.3f'])
+    parms.append(['Force',' Applied load(MN): ','%.3f'])
     objList = {}
 
@@ -657 +668 @@
                 
     mainSizer = wx.BoxSizer(wx.VERTICAL)
-    mainSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' Sample parameters: '),0,wx.ALIGN_CENTER_VERTICAL)
+    mainSizer.Add(wx.StaticText(self.dataDisplay,label=' Sample parameters: '),0,wx.ALIGN_CENTER_VERTICAL)
     mainSizer.Add((5,5),0)
-    scaleSizer = wx.BoxSizer(wx.HORIZONTAL)
+    parmSizer = wx.FlexGridSizer(9,2,5,0)
     scaleRef = wx.CheckBox(self.dataDisplay,label=' Histogram scale factor: ')
     scaleRef.SetValue(data['Scale'][1])
     scaleRef.Bind(wx.EVT_CHECKBOX, OnScaleRef)
-    scaleSizer.Add(scaleRef,0,wx.ALIGN_CENTER_VERTICAL)
+    parmSizer.Add(scaleRef,0,wx.ALIGN_CENTER_VERTICAL)
     scaleVal = wx.TextCtrl(self.dataDisplay,wx.ID_ANY,
         '%.4f'%(data['Scale'][0]),style=wx.TE_PROCESS_ENTER)
     scaleVal.Bind(wx.EVT_TEXT_ENTER,OnScaleVal)
     scaleVal.Bind(wx.EVT_KILL_FOCUS,OnScaleVal)
-    scaleSizer.Add(scaleVal,0,wx.ALIGN_CENTER_VERTICAL)
-    mainSizer.Add(scaleSizer)
-    mainSizer.Add((0,5),0)
+    parmSizer.Add(scaleVal,0,wx.ALIGN_CENTER_VERTICAL)
     typeSizer = wx.BoxSizer(wx.HORIZONTAL)
     choices = ['Debye-Scherrer','Bragg-Brentano',]
@@ -676 +685 @@
         style=wx.CB_READONLY|wx.CB_DROPDOWN)
     histoType.Bind(wx.EVT_COMBOBOX, OnHistoType)
-    typeSizer.Add(histoType)
-    mainSizer.Add(typeSizer)
-    mainSizer.Add((0,5),0)
+    parmSizer.Add(histoType)
+    parmSizer.Add((5,5),0)
     
     for parm in parms:
-        parmSizer = wx.BoxSizer(wx.HORIZONTAL)
         if 'list' in str(type(data[parm[0]])):
             parmRef = wx.CheckBox(self.dataDisplay,label=parm[1])
@@ -698 +705 @@
         parmVal.Bind(wx.EVT_TEXT_ENTER,OnParmVal)
         parmVal.Bind(wx.EVT_KILL_FOCUS,OnParmVal)
-        parmSizer.Add(parmVal,0,wx.ALIGN_CENTER_VERTICAL)
-        mainSizer.Add(parmSizer)
-        mainSizer.Add((0,5),0)    
+        parmSizer.Add(parmVal,1,wx.EXPAND)
+    mainSizer.Add(parmSizer)
+    mainSizer.Add((0,5),0)    
     
     mainSizer.Layout()    
@@ -927 +934 @@
             G2plt.PlotPatterns(self)
         
-        
+    def CopyUnitCell(event):
+        controls,bravais,cells,dmin = self.PatternTree.GetItemPyData(UnitCellsId)
+        for Cell in cells:
+            if Cell[-1]:
+                break
+        cell = Cell[2:9]
+        controls[4] = 1
+        controls[5] = bravaisSymb[cell[0]]
+        controls[6:12] = cell[1:8]
+        controls[12] = G2lat.calc_V(G2lat.cell2A(controls[6:12]))
+        self.PatternTree.SetItemPyData(UnitCellsId,[controls,bravais,cells,dmin])
+        UpdateUnitCellsGrid(self,data)
+        
+        self.dataFrame.RefineCell.Enable(True)
+                
     def RefineCell(event):
         def cellPrint(ibrav,A):
@@ -955 +976 @@
         cell = controls[6:12]
         A = G2lat.cell2A(cell)
-        print controls[5]
         ibrav = bravaisSymb.index(controls[5])
         dmin = G2indx.getDmin(peaks)-0.005
+        self.HKL = G2lat.GenHBravais(dmin,ibrav,A)
+        G2indx.IndexPeaks(peaks,self.HKL)
         if controls[0]:
             Lhkl,M20,X20,Aref,Zero = G2indx.refinePeaksZ(peaks,inst[1],ibrav,A,controls[1])            
@@ -1017 +1039 @@
             UpdateUnitCellsGrid(self,data)
                 
-    def CopyUnitCell(event):
-        controls,bravais,cells,dmin = self.PatternTree.GetItemPyData(UnitCellsId)
-        for Cell in cells:
-            if Cell[-1]:
-                break
-        cell = Cell[2:9]
-        controls[4] = 1
-        controls[5] = bravaisSymb[cell[0]]
-        controls[6:12] = cell[1:8]
-        controls[12] = G2lat.calc_V(G2lat.cell2A(controls[6:12]))
-        self.PatternTree.SetItemPyData(UnitCellsId,[controls,bravais,cells,dmin])
-        UpdateUnitCellsGrid(self,data)
-        self.dataFrame.RefineCell.Enable(True)
-        
-    def MakeNewPhase(event):
-        if not G2gd.GetPatternTreeItemId(self,self.root,'Phases'):
-            sub = self.PatternTree.AppendItem(parent=self.root,text='Phases')
-        else:
-            sub = G2gd.GetPatternTreeItemId(self,self.root,'Phases')
-        PhaseName = ''
-        dlg = wx.TextEntryDialog(None,'Enter a name for this phase','Phase Name Entry','New phase',
-            style=wx.OK)
-        if dlg.ShowModal() == wx.ID_OK:
-            PhaseName = dlg.GetValue()
-        dlg.Destroy()
-        cells = self.PatternTree.GetItemPyData(UnitCellsId)[2]
-        for Cell in cells:
-            if Cell[-1]:
-                break
-        cell = Cell[2:10]        
-        sub = self.PatternTree.AppendItem(parent=sub,text=PhaseName)
-        E,SGData = G2spc.SpcGroup(spaceGroups[cell[0]])
-        self.PatternTree.SetItemPyData(sub, \
-            {'General':{'Name':'phase name','Type':'nuclear','SGData':SGData,
-            'Cell':[False,]+cell[1:],
-            'Pawley dmin':0.25},'Atoms':[],'Drawing':{},'Histograms':{}})
-        Status.SetStatusText('Change space group if needed')
-            
     def RefreshUnitCellsGrid(event):
         cells,dmin = self.PatternTree.GetItemPyData(UnitCellsId)[2:]
@@ -1076 +1060 @@
                     G2plt.PlotPatterns(self)
         
+    def MakeNewPhase(event):
+        if not G2gd.GetPatternTreeItemId(self,self.root,'Phases'):
+            sub = self.PatternTree.AppendItem(parent=self.root,text='Phases')
+        else:
+            sub = G2gd.GetPatternTreeItemId(self,self.root,'Phases')
+        PhaseName = ''
+        dlg = wx.TextEntryDialog(None,'Enter a name for this phase','Phase Name Entry','New phase',
+            style=wx.OK)
+        try:
+            if dlg.ShowModal() == wx.ID_OK:
+                PhaseName = dlg.GetValue()
+                cells = self.PatternTree.GetItemPyData(UnitCellsId)[2]
+                for Cell in cells:
+                    if Cell[-1]:
+                        break
+                cell = Cell[2:10]        
+                sub = self.PatternTree.AppendItem(parent=sub,text=PhaseName)
+                E,SGData = G2spc.SpcGroup(spaceGroups[cell[0]])
+                self.PatternTree.SetItemPyData(sub, \
+                    {'General':{'Name':PhaseName,'Type':'nuclear','SGData':SGData,
+                    'Cell':[False,]+cell[1:],
+                    'Pawley dmin':1.00},'Atoms':[],'Drawing':{},'Histograms':{}})
+                Status.SetStatusText('Change space group if needed')
+        finally:
+            dlg.Destroy()
+            
     if self.dataDisplay:
         self.dataFrame.Clear()
@@ -1170 +1180 @@
     zeroVar.Bind(wx.EVT_CHECKBOX,OnZeroVar)
     littleSizer.Add(zeroVar,0,wx.ALIGN_CENTER_VERTICAL)
-    hklShow = wx.CheckBox(self.dataDisplay,label="  Show starting hkl positions")
+    hklShow = wx.CheckBox(self.dataDisplay,label="  Show hkl positions")
     hklShow.Bind(wx.EVT_CHECKBOX,OnHklShow)
     littleSizer.Add(hklShow,0,wx.ALIGN_CENTER_VERTICAL)
@@ -1217 +1227 @@
         rowLabels = []
         colLabels = ['M20','X20','use','Bravais','a','b','c','alpha','beta','gamma','Volume']
-        Types = [wg.GRID_VALUE_FLOAT+':10,2',wg.GRID_VALUE_NUMBER,wg.GRID_VALUE_BOOL,wg.GRID_VALUE_STRING,
-            wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,5',
-            wg.GRID_VALUE_FLOAT+':10,3',wg.GRID_VALUE_FLOAT+':10,3',wg.GRID_VALUE_FLOAT+':10,3',
-            wg.GRID_VALUE_FLOAT+':10,2']
+        Types = [wg.GRID_VALUE_FLOAT+':10,2',wg.GRID_VALUE_NUMBER,wg.GRID_VALUE_BOOL,wg.GRID_VALUE_STRING,]+ \
+            3*[wg.GRID_VALUE_FLOAT+':10,5',]+3*[wg.GRID_VALUE_FLOAT+':10,3',]+ \
+            [wg.GRID_VALUE_FLOAT+':10,2']
         numRows = len(cells)
         table = []
@@ -1237 +1246 @@
         gridDisplay.SetTable(UnitCellsTable, True)
         self.dataFrame.CopyCell.Enable(True)
-#        gridDisplay.Bind(wg.EVT_GRID_CELL_CHANGE, RefreshUnitCellsGrid)
         gridDisplay.Bind(wg.EVT_GRID_CELL_LEFT_CLICK,RefreshUnitCellsGrid)
         gridDisplay.SetMargins(0,0)
@@ -1249 +1257 @@
                     gridDisplay.SetReadOnly(r,c,isReadOnly=True)
         gridDisplay.SetSize(bottomSize)
+
+def UpdateReflectionGrid(self,data):
+    if not data:
+        print 'No phases, no reflections'
+        return
+    phases = data.keys()
+    
+    def OnSelectPhase(event):
+        dlg = wx.SingleChoiceDialog(self,'Select','Phase',phases)
+        try:
+            if dlg.ShowModal() == wx.ID_OK:
+                sel = dlg.GetSelection()
+                self.RefList = phases[sel]
+                UpdateReflectionGrid(self,data)
+        finally:
+            dlg.Destroy()
+        G2plt.PlotPatterns(self)
+        
+        
+    if self.dataDisplay:
+        self.dataFrame.Clear()
+    self.dataFrame.SetMenuBar(self.dataFrame.ReflMenu)
+    if not self.dataFrame.GetStatusBar():
+        Status = self.dataFrame.CreateStatusBar()    
+    self.dataDisplay = wx.Panel(self.dataFrame)
+    self.Bind(wx.EVT_MENU, OnSelectPhase, id=G2gd.wxID_SELECTPHASE)
+    self.dataFrame.SelectPhase.Enable(False)
+    if len(data) > 1:
+        self.dataFrame.SelectPhase.Enable(True)
+    rowLabels = []
+    refList = []
+    for h,k,l,m,d,pos,sig,gam,fo,fc,x in data[self.RefList]:
+        refList.append([h,k,l,m,d,pos,sig,gam,fo,fc])
+    for i in range(len(refList)): rowLabels.append(str(i+1))
+    colLabels = ['H','K','L','mul','d','pos','sig','gam','Fosq','Fcsq',]
+    Types = 4*[wg.GRID_VALUE_LONG,]+6*[wg.GRID_VALUE_FLOAT+':10,4',]
+    self.PeakTable = G2gd.Table(refList,rowLabels=rowLabels,colLabels=colLabels,types=Types)
+    self.dataFrame.SetLabel('Reflection List for '+self.RefList)
+    self.dataDisplay = G2gd.GSGrid(parent=self.dataFrame)
+    self.dataDisplay.SetTable(self.PeakTable, True)
+    self.dataDisplay.SetMargins(0,0)
+    self.dataDisplay.AutoSizeColumns(False)
+    self.dataFrame.setSizePosLeft([535,350])
 
 def UpdatePDFGrid(self,data):

trunk/GSASIIspc.py

@@ -299 +299 @@
     else:
         return zip(XYZEquiv,Idup,Cell)
+
+def GenHKLf(HKL,SGData,Friedel=False):
+    '''
+    Uses old GSAS Fortran routine genhkl.for
+    input:
+        HKL - [h,k,l]
+        SGData - space group data obtained from SpcGroup
+        Friedel = True to retain Friedel pairs for centrosymmetric case
+    returns:
+        iabsnt = True is reflection is forbidden by symmetry
+        mulp = reflection multiplicity including Fridel pairs
+        Uniq = numpy array of equivalent hkl in descending order of h,k,l
+    '''
+    hklf = HKL+[0,]
+    Ops = SGData['SGOps']
+    OpM = np.array([op[0] for op in Ops])
+    OpT = np.array([op[1] for op in Ops])
+    Inv = SGData['SGInv']
+    Cen = np.array([cen for cen in SGData['SGCen']])
+    
+    Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
+    h,k,l,f = Uniq
+    Uniq=np.array(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq]))
+    Uniq = np.array(Uniq)
+    
+    return iabsnt,2*mulp,Uniq       #include Friedel pairs in powder mulp
+
         
-def GenHKL(HKL,SGData,Friedel=False):
+def GenHKL(HKL,SGData,Friedel=False):          #gives wrong answers!
     '''
     Generate set of equivalent reflections
@@ -312 +339 @@
         Uniq = numpy array of equivalent hkl in descending order of h,k,l
         Phs = numpy array of corresponding phase factors (multiples of 2pi)
+        
+    NB: only works on one HKL at a time - 
+        it can be made to do an array of HKLs - not any faster!
     ''' 
     iabsnt = False
@@ -322 +352 @@
     # get operators & expand if centrosymmetric
     Ops = SGData['SGOps']
-    opM = np.array([op[0].T for op in Ops])
+    opM = np.array([op[0] for op in Ops])
     opT = np.array([op[1] for op in Ops])
     if Inv:
@@ -330 +360 @@
     eqH = np.inner(opM,H)
     HT = np.dot(opT,H)%1.
-    dup = len(ma.nonzero([ma.allclose(HKL,hkl,atol=0.001) for hkl in eqH])[0])
+    dup = len(ma.nonzero([ma.allclose(H,hkl,atol=0.001) for hkl in eqH])[0])
     mulp = len(eqH)/dup
     # generate unique reflection set (with or without Friedel pairs)
@@ -347 +377 @@
     HPT = np.around(HPT-HPT[0],4)%1.
     iabsnt |= np.any(HPT)
-    return iabsnt,mulp,Uniq,Phs
+    return iabsnt,mulp,Uniq
+                                   
+#def GenHKLs(H,SGData,Friedel=False):
+#    '''
+#    Generate set of equivalent reflections
+#    input:
+#        H - np.array of [h,k,l] assumed to exclude lattice centering extinctions
+#        SGData - space group data obtained from SpcGroup
+#        Friedel = True to retain Friedel pairs for centrosymmetric case
+#    returns:
+#        iabsnt = True is reflection is forbidden by symmetry
+#        mulp = reflection multiplicity including Fridel pairs
+#        Uniq = numpy array of equivalent hkl in descending order of h,k,l
+#        Phs = numpy array of corresponding phase factors (multiples of 2pi)
+#    this is not any faster than GenHKL above - oh well
+#    ''' 
+#    H = np.array(H,dtype=float)
+#    Inv = SGData['SGInv']
+#    Ops = SGData['SGOps']
+#    opM = np.array([op[0] for op in Ops])
+#    opT = np.array([op[1] for op in Ops])
+#    if Inv:
+#        opM = np.concatenate((opM,-opM))
+#        opT = np.concatenate((opT,-opT))
+#    # generate full hkl table and phase factors; find duplicates for multiplicity
+#    eqH = np.swapaxes(np.swapaxes(np.inner(opM,H),0,2),1,2)
+#    HeqH = zip(H,eqH)
+#    dup = np.array([sum([ma.allclose(h,hkl,atol=0.001) for hkl in HKL]) for h,HKL in HeqH])
+#    mulp = len(eqH[0])/dup
+#    # generate unique reflection set (with or without Friedel pairs)
+#    HT = np.dot(H,opT.T)%1.
+#    HPT = zip(mulp,dup,[[tuple(hp) for hp in HP] for HP in np.dstack((eqH,HT))])
+#    Uniq = []
+#    Phs = []
+#    iabsnt = []
+#    for mp,dp,hpt in HPT:
+#        hpt.sort()
+#        upt = hpt[::dp]
+#        upt.reverse()
+#        if Inv and not Friedel:
+#            upt = upt[:len(upt)/2]            
+#        uniq = np.split(upt,[3,4],axis=1)
+#        Phs.append(uniq[1].T[0])
+#        Uniq.append(uniq[0])
+#        # determine if extinct
+#        hpt = np.split(hpt,[3,4],axis=1)[1].T[0]
+#        hpt = np.reshape(hpt,(mp,-1)).T
+#        hpt = np.around(hpt-hpt[0],4)%1.
+#        iabsnt.append(np.any(hpt))
+#    return iabsnt,mulp,Uniq,Phs
                                    
 def GetOprPtrName(key):            

trunk/GSASIIstruct.py

@@ -8 +8 @@
 ########### SVN repository information ###################
 import sys
+import os
 import os.path as ospath
 import numpy as np
@@ -13 +14 @@
 import time
 import math
+import wx
 import GSASIIpath
 import GSASIIElem as G2el
 import GSASIIlattice as G2lat
 import GSASIIspc as G2spc
+import GSASIIpwd as G2pwd
+import GSASIImapvars as G2mv
+import scipy.optimize as so
+
+sind = lambda x: np.sin(x*np.pi/180.)
+cosd = lambda x: np.cos(x*np.pi/180.)
+tand = lambda x: np.tan(x*np.pi/180.)
+asind = lambda x: 180.*np.arcsin(x)/np.pi
+
 
 def ShowBanner():
@@ -34 +45 @@
     '''
     file = open(GPXfile,'rb')
-    HistogramNames = []
     while True:
         try:
@@ -82 +92 @@
     return PhaseNames
 
-def GetPhaseData(GPXfile,PhaseName):
-    ''' Returns the 'General' and 'Atoms' objects for PhaseName from GSASII gpx file
+def GetAllPhaseData(GPXfile,PhaseName):
+    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
     input:
         GPXfile = gpx full file name
         PhaseName = phase name
     return:
-        General = dictionary of general phase info.
-        Atoms = list of atom parameters
-        these are returned empty if PhaseName not found
+        phase dictionary
     '''        
     file = open(GPXfile,'rb')
@@ -104 +112 @@
             for datus in data[1:]:
                 if datus[0] == PhaseName:
-                    General = datus[1]['General']
-                    Atoms = datus[1]['Atoms']                                            
+                    break
     file.close()
-    return {'General':General,'Atoms':Atoms}
-    
-def ShowPhaseData(phaseNames,PhaseData):
-    print ' Phases:'
-    for name in phaseNames:
-        General = PhaseData[name]['General']
-        Atoms = PhaseData[name]['Atoms']
-        print '\n Phase name: ',General['Name']
-        SGtext = G2spc.SGPrint(General['SGData'])
-        for line in SGtext: print line
-        cell = General['Cell']
-        print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
-            ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =','%.3f'%(cell[6]),' volume =','%.3f'%(cell[7])
-        print ' Refine?',cell[0]
-        print '\n Atoms:'
-        line = '   name    type  refine?   x         y         z    '+ \
-            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
-        if General['Type'] == 'magnetic':
-            line += '   Mx     My     Mz'
-        elif General['Type'] == 'macromolecular':
-            line = ' res no  residue  chain '+line
-        print line
-        if General['Type'] == 'nuclear':
-            print 135*'-'
-            for at in Atoms:
-                line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
-                    '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
-                if at[9] == 'I':
-                    line += '%8.4f'%(at[10])+48*' '
-                else:
-                    line += 8*' '
-                    for i in range(6):
-                        line += '%8.4f'%(at[11+i])
-                print line
-#        elif General['Type'] == 'magnetic':
-#        elif General['Type'] == 'macromolecular':
-        
-
+    return datus[1]
+    
 def GetHistogramNames(GPXfile):
     ''' Returns a list of histogram names found in GSASII gpx file
@@ -165 +136 @@
     return HistogramNames
     
+def GetUsedHistogramsAndPhases(GPXfile):
+    ''' Returns all histograms that are found in any phase
+    and any phase that uses a histogram
+    input:
+        GPXfile = .gpx full file name
+    return:
+        Histograms = dictionary of histograms as {name:data,...}
+        Phases = dictionary of phases that use histograms
+    '''
+    phaseNames = GetPhaseNames(GPXfile)
+    phaseData = {}
+    for name in phaseNames: 
+        phaseData[name] =  GetAllPhaseData(GPXfile,name)
+    Histograms = {}
+    Phases = {}
+    pId = 0
+    hId = 0
+    for phase in phaseData:
+        Phase = phaseData[phase]
+        if Phase['Histograms']:
+            if phase not in Phases:
+                Phase['pId'] = pId
+                pId += 1
+                Phases[phase] = Phase
+            for hist in Phase['Histograms']:
+                if hist not in Histograms:
+                    if 'PWDR' in hist[:4]: 
+                        Histograms[hist] = GetPWDRdata(GPXfile,hist)
+                    elif 'HKLF' in hist[:4]:
+                        Histograms[hist] = GetHKLFdata(GPXfile,hist)
+                    #future restraint, etc. histograms here            
+                    Histograms[hist]['hId'] = hId
+                    hId += 1
+    return Histograms,Phases
+    
+def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases):
+    ''' Updates gpxfile from all histograms that are found in any phase
+    and any phase that used a histogram
+    input:
+        GPXfile = .gpx full file name
+        Histograms = dictionary of histograms as {name:data,...}
+        Phases = dictionary of phases that use histograms
+    '''
+                        
+    def GPXBackup(GPXfile):
+        import distutils.file_util as dfu
+        GPXpath,GPXname = ospath.split(GPXfile)
+        Name = ospath.splitext(GPXname)[0]
+        files = os.listdir(GPXpath)
+        last = 0
+        for name in files:
+            name = name.split('.')
+            if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
+                last = max(last,int(name[1].strip('bak'))+1)
+        GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
+        dfu.copy_file(GPXfile,GPXback)
+        return GPXback
+        
+    GPXback = GPXBackup(GPXfile)
+    print '\n',130*'-'
+    print 'Read from file:',GPXback
+    print 'Save to file  :',GPXfile
+    infile = open(GPXback,'rb')
+    outfile = open(GPXfile,'wb')
+    while True:
+        try:
+            data = cPickle.load(infile)
+        except EOFError:
+            break
+        datum = data[0]
+        print 'read: ',datum[0]
+        if datum[0] == 'Phases':
+            for iphase in range(len(data)):
+                if data[iphase][0] in Phases:
+                    phaseName = data[iphase][0]
+                    data[iphase][1] = Phases[phaseName]
+        try:
+            histogram = Histograms[datum[0]]
+            print 'found ',datum[0]
+            data[0][1][1] = histogram['Data']
+            for datus in data[1:]:
+                print '    read: ',datus[0]
+                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
+                    datus[1] = histogram[datus[0]]
+        except KeyError:
+            pass
+                                
+        cPickle.dump(data,outfile,1)
+    infile.close()
+    outfile.close()
+    print 'refinement save successful'
+                    
 def GetPWDRdata(GPXfile,PWDRname):
     ''' Returns powder data from GSASII gpx file
@@ -171 +234 @@
         PWDRname = powder histogram name as obtained from GetHistogramNames
     return: 
-        PWDRdata = powder data list:
+        PWDRdata = powder data dictionary with:
+            Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
         
     '''
     file = open(GPXfile,'rb')
-    PWDRdata = []
+    PWDRdata = {}
     while True:
         try:
@@ -183 +247 @@
         datum = data[0]
         if datum[0] == PWDRname:
-            PWDRdata = datum[1:][0]
-                                                
-            
+            PWDRdata['Data'] = datum[1][1]          #powder data arrays
+            PWDRdata[data[2][0]] = data[2][1]       #Limits
+            PWDRdata[data[3][0]] = data[3][1]       #Background
+            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
+            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
+            try:
+                PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
+            except IndexError:
+                PWDRdata['Reflection lists'] = {}
     file.close()
-    return PWDRdata#,Limits,InstrumentParms
+    return PWDRdata
     
 def GetHKLFdata(GPXfile,HKLFname):
@@ -226 +296 @@
                 FFtable[El] = item
     return FFtable
-        
+    
+def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
+    if SGLaue in ['-1','2/m','mmm']:
+        return                      #no Pawley symmetry required constraints
+    for varyI in pawleyVary:
+        refI = int(varyI.split(':')[-1])
+        ih,ik,il = PawleyRef[refI][:3]
+        for varyJ in pawleyVary[:refI]:
+            refJ = int(varyJ.split(':')[-1])
+            jh,jk,jl = PawleyRef[refJ][:3]
+            if SGLaue in ['4/m','4/mmm']:
+                isum = ih**2+ik**2
+                jsum = jh**2+jk**2
+                if abs(il) == abs(jl) and isum == jsum:
+                    G2mv.StoreEquivalence(varyJ,(varyI,))
+            elif SGLaue in ['3R','3mR']:
+                isum = ih**2+ik**2+il**2
+                jsum = jh**2+jk**2*jl**2
+                isum2 = ih*ik+ih*il+ik*il
+                jsum2 = jh*jk+jh*jl+jk*jl
+                if isum == jsum and isum2 == jsum2:
+                    G2mv.StoreEquivalence(varyJ,(varyI,))
+            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
+                isum = ih**2+ik**2+ih*ik
+                jsum = jh**2+jk**2+jh*jk
+                if abs(il) == abs(jl) and isum == jsum:
+                    G2mv.StoreEquivalence(varyJ,(varyI,))
+            elif SGLaue in ['m3','m3m']:
+                isum = ih**2+ik**2+il**2
+                jsum = jh**2+jk**2+jl**2
+                if isum == jsum:
+                    G2mv.StoreEquivalence(varyJ,(varyI,))
+        
+def GetPhaseData(PhaseData):
+    
+    def cellVary(pfx,SGData): 
+        if SGData['SGLaue'] in ['-1',]:
+            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
+        elif SGData['SGLaue'] in ['2/m',]:
+            if SGData['SGUniq'] == 'a':
+                return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
+            elif SGData['SGUniq'] == 'b':
+                return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
+            else:
+                return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
+        elif SGData['SGLaue'] in ['mmm',]:
+            return [pfx+'A0',pfx+'A1',pfx+'A2']
+        elif SGData['SGLaue'] in ['4/m','4/mmm']:
+            return [pfx+'A0',pfx+'A2']
+        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
+            return [pfx+'A0',pfx+'A2']
+        elif SGData['SGLaue'] in ['3R', '3mR']:
+            return [pfx+'A0',pfx+'A3']                       
+        elif SGData['SGLaue'] in ['m3m','m3']:
+            return [pfx+'A0']
+        
+        
+    print ' Phases:'
+    phaseVary = []
+    phaseDict = {}
+    phaseConstr = {}
+    pawleyLookup = {}
+    for name in PhaseData:
+        General = PhaseData[name]['General']
+        pId = PhaseData[name]['pId']
+        pfx = str(pId)+'::'
+        Atoms = PhaseData[name]['Atoms']
+        try:
+            PawleyRef = PhaseData[name]['Pawley ref']
+        except KeyError:
+            PawleyRef = []
+        print '\n Phase name: ',General['Name']
+        SGData = General['SGData']
+        SGtext = G2spc.SGPrint(SGData)
+        for line in SGtext: print line
+        cell = General['Cell']
+        A = G2lat.cell2A(cell[1:7])
+        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]})
+        if cell[0]:
+            phaseVary = cellVary(pfx,SGData)
+        print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
+            ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
+            '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
+        if Atoms:
+            print '\n Atoms:'
+            line = '   name    type  refine?   x         y         z    '+ \
+                '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
+            if General['Type'] == 'magnetic':
+                line += '   Mx     My     Mz'
+            elif General['Type'] == 'macromolecular':
+                line = ' res no  residue  chain '+line
+            print line
+            if General['Type'] == 'nuclear':
+                print 135*'-'
+                for at in Atoms:
+                    line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
+                        '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
+                    if at[9] == 'I':
+                        line += '%8.4f'%(at[10])+48*' '
+                    else:
+                        line += 8*' '
+                        for i in range(6):
+                            line += '%8.4f'%(at[11+i])
+                    print line
+                    if 'X' in at[2]:
+                        xId,xCoef = G2spc.GetCSxinel(at[7])
+                    if 'U' in at[2]:
+                        uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
+#        elif General['Type'] == 'magnetic':
+#        elif General['Type'] == 'macromolecular':
+#       PWDR: harmonics texture, unit cell, atom parameters
+        elif PawleyRef:
+            pawleyVary = []
+            for i,refl in enumerate(PawleyRef):
+                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
+                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
+                if refl[5]:
+                    pawleyVary.append(pfx+'PWLref:'+str(i))
+            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
+            phaseVary += pawleyVary    
+                
+    return phaseVary,phaseDict,pawleyLookup
+    
+def SetPhaseData(parmDict,sigDict,Phases):
+    
+    def cellFill(pfx,SGData,parmDict,sigDict): 
+        if SGData['SGLaue'] in ['-1',]:
+            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
+                parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
+            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
+                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
+        elif SGData['SGLaue'] in ['2/m',]:
+            if SGData['SGUniq'] == 'a':
+                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
+                    parmDict[pfx+'A3'],0,0]
+                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
+                    sigDict[pfx+'A3'],0,0]
+            elif SGData['SGUniq'] == 'b':
+                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
+                    0,parmDict[pfx+'A4'],0]
+                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
+                    0,sigDict[pfx+'A4'],0]
+            else:
+                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
+                    0,0,parmDict[pfx+'A5']]
+                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
+                    0,0,sigDict[pfx+'A5']]
+        elif SGData['SGLaue'] in ['mmm',]:
+            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
+            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
+        elif SGData['SGLaue'] in ['4/m','4/mmm']:
+            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
+            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
+        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
+            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
+                parmDict[pfx+'A0'],0,0]
+            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
+        elif SGData['SGLaue'] in ['3R', '3mR']:
+            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
+                parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
+            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
+        elif SGData['SGLaue'] in ['m3m','m3']:
+            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
+            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
+        return A,sigA
+            
+    print '\n Phases:'
+    for phase in Phases:
+        print ' Result for phase: ',phase
+        Phase = Phases[phase]
+        General = Phase['General']
+        SGData = General['SGData']
+        cell = General['Cell']
+        pId = Phase['pId']
+        pfx = str(pId)+'::'
+        if cell[0]:
+            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
+            print ' Reciprocal metric tensor: '
+            ptfmt = "%15.9f"
+            names = ['A11','A22','A33','A12','A13','A23']
+            namstr = '  names :'
+            ptstr =  '  values:'
+            sigstr = '  esds  :'
+            for name,a,siga in zip(names,A,sigA):
+                namstr += '%15s'%(name)
+                ptstr += ptfmt%(a)
+                if siga:
+                    sigstr += ptfmt%(siga)
+                else:
+                    sigstr += 15*' '
+            print namstr
+            print ptstr
+            print sigstr
+            cell[1:7] = G2lat.A2cell(A)
+            cell[7] = G2lat.calc_V(A)
+            print ' New unit cell: a =%.5f'%(cell[1]),' b =%.5f'%(cell[2]), \
+                ' c =%.5f'%(cell[3]),' alpha =%.3f'%(cell[4]),' beta =%.3f'%(cell[5]), \
+                ' gamma =%.3f'%(cell[6]),' volume =%.3f'%(cell[7])
+        if 'Pawley' in Phase['General']['Type']:
+            pawleyRef = Phase['Pawley ref']
+            for i,refl in enumerate(pawleyRef):
+                key = pfx+'PWLref:'+str(i)
+                refl[6] = parmDict[key]
+                if key in sigDict:
+                    refl[7] = sigDict[key]
+                else:
+                    refl[7] = 0
+
+def GetHistogramPhaseData(Phases,Histograms):
+    
+    def PrintSize(hapData):
+        line = '\n Size model    : '+hapData[0]
+        if hapData[0] in ['isotropic','uniaxial']:
+            line += ' equatorial:'+'%12.2f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
+            if hapData[0] == 'uniaxial':
+                line += ' axial:'+'%12.2f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
+            print line
+        else:
+            print line
+            Snames = ['S11','S22','S33','S12','S13','S23']
+            ptlbls = ' names :'
+            ptstr =  ' values:'
+            varstr = ' refine:'
+            for i,name in enumerate(Snames):
+                ptlbls += '%12s' % (name)
+                ptstr += '%12.6f' % (hapData[4][i])
+                varstr += '%12s' % (str(hapData[5][i]))
+            print ptlbls
+            print ptstr
+            print varstr
+        
+    def PrintMuStrain(hapData,SGData):
+        line = '\n Mustrain model: '+hapData[0]
+        if hapData[0] in ['isotropic','uniaxial']:
+            line += ' equatorial:'+'%12.4f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
+            if hapData[0] == 'uniaxial':
+                line += ' axial:'+'%12.4f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
+            print line
+        else:
+            print line
+            Snames = G2spc.MustrainNames(SGData)
+            ptlbls = ' names :'
+            ptstr =  ' values:'
+            varstr = ' refine:'
+            for i,name in enumerate(Snames):
+                ptlbls += '%12s' % (name)
+                ptstr += '%12.6f' % (hapData[4][i])
+                varstr += '%12s' % (str(hapData[5][i]))
+            print ptlbls
+            print ptstr
+            print varstr
+        
+    
+    hapDict = {}
+    hapVary = []
+    controlDict = {}
+    poType = {}
+    poAxes = {}
+    spAxes = {}
+    spType = {}
+    
+    for phase in Phases:
+        HistoPhase = Phases[phase]['Histograms']
+        SGData = Phases[phase]['General']['SGData']
+        cell = Phases[phase]['General']['Cell'][1:7]
+        A = G2lat.cell2A(cell)
+        pId = Phases[phase]['pId']
+        for histogram in HistoPhase:
+            print '\n Phase: ',phase,' in histogram: ',histogram
+            hapData = HistoPhase[histogram]
+            Histogram = Histograms[histogram]
+            hId = Histogram['hId']
+            limits = Histogram['Limits'][1]
+            inst = Histogram['Instrument Parameters']
+            inst = dict(zip(inst[3],inst[1]))
+            Zero = inst['Zero']
+            if 'C' in inst['Type']:
+                try:
+                    wave = inst['Lam']
+                except KeyError:
+                    wave = inst['Lam1']
+                dmin = wave/(2.0*sind(limits[1]/2.0))
+            pfx = str(pId)+':'+str(hId)+':'
+            for item in ['Scale','Extinction']:
+                hapDict[pfx+item] = hapData[item][0]
+                if hapData[item][1]:
+                    hapVary.append(pfx+item)
+            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
+            if hapData['Pref.Ori.'][0] == 'MD':
+                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
+                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
+                if hapData['Pref.Ori.'][2]:
+                    hapVary.append(pfx+'MD')
+            else:                           #'SH' spherical harmonics
+                for item in hapData['Pref.Ori.'][5]:
+                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
+                    if hapData['Pref.Ori.'][2]:
+                        hapVary.append(pfx+item)                    
+            for item in ['Mustrain','Size']:
+                controlDict[pfx+item+'Type'] = hapData[item][0]
+                if hapData[item][0] in ['isotropic','uniaxial']:
+                    hapDict[pfx+item+':0'] = hapData[item][1][0]
+                    if hapData[item][2][0]:
+                        hapVary.append(pfx+item+':0')
+                    if hapData[item][0] == 'uniaxial':
+                        controlDict[pfx+item+'Axis'] = hapData[item][3]
+                        hapDict[pfx+item+':1'] = hapData[item][1][1]
+                        if hapData[item][2][1]:
+                            hapVary.append(pfx+item+':1')
+                else:       #generalized for mustrain or ellipsoidal for size
+                    if item == 'Mustrain':
+                        names = G2spc.MustrainNames(SGData)
+                        pwrs = []
+                        for name in names:
+                            h,k,l = name[1:]
+                            pwrs.append([int(h),int(k),int(l)])
+                        controlDict[pfx+'MuPwrs'] = pwrs
+                    for i in range(len(hapData[item][4])):
+                        sfx = ':'+str(i)
+                        hapDict[pfx+item+sfx] = hapData[item][4][i]
+                        if hapData[item][5][i]:
+                            hapVary.append(pfx+item+sfx)
+                            
+            print '\n Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
+            print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
+            if hapData['Pref.Ori.'][0] == 'MD':
+                Ax = hapData['Pref.Ori.'][3]
+                print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
+                    ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])                
+            PrintSize(hapData['Size'])
+            PrintMuStrain(hapData['Mustrain'],SGData)
+            HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
+            refList = []
+            for h,k,l,d in HKLd:
+                ext,mul,Uniq = G2spc.GenHKLf([h,k,l],SGData)
+                if ext:
+                    continue
+                if 'C' in inst['Type']:
+                    pos = 2.0*asind(wave/(2.0*d))
+                    if limits[0] < pos < limits[1]:
+                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,Uniq])
+                else:
+                    raise ValueError 
+            Histogram['Reflection Lists'][phase] = refList
+    return hapVary,hapDict,controlDict
+    
+def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms):
+    
+    def PrintSizeAndSig(hapData,sizeSig):
+        line = '\n Size model: '+hapData[0]
+        if hapData[0] in ['isotropic','uniaxial']:
+            line += ' equatorial:%12.2f'%(hapData[1][0])
+            if sizeSig[0][0]:
+                line += ', sig: %8.2f'%(sizeSig[0][0])
+            if hapData[0] == 'uniaxial':
+                line += ' axial:%12.2f'%(hapData[1][1])
+                if sizeSig[0][1]:
+                    line += ', sig: %8.2f'%(sizeSig[0][1])
+            print line
+        else:
+            print line
+            Snames = ['S11','S22','S33','S12','S13','S23']
+            ptlbls = ' name  :'
+            ptstr =  ' value :'
+            sigstr = ' sig   :'
+            for i,name in enumerate(Snames):
+                ptlbls += '%12s' % (name)
+                ptstr += '%12.6f' % (hapData[4][i])
+                if sizeSig[1][i]:
+                    sigstr += '%12.6f' % (sizeSig[1][i])
+                else:
+                    sigstr += 12*' '
+            print ptlbls
+            print ptstr
+            print sigstr
+        
+    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
+        line = '\n Mustrain model: '+hapData[0]
+        if hapData[0] in ['isotropic','uniaxial']:
+            line += ' equatorial:%12.4f'%(hapData[1][0])
+            if mustrainSig[0][0]:
+                line += ',sig: %12.4f'%(mustrainSig[0][0])
+            if hapData[0] == 'uniaxial':
+                line += ' axial:%12.4f'%(hapData[1][1])
+                if mustrainSig[0][1]:
+                     line += ',sig: %12.4f'%(mustrainSig[0][1])
+            print line
+        else:
+            print line
+            Snames = G2spc.MustrainNames(SGData)
+            ptlbls = ' name  :'
+            ptstr =  ' value :'
+            sigstr = ' sig   :'
+            for i,name in enumerate(Snames):
+                ptlbls += '%12s' % (name)
+                ptstr += '%12.6f' % (hapData[4][i])
+                sigstr += '%12.6f' % (mustrainSig[1][i])
+            print ptlbls
+            print ptstr
+            print sigstr
+        
+    for phase in Phases:
+        HistoPhase = Phases[phase]['Histograms']
+        SGData = Phases[phase]['General']['SGData']
+        pId = Phases[phase]['pId']
+        for histogram in HistoPhase:
+            print '\n Phase: ',phase,' in histogram: ',histogram
+            hapData = HistoPhase[histogram]
+            Histogram = Histograms[histogram]
+            hId = Histogram['hId']
+            pfx = str(pId)+':'+str(hId)+':'
+            
+# Add this stuff for Rietveld refinement!!!!            
+#            for item in ['Scale','Extinction']:
+#                hapDict[pfx+item] = hapData[item][0]
+#                if hapData[item][1]:
+#                    hapVary.append(pfx+item)
+#            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
+#            if hapData['Pref.Ori.'][0] == 'MD':
+#                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
+#                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
+#                if hapData['Pref.Ori.'][2]:
+#                    hapVary.append(pfx+'MD')
+#            else:                           #'SH' spherical harmonics
+#                for item in hapData['Pref.Ori.'][5]:
+#                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
+#                    if hapData['Pref.Ori.'][2]:
+#                        hapVary.append(pfx+item)
+#            print '\n Phase fraction  : %10.4f'%(hapData['Scale'][0])
+#            print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0])
+#            if hapData['Pref.Ori.'][0] == 'MD':
+#                Ax = hapData['Pref.Ori.'][3]
+#                print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
+#                    ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2]) 
+               
+            SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
+                'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]]}                  
+            for item in ['Mustrain','Size']:
+                if hapData[item][0] in ['isotropic','uniaxial']:
+                    hapData[item][1][0] = parmDict[pfx+item+':0']
+                    if hapData[item][2][0]: 
+                        SizeMuStrSig[item][0][0] = sigDict[pfx+item+':0']
+                    if hapData[item][0] == 'uniaxial':
+                        hapData[item][1][1] = parmDict[pfx+item+':1']
+                        if hapData[item][2][1]:
+                            SizeMuStrSig[item][0][1] = sigDict[pfx+item+':1']
+                else:       #generalized for mustrain or ellipsoidal for size
+                    for i in range(len(hapData[item][4])):
+                        sfx = ':'+str(i)
+                        hapData[item][4][i] = parmDict[pfx+item+sfx]
+                        if hapData[item][5][i]:
+                            SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
+                            
+            PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
+            PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
+    
+def GetHistogramData(Histograms):
+    
+    def GetBackgroundParms(hId,Background):
+        bakType,bakFlag = Background[:2]
+        backVals = Background[3:]
+        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
+        if bakFlag:                                 #returns backNames as varyList = backNames
+            return bakType,dict(zip(backNames,backVals)),backNames
+        else:                                       #no background varied; varyList = []
+            return bakType,dict(zip(backNames,backVals)),[]
+        
+    def GetInstParms(hId,Inst):
+        insVals,insFlags,insNames = Inst[1:4]
+        dataType = insVals[0]
+        instDict = {}
+        insVary = []
+        pfx = ':'+str(hId)+':'
+        for i,flag in enumerate(insFlags):
+            insName = pfx+insNames[i]
+            instDict[insName] = insVals[i]
+            if flag:
+                insVary.append(insName)
+        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.1)
+        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.1)
+        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.001)
+        return dataType,instDict,insVary
+        
+    def GetSampleParms(hId,Sample):
+        sampVary = []
+        pfx = ':'+str(hId)+':'        
+        sampDict = {pfx+'Gonio. radius':Sample['Gonio. radius']}
+        Type = Sample['Type']
+        if 'Bragg' in Type:             #Bragg-Brentano
+            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
+                sampDict[pfx+item] = Sample[item][0]
+                if Sample[item][1]:
+                    sampVary.append(pfx+item)
+        elif 'Debye' in Type:        #Debye-Scherrer
+            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
+                sampDict[pfx+item] = Sample[item][0]
+                if Sample[item][1]:
+                    sampVary.append(pfx+item)
+        return Type,sampDict,sampVary
+        
+    def PrintBackground(Background):
+        print '\n Background function: ',Background[0],' Refine?',bool(Background[1])
+        line = ' Coefficients: '
+        for back in Background[3:]:
+            line += '%10.3f'%(back)
+        print line 
+        
+    def PrintInstParms(Inst):
+        print '\n Instrument Parameters:'
+        ptlbls = ' name  :'
+        ptstr =  ' value :'
+        varstr = ' refine:'
+        instNames = Inst[3][1:]
+        for i,name in enumerate(instNames):
+            ptlbls += '%12s' % (name)
+            ptstr += '%12.6f' % (Inst[1][i+1])
+            if name in ['Lam1','Lam2','Azimuth']:
+                varstr += 12*' '
+            else:
+                varstr += '%12s' % (str(bool(Inst[2][i+1])))
+        print ptlbls
+        print ptstr
+        print varstr
+        
+    def PrintSampleParms(Sample):
+        print '\n Sample Parameters:'
+        ptlbls = ' name  :'
+        ptstr =  ' value :'
+        varstr = ' refine:'
+        if 'Bragg' in Sample['Type']:
+            for item in ['Scale','Shift','Transparency']:
+                ptlbls += '%14s'%(item)
+                ptstr += '%14.4f'%(Sample[item][0])
+                varstr += '%14s'%(str(bool(Sample[item][1])))
+            
+        elif 'Debye' in Type:        #Debye-Scherrer
+            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
+                ptlbls += '%14s'%(item)
+                ptstr += '%14.4f'%(Sample[item][0])
+                varstr += '%14s'%(str(bool(Sample[item][1])))
+
+        print ptlbls
+        print ptstr
+        print varstr
+        
+
+    histDict = {}
+    histVary = []
+    controlDict = {}
+    for histogram in Histograms:
+        Histogram = Histograms[histogram]
+        hId = Histogram['hId']
+        pfx = ':'+str(hId)+':'
+        controlDict[pfx+'Limits'] = Histogram['Limits'][1]
+        
+        Background = Histogram['Background'][0]
+        Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
+        controlDict[pfx+'bakType'] = Type
+        histDict.update(bakDict)
+        histVary += bakVary
+        
+        Inst = Histogram['Instrument Parameters']
+        Type,instDict,insVary = GetInstParms(hId,Inst)
+        controlDict[pfx+'histType'] = Type
+        histDict.update(instDict)
+        histVary += insVary
+        
+        Sample = Histogram['Sample Parameters']
+        Type,sampDict,sampVary = GetSampleParms(hId,Sample)
+        controlDict[pfx+'instType'] = Type
+        histDict.update(sampDict)
+        histVary += sampVary
+
+        print '\n Histogram: ',histogram,' histogram Id: ',hId
+        print 135*'-'
+        Units = {'C':' deg','T':' msec'}
+        units = Units[controlDict[pfx+'histType'][2]]
+        Limits = controlDict[pfx+'Limits']
+        print ' Instrument type: ',Sample['Type']
+        print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
+        PrintSampleParms(Sample)
+        PrintInstParms(Inst)
+        PrintBackground(Background)
+        
+    return histVary,histDict,controlDict
+    
+def SetHistogramData(parmDict,sigDict,Histograms):
+    
+    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
+        backVals = Background[3:]
+        backSig = [0 for i in range(len(Background[3:]))]
+        for i in range(len(backVals)):
+            backVals[i] = parmDict[pfx+'Back:'+str(i)]
+            if Background[1]:
+                backSig[i] = sigDict[pfx+'Back:'+str(i)]
+        return backSig
+        
+    def SetInstParms(pfx,Inst,parmDict,sigDict):
+        insVals,insFlags,insNames = Inst[1:4]
+        instSig = [0 for i in range(len(insVals))]
+        for i,flag in enumerate(insFlags):
+            insName = pfx+insNames[i]
+            insVals[i] = parmDict[insName]
+            if flag:
+                instSig[i] = sigDict[insName]
+        return instSig
+        
+    def SetSampleParms(pfx,Sample,parmDict,sigDict):
+        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
+            sampSig = [0 for i in range(3)]
+            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
+                Sample[item][0] = parmDict[pfx+item]
+                if Sample[item][1]:
+                    sampSig[i] = sigDict[pfx+item]
+        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
+            sampSig = [0 for i in range(4)]
+            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
+                Sample[item][0] = parmDict[pfx+item]
+                if Sample[item][1]:
+                    sampSig[i] = sigDict[pfx+item]
+        return sampSig
+        
+    def PrintBackgroundSig(Background,backSig):
+        print '\n Background function: ',Background[0]
+        valstr = ' value : '
+        sigstr = ' sig   : '
+        for i,back in enumerate(Background[3:]):
+            valstr += '%10.4f'%(back)
+            if Background[1]:
+                sigstr += '%10.4f'%(backSig[i])
+            else:
+                sigstr += 10*' '
+        print valstr
+        print sigstr 
+        
+    def PrintInstParmsSig(Inst,instSig):
+        print '\n Instrument Parameters:'
+        ptlbls = ' names :'
+        ptstr =  ' value :'
+        sigstr = ' sig   :'
+        instNames = Inst[3][1:]
+        for i,name in enumerate(instNames):
+            ptlbls += '%12s' % (name)
+            ptstr += '%12.6f' % (Inst[1][i+1])
+            if instSig[i+1]:
+                sigstr += '%12.6f' % (instSig[i+1])
+            else:
+                sigstr += 12*' '
+        print ptlbls
+        print ptstr
+        print sigstr
+        
+    def PrintSampleParmsSig(Sample,sampleSig):
+        print '\n Sample Parameters:'
+        ptlbls = ' names :'
+        ptstr =  ' values:'
+        sigstr = ' sig   :'
+        if 'Bragg' in Sample['Type']:
+            for i,item in enumerate(['Scale','Shift','Transparency']):
+                ptlbls += '%14s'%(item)
+                ptstr += '%14.4f'%(Sample[item][0])
+                sigstr += '%14.4f'%(sampleSig[i])
+            
+        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
+            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
+                ptlbls += '%14s'%(item)
+                ptstr += '%14.4f'%(Sample[item][0])
+                sigstr += '%14.4f'%(sampleSig[i])
+
+        print ptlbls
+        print ptstr
+        print sigstr
+        
+    for histogram in Histograms:
+        if 'PWDR' in histogram:
+            Histogram = Histograms[histogram]
+            hId = Histogram['hId']
+            pfx = ':'+str(hId)+':'
+            
+            Background = Histogram['Background'][0]
+            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
+            
+            Inst = Histogram['Instrument Parameters']
+            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
+        
+            Sample = Histogram['Sample Parameters']
+            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
+
+            print '\n Histogram: ',histogram,' histogram Id: ',hId
+            print 135*'-'
+            print ' Instrument type: ',Sample['Type']
+            PrintSampleParmsSig(Sample,sampSig)
+            PrintInstParmsSig(Inst,instSig)
+            PrintBackgroundSig(Background,backSig)
+
+
 def SetupSFcalc(General,Atoms):
     ''' setup data for use in StructureFactor; mostly rearranging arrays
@@ -303 +1068 @@
         print 'Absent'
         
+def Dict2Values(parmdict, varylist):
+    '''Use before call to leastsq to setup list of values for the parameters 
+    in parmdict, as selected by key in varylist'''
+    return [parmdict[key] for key in varylist] 
+    
+def Values2Dict(parmdict, varylist, values):
+    ''' Use after call to leastsq to update the parameter dictionary with 
+    values corresponding to keys in varylist'''
+    parmdict.update(zip(varylist,values))
+    
+def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
+    
+    def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
+        if calcControls[phfx+'SizeType'] == 'isotropic':
+            gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*cosd(refl[5]/2.))
+        elif calcControls[phfx+'SizeType'] == 'uniaxial':
+            H = np.array(refl[:3])
+            P = np.array(calcControls[phfx+'SizeAxis'])
+            cosP,sinP = G2lat.CosSinAngle(H,V,P)
+            gam = (1.8*wave/np.pi)*parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1']
+            gam /= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2)*cosd(refl[5]/2.)
+        else:           #ellipsoidal crystallites
+            H = np.array(refl[:3])
+            gam += 1.8*wave/(np.pi*cosd(refl[5])*np.inner(H,np.inner(sizeEllipse,H)))            
+        if calcControls[phfx+'MustrainType'] == 'isotropic':
+            gam += 0.018*parmDict[phfx+'Mustrain:0']*tand(refl[5]/2.)/np.pi
+        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
+            H = np.array(refl[:3])
+            P = np.array(calcControls[phfx+'MustrainAxis'])
+            cosP,sinP = G2lat.CosSinAngle(H,V,P)
+            Si = parmDict[phfx+'Mustrain:0']
+            Sa = parmDict[phfx+'Mustrain:1']
+            gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
+        else:       #generalized - P.W. Stephens model
+            pwrs = calcControls[phfx+'MuPwrs']
+            sum = 0
+            for i,pwr in enumerate(pwrs):
+                sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
+            gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum
+            
+        return gam  
+        
+        
+    hId = Histogram['hId']
+    hfx = ':%d:'%(hId)
+    bakType = calcControls[hfx+'bakType']
+    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
+    yc = np.zeros_like(yb)
+        
+    if 'C' in calcControls[hfx+'histType']:    
+        U = parmDict[hfx+'U']
+        V = parmDict[hfx+'V']
+        W = parmDict[hfx+'W']
+        X = parmDict[hfx+'X']
+        Y = parmDict[hfx+'Y']
+        shl = max(parmDict[hfx+'SH/L'],0.001)
+        Zero = parmDict[hfx+'Zero']
+        pola = parmDict[hfx+'Polariz.']
+        Ka2 = False
+        if hfx+'Lam1' in parmDict.keys():
+            wave = parmDict[hfx+'Lam1']
+            Ka2 = True
+            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
+            kRatio = parmDict[hfx+'I(L2)/I(L1)']
+        else:
+            wave = parmDict[hfx+'Lam']
+    else:
+        print 'TOF Undefined at present'
+        raise ValueError
+    for phase in Histogram['Reflection Lists']:
+        refList = Histogram['Reflection Lists'][phase]
+        Phase = Phases[phase]
+        pId = Phase['pId']
+        pfx = '%d::'%(pId)
+        phfx = '%d:%d:'%(pId,hId)
+        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
+        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
+        sizeEllipse = []
+        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
+            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
+        for refl in refList:
+            if 'C' in calcControls[hfx+'histType']:
+                h,k,l,m = refl[:4]
+                dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
+                d = np.sqrt(dsq)
+                pos = 2.0*asind(wave/(2.0*d))+Zero
+                const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
+                if 'Bragg' in calcControls[hfx+'instType']:
+                    pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
+                        1.e-7*parmDict[hfx+'Transparency']*sind(pos))            #trans(=1/mueff) in Angstroms
+                else:               #Debye-Scherrer - simple but maybe not right
+                    pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
+                refl[5] = pos
+                tanPos = tand(pos/2.0)
+                refl[6] = U*tanPos**2+V*tanPos+W        #save peak sigma
+                refl[7] = X/cosd(pos/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma
+                if 'Pawley' in Phase['General']['Type']:
+                    try:
+                        refl[8] = parmDict[hfx+'Scale']*m*parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]
+                    except KeyError:
+                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
+                        continue
+                else:
+                    raise ValueError       #wants strctrfacr calc here
+                Wd,fmin,fmax = G2pwd.getWidths(pos,refl[6],refl[7],shl)
+                if pos > 90:
+                    fmin,fmax = [fmax,fmin]          
+                iBeg = np.searchsorted(x,pos-fmin)
+                iFin = np.searchsorted(x,pos+fmax)
+                if not iBeg+iFin:       #peak below low limit - skip peak
+                    continue
+                elif not iBeg-iFin:     #peak above high limit - done
+                    return yc,yb
+                yc[iBeg:iFin] += G2pwd.getFCJVoigt(pos,refl[8],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
+                if Ka2:
+                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
+                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
+                    if pos > 90:
+                        fmin,fmax = [fmax,fmin]          
+                    iBeg = np.searchsorted(x,pos2-fmin)
+                    iFin = np.searchsorted(x,pos2+fmax)
+                    yc[iBeg:iFin] += G2pwd.getFCJVoigt(pos2,refl[8]*kRatio,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
+            else:
+                raise ValueError
+    return yc,yb    
+            
+
+    
 def Refine(GPXfile):
+    
+    def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):        
+        parmdict.update(zip(varylist,values))
+        G2mv.Dict2Map(parmDict)
+        Histograms,Phases = HistoPhases
+        M = np.empty(0)
+        sumwYo = 0
+        Nobs = 0
+        for histogram in Histograms:
+            if 'PWDR' in histogram[:4]:
+                Histogram = Histograms[histogram]
+                hId = Histogram['hId']
+                hfx = ':%d:'%(hId)
+                Limits = calcControls[hfx+'Limits']
+                x,y,w,yc,yb,yd = Histogram['Data']
+                yc *= 0.0                           #zero full calcd profiles
+                yb *= 0.0
+                yd *= 0.0
+                xB = np.searchsorted(x,Limits[0])
+                xF = np.searchsorted(x,Limits[1])
+                Nobs += xF-xB
+                Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
+                sumwYo += Histogram['sumwYo']
+                yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
+                    varylist,Histogram,Phases,calcControls,pawleyLookup)
+                yc[xB:xF] *= parmDict[hfx+'Scale']
+                yc[xB:xF] += yb[xB:xF]
+                yd[xB:xF] = y[xB:xF]-yc[xB:xF]
+                Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
+                M = np.concatenate((M,np.sqrt(w[xB:xF])*(yd[xB:xF])))
+#                print 'sum M^2,y,yb,yc',np.sum(M**2),np.sum(y[xB:xF]),np.sum(yb[xB:xF]),np.sum(yc[xB:xF])
+        Histograms['sumwYo'] = sumwYo
+        Histograms['Nobs'] = Nobs
+        Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
+        if dlg:
+            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
+            if not GoOn:
+                return -M           #abort!!
+        return M
+    
     ShowBanner()
-    Controls = GetControls(GPXfile)
-    ShowControls(Controls)
-    
-    phaseNames = GetPhaseNames(GPXfile)
-    phaseData = {}
-    for name in phaseNames: 
-        phaseData[name] =  GetPhaseData(GPXfile,name)
-    ShowPhaseData(phaseNames,phaseData)
-       
-    histograms = GetHistogramNames(GPXfile)
-    for hist in histograms:
-        if 'PWDR' in hist[:4]: 
-            print hist
-            PWDRdata = GetPWDRdata(GPXfile,hist)
-            print PWDRdata
+    varyList = []
+    parmDict = {}
+    calcControls = {}    
+#    Controls = GetControls(GPXfile)
+#    ShowControls(Controls)            
+    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
+    if not Phases:
+        print ' *** ERROR - you have no phases to refine! ***'
+        print ' *** Refine aborted ***'
+        raise Exception
+    if not Histograms:
+        print ' *** ERROR - you have no data to refine with! ***'
+        print ' *** Refine aborted ***'
+        raise Exception
+    phaseVary,phaseDict,pawleyLookup = GetPhaseData(Phases)
+    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
+    calcControls.update(controlDict)
+    histVary,histDict,controlDict = GetHistogramData(Histograms)
+    calcControls.update(controlDict)
+    varyList = phaseVary+histVary+hapVary
+    parmDict.update(phaseDict)
+    parmDict.update(hapDict)
+    parmDict.update(histDict)
+    constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
+    groups,parmlist = G2mv.GroupConstraints(constrDict)
+    G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList)
+    G2mv.Map2Dict(parmDict,varyList)
+
+    while True:
+        begin = time.time()
+        values =  np.array(Dict2Values(parmDict, varyList))
+        dlg = wx.ProgressDialog('Residual','Powder profile Rwp =',101.0, 
+            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
+        screenSize = wx.ClientDisplayRect()
+        Size = dlg.GetSize()
+        dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
+        try:
+            result = so.leastsq(errRefine,values,full_output=True,  #factor=1.,epsfcn=0.00001,ftol=0.0001,
+                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
+            print len(result)
+        finally:
+            dlg.Destroy()
+        runtime = time.time()-begin    
+        chisq = np.sum(result[2]['fvec']**2)
+        ncyc = int(result[2]['nfev']/len(varyList))
+        Values2Dict(parmDict, varyList, result[0])
+        G2mv.Dict2Map(parmDict)
+        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
+        GOF = chisq/(Histograms['Nobs']-len(varyList))
+        print '\n Refinement results:'
+        print 135*'-'
+        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
+        print 'Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
+        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
+        try:
+            sig = np.sqrt(np.diag(result[1])*GOF)
+            if np.any(np.isnan(sig)):
+                print '*** Least squares aborted - some invalid esds possible ***'
+            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
+#            for item in table: print item,table[item]               #useful debug - are things shifting?
+            break                   #refinement succeeded - finish up!
+        except ValueError:          #result[1] is None on singular matrix
+            print '**** Refinement failed - singular matrix ****'
+            Ipvt = result[2]['ipvt']
+            for i,ipvt in enumerate(Ipvt):
+                if not np.sum(result[2]['fjac'],axis=1)[i]:
+                    print 'Removing parameter: ',varyList[ipvt-1]
+                    del(varyList[ipvt-1])
+                    break
+
+    sigDict = dict(zip(varyList,sig))
+    SetPhaseData(parmDict,sigDict,Phases)
+    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
+    SetHistogramData(parmDict,sigDict,Histograms)
+    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases)
 
 def main():