Changeset 342


Ignore:
Timestamp:
Aug 5, 2011 2:35:43 PM (10 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:
1 added
12 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASII.py

    r336 r342  
    147147        self.Refine = parent.Append(help='', id=wxID_REFINE, kind=wx.ITEM_NORMAL,
    148148            text='Refine')
    149         self.Refine.Enable(False)
     149        self.Refine.Enable(True)
    150150        self.Bind(wx.EVT_MENU, self.OnRefine, id=wxID_REFINE)
    151151        self.Solve = parent.Append(help='', id=wxID_SOLVE, kind=wx.ITEM_NORMAL,
     
    391391                            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Peak List'),[])
    392392                            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
    393                             self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])             
     393                            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])
     394                            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Reflection Lists'),{})             
    394395                            self.PatternId = G2gd.GetPatternTreeItemId(self,Id,'Limits')
    395396                    finally:
     
    746747                    Id = self.PatternTree.AppendItem(parent=self.root,text=outname)
    747748                    if Id:
    748                         Sample = {'Scale':[1.0,True],'Type':'Debye-Scherrer','Absorption':[0.0,False],'DisplaceX':[0.0,False],
    749                             'DisplaceY':[0.0,False],'Diffuse':[],'Temperature':300.,'Pressure':1.0,'Humidity':0.0,'Voltage':0.0,'Force':0.0}
     749                        Sample = {'Scale':[1.0,True],'Type':'Debye-Scherrer','Absorption':[0.0,False],
     750                            'DisplaceX':[0.0,False],'DisplaceY':[0.0,False],'Diffuse':[],
     751                            'Temperature':300.,'Pressure':1.0,'Humidity':0.0,
     752                            'Voltage':0.0,'Force':0.0,'Gonio. radius':200.0}
    750753                        self.PatternTree.SetItemPyData(Id,[[''],[Xsum,Ysum,Wsum,YCsum,YBsum,YDsum]])
    751754                        self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Comments'),Comments)                   
     
    757760                        self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
    758761                        self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])             
     762                        self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Reflection Lists'),{})             
    759763                        self.PatternTree.SelectItem(Id)
    760764                        self.PatternTree.Expand(Id)
     
    878882        sub = self.PatternTree.AppendItem(parent=sub,text=PhaseName)
    879883        E,SGData = G2spc.SpcGroup('P 1')
    880         self.PatternTree.SetItemPyData(sub, \
    881             {'General':{'Name':PhaseName,'Type':'nuclear','SGData':SGData,
    882             'Cell':[False,10.,10.,10.,90.,90.,90,1000.],
    883             'Pawley dmin':1.0},'Atoms':[],'Drawing':{},'Histograms':{},'Pawley ref':[],
    884             'Models':{},'SH Texture':{'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
    885             'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
    886             'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}})
     884        self.PatternTree.SetItemPyData(sub, {
     885            'General':{
     886                'Name':PhaseName,
     887                'Type':'nuclear',
     888                'SGData':SGData,
     889                'Cell':[False,10.,10.,10.,90.,90.,90,1000.],
     890                'Pawley dmin':1.0},
     891            'Atoms':[],
     892            'Drawing':{},
     893            'Histograms':{},
     894            'Pawley ref':[],
     895            'Models':{},
     896            'SH Texture':{
     897                'Order':0,
     898                'Model':'cylindrical',
     899                'Sample omega':[False,0.0],
     900                'Sample chi':[False,0.0],
     901                'Sample phi':[False,0.0],
     902                'SH Coeff':[False,{}],
     903                'SHShow':False,
     904                'PFhkl':[0,0,1],
     905                'PFxyz':[0,0,1],
     906                'PlotType':'Pole figure'}
     907            })
    887908       
    888909    def OnDeletePhase(self,event):
     910        #Hmm, also need to delete this phase from Reflection Lists for each PWDR histogram
    889911        if self.dataFrame:
    890912            self.dataFrame.Clear()
     
    13231345       
    13241346    def OnRefine(self,event):
     1347        self.OnFileSave(event)
    13251348        #works - but it'd be better if it could restore plots
    13261349        G2str.Refine(self.GSASprojectfile)
  • trunk/GSASIIIO.py

    r303 r342  
    218218        sig = float(item[6])/rt2ln2x2
    219219        Iobs = float(item[7])*mult
    220         PawleyPeaks.append([h,k,l,mult,tth,sig,False,Iobs,0.0,[]])
     220        PawleyPeaks.append([h,k,l,mult,tth,False,Iobs,0.0])
    221221        S = File.readline()
    222222        item = S.split()
     
    801801            for datus in data[1:]:
    802802                print '    load: ',datus[0]
    803                                
    804803                sub = self.PatternTree.AppendItem(Id,datus[0])
    805804                self.PatternTree.SetItemPyData(sub,datus[1])
    806                                                
    807805            if 'IMG' in datum[0]:                   #retreive image default flag & data if set
    808806                Data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Image Controls'))
     
    10171015    textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
    10181016        'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
    1019         'SHShow':False,'PFhkl':[0,0,1]}
     1017        'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}
    10201018    shNcof = 0
    10211019    file = open(filename, 'Ur')
    1022     Phase = {}
    10231020    S = 1
    10241021    Expr = [{},{},{},{},{},{},{},{},{}]
     
    11211118                shCoef[key] = float(val)
    11221119        textureData['SH Coeff'] = [False,shCoef]
     1120       
     1121    Phase = {
     1122            'General':{
     1123                'Name':PhaseName,
     1124                'Type':Ptype,
     1125                'SGData':SGData,
     1126                'Cell':[False,]+abc+angles+[Volume,],
     1127                'Pawley dmin':1.0},
     1128            'Atoms':Atoms,
     1129            'Drawing':{},
     1130            'Histograms':{},
     1131            'Pawley ref':[],
     1132            'Models':{},
     1133            'SH Texture':textureData
     1134            }
    11231135           
    1124     Phase['General'] = {'Name':PhaseName,'Type':Ptype,'SGData':SGData,
    1125         'Cell':[False,]+abc+angles+[Volume,],'Pawley ref':[],'Models':{},'SH Texture':textureData}
    1126     Phase['Atoms'] = Atoms
    1127     Phase['Drawing'] = {}
    1128     Phase['Histograms'] = {}
    11291136    return Phase
    11301137       
  • trunk/GSASIIgrid.py

    r335 r342  
    4646] = [wx.NewId() for _init_coll_IndPeaks_Items in range(1)]
    4747
    48 [ wxID_UNDO,wxID_LSQPEAKFIT,wxID_BFGSPEAKFIT,wxID_RESETSIGGAM,
    49 ] = [wx.NewId() for _init_coll_PEAK_Items in range(4)]
     48[ wxID_UNDO,wxID_LSQPEAKFIT,wxID_LSQONECYCLE,wxID_BFGSPEAKFIT,wxID_RESETSIGGAM,
     49] = [wx.NewId() for _init_coll_PEAK_Items in range(5)]
    5050
    5151[  wxID_INDEXPEAKS, wxID_REFINECELL, wxID_COPYCELL, wxID_MAKENEWPHASE,
    5252] = [wx.NewId() for _init_coll_INDEX_Items in range(4)]
     53
     54[ wxID_SELECTPHASE,
     55] = [wx.NewId() for _init_coll_Refl_Items in range(1)]
    5356
    5457[ wxID_PDFCOPYCONTROLS, wxID_PDFSAVECONTROLS, wxID_PDFLOADCONTROLS,
     
    9295        parent.Append(menu=self.IndexEdit, title='Cell Index/Refine')
    9396       
     97    def _init_coll_ReflMenu(self,parent):
     98        parent.Append(menu=self.ReflEdit, title='Reflection List')
     99
    94100    def _init_coll_PDFMenu(self,parent):
    95101        parent.Append(menu=self.PDFEdit, title='PDF Controls')
     
    148154        parent.Append(id=wxID_PAWLEYLOAD, kind=wx.ITEM_NORMAL,text='Pawley create',
    149155            help='Initialize Pawley reflection list')
    150         parent.Append(id=wxID_PAWLEYIMPORT, kind=wx.ITEM_NORMAL,text='Pawley import',
    151             help='Import Pawley reflection list')
     156#        parent.Append(id=wxID_PAWLEYIMPORT, kind=wx.ITEM_NORMAL,text='Pawley import',
     157#            help='Import Pawley reflection list')
    152158        parent.Append(id=wxID_PAWLEYDELETE, kind=wx.ITEM_NORMAL,text='Pawley delete',
    153159            help='Delete Pawley reflection list')
     
    156162        parent.Append(help='Load/Reload index peaks from peak list',id=wxID_INDXRELOAD,
    157163            kind=wx.ITEM_NORMAL,text='Load/Reload')
     164           
     165    def _init_coll_Refl_Items(self,parent):
     166        self.SelectPhase = parent.Append(help='Select phase for reflection list',id=wxID_SELECTPHASE,
     167            kind=wx.ITEM_NORMAL,text='Select phase')
    158168           
    159169    def _init_coll_Image_Items(self,parent):
     
    191201        self.PeakFit = parent.Append(id=wxID_LSQPEAKFIT, kind=wx.ITEM_NORMAL,text='LSQ PeakFit',
    192202            help='Peak fitting via least-squares' )
     203        self.PFOneCycle = parent.Append(id=wxID_LSQONECYCLE, kind=wx.ITEM_NORMAL,text='LSQ one cycle',
     204            help='One cycle of Peak fitting via least-squares' )
    193205#        self.PeakFit = parent.Append(id=wxID_BFGSPEAKFIT, kind=wx.ITEM_NORMAL,text='BFGS PeakFit',
    194206#            help='Peak fitting via BFGS algorithm' )
     
    236248        self.IndPeaksMenu = wx.MenuBar()
    237249        self.IndexMenu = wx.MenuBar()
     250        self.ReflMenu = wx.MenuBar()
    238251        self.PDFMenu = wx.MenuBar()
    239252        self.AtomEdit = wx.Menu(title='')
     
    247260        self.IndPeaksEdit = wx.Menu(title='')
    248261        self.IndexEdit = wx.Menu(title='')
     262        self.ReflEdit = wx.Menu(title='')
    249263        self.PDFEdit = wx.Menu(title='')
    250264        self._init_coll_AtomsMenu(self.AtomsMenu)
     
    268282        self._init_coll_IndexMenu(self.IndexMenu)
    269283        self._init_coll_Index_Items(self.IndexEdit)
     284        self._init_coll_ReflMenu(self.ReflMenu)
     285        self._init_coll_Refl_Items(self.ReflEdit)
    270286        self._init_coll_PDFMenu(self.PDFMenu)
    271287        self._init_coll_PDF_Items(self.PDFEdit)
    272288        self.UnDo.Enable(False)
    273289        self.PeakFit.Enable(False)
     290        self.PFOneCycle.Enable(False)
    274291        self.IndexPeaks.Enable(False)
    275292        self.CopyCell.Enable(False)
     
    879896        else:
    880897            G2plt.PlotPatterns(self)
     898    elif self.PatternTree.GetItemText(item) == 'Reflection Lists':
     899        self.PatternId = self.PatternTree.GetItemParent(item)
     900        self.PickId = item
     901        data = self.PatternTree.GetItemPyData(item)
     902        self.RefList = ''
     903        if len(data):
     904            self.RefList = data.keys()[0]
     905        G2pdG.UpdateReflectionGrid(self,data)
     906        G2plt.PlotPatterns(self)
    881907     
  • trunk/GSASIIimgGUI.py

    r296 r342  
    88########### SVN repository information ###################
    99import wx
    10 import wx.grid as wg
    1110import matplotlib as mpl
    1211import math
  • trunk/GSASIIindex.py

    r312 r342  
    268268   
    269269    values = A2values(ibrav,A)
    270     result = so.leastsq(errFit,values,args=(ibrav,Peaks[7],Peaks[4:7],Pwr),
    271         full_output=True,factor=0.1)
     270    result = so.leastsq(errFit,values,full_output=True,ftol=0.0001,
     271        args=(ibrav,Peaks[7],Peaks[4:7],Pwr))
    272272    A = Values2A(ibrav,result[0])
    273273    return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Pwr)**2),A,result
     
    286286    values = A2values(ibrav,A)
    287287    values.append(Z)
    288     result = so.leastsq(errFit,values,args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Pwr),
    289         full_output=True)
     288    result = so.leastsq(errFit,values,full_output=True,ftol=0.001,
     289        args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Pwr))
    290290    A = Values2A(ibrav,result[0][:-1])
    291291    Z = result[0][-1]
     
    348348    dmin = getDmin(peaks)
    349349    smin = 1.0e10
    350     pwr = 4
     350    pwr = 8
    351351    maxTries = 10
    352352    OK = False
     
    380380        Peaks = np.array(peaks).T
    381381        H = Peaks[4:7]
    382         Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A))
     382        try:
     383            Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A))
     384        except FloatingPointError:
     385            print G2lat.calc_rDsq(H,A)
     386            Peaks[8] = 1.0
    383387        peaks = Peaks.T
    384388       
  • trunk/GSASIIlattice.py

    r326 r342  
    162162    B = nl.inv(A)
    163163    return A,B
     164   
     165def U6toUij(U6):
     166    '''Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23]
     167    returns
     168    '''
     169    U = np.zeros(shape=(3,3))
     170    U = [
     171        [U6[0],  U6[3],  U6[4]],
     172        [U6[3],  U6[1],  U6[5]],
     173        [U6[4],  U6[5],  U6[2]]]
     174    return U
    164175       
    165176def Uij2betaij(Uij,G):
     
    173184    '''
    174185    pass
     186   
     187def CosSinAngle(U,V,G):
     188    ''' calculate sin & cos of angle betwee U & V in generalized coordinates
     189    defined by metric tensor G
     190    input:
     191        U & V - 3-vectors assume numpy arrays
     192        G - metric tensor for U & V defined space assume numpy array
     193    return:
     194        cos(phi) & sin(phi)
     195    '''
     196    u = U/nl.norm(U)
     197    v = V/nl.norm(V)
     198    cosP = np.inner(u,np.inner(G,v))
     199    sinP = np.sqrt(1.0-cosP**2)
     200    return cosP,sinP
    175201   
    176202def criticalEllipse(prob):
     
    485511    return sortHKLd(HKL,True,False)
    486512   
    487 def GenHLaue(dmin,SGLaue,SGLatt,SGUniq,A):
     513def GenHLaue(dmin,SGData,A):
    488514    '''Generate the crystallographically unique powder diffraction reflections
    489515    for a lattice and Bravais type
    490     '''
    491 # dmin - minimum d-spacing
    492 # SGLaue - Laue group symbol = '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm',
    493 #                            '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
    494 # SGLatt - lattice centering = 'P','A','B','C','I','F'
    495 # SGUniq - code for unique monoclinic axis = 'a','b','c'
    496 # A - 6 terms as defined in calc_rDsq
    497 # returns - HKL = list of [h,k,l,d] sorted with largest d first and is unique
    498 # part of reciprocal space ignoring anomalous dispersion
     516    Input:
     517        dmin - minimum d-spacing
     518        SGData - space group dictionary with at least:
     519            SGLaue - Laue group symbol = '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm',
     520                     '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
     521            SGLatt - lattice centering = 'P','A','B','C','I','F'
     522            SGUniq - code for unique monoclinic axis = 'a','b','c'
     523        A - 6 terms as defined in calc_rDsq
     524    Return;
     525        HKL = list of [h,k,l,d] sorted with largest d first and is unique
     526            part of reciprocal space ignoring anomalous dispersion
     527    '''
    499528    import math
     529    SGLaue = SGData['SGLaue']
     530    SGLatt = SGData['SGLatt']
     531    SGUniq = SGData['SGUniq']
    500532    #finds maximum allowed hkl for given A within dmin
    501533    if SGLaue in ['3R','3mR']:        #Rhombohedral axes
     
    645677                    if OdfChk(SGLaue,iord,n):
    646678                        coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
     679            else:
     680                for n in [i-iord for i in range(2*iord+1)]:
     681                    if OdfChk(SGLaue,iord,n):
     682                        coeffNames.append('C(%d,%d)'%(iord,n))
    647683    return coeffNames
    648684
     
    11301166        system = spdict['SGSys']
    11311167
    1132         g2list = GenHLaue(dmin,Laue,center,Axis,cell2A(cell))
     1168        g2list = GenHLaue(dmin,spdict,cell2A(cell))
    11331169        #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
    11341170        #    print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
  • trunk/GSASIIphsGUI.py

    r335 r342  
    197197            self.dataFrame.SetLabel('Phase Data for '+generalData['Name'])
    198198            self.PatternTree.SetItemText(Item,generalData['Name'])
     199            #Hmm, need to change phase name key in Reflection Lists for each histogram
    199200                       
    200201        def OnPhaseType(event):
     
    208209                        self.dataDisplay.DeletePage(self.dataDisplay.FindPage('Draw Options'))
    209210                        self.dataDisplay.DeletePage(self.dataDisplay.FindPage('Draw Atoms'))
    210                         self.dataDisplay.AdvanceSelection()
    211                     if not self.dataDisplay.FindPage('Pawley reflections'):     
    212                         self.dataDisplay.AddPage(G2gd.GSGrid(self.dataDisplay),'Pawley reflections')
     211                    if not self.dataDisplay.FindPage('Pawley reflections'):
     212                        self.PawleyRefl = G2gd.GSGrid(self.dataDisplay)     
     213                        self.dataDisplay.AddPage(self.PawleyRefl,'Pawley reflections')
    213214            else:
    214                 TypeTxt.SetValue(generalData['Type'])
    215                
     215                TypeTxt.SetValue(generalData['Type'])               
    216216               
    217217        def OnSpaceGroup(event):
     
    482482           
    483483        AAchoice = ": ,ALA,ARG,ASN,ASP,CYS,GLN,GLU,GLY,HIS,ILE,LEU,LYS,MET,PHE,PRO,SER,THR,TRP,TYR,VAL,MSE,HOH,UNK"
    484         Types = [wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,wg.GRID_VALUE_CHOICE+": ,X,XU,U,F,FX,FXU,FU",
    485             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
     484        Types = [wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,wg.GRID_VALUE_CHOICE+": ,X,XU,U,F,FX,FXU,FU",]+ \
     485            3*[wg.GRID_VALUE_FLOAT+':10,5',]+[wg.GRID_VALUE_FLOAT+':10,4', #x,y,z,frac
    486486            wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,wg.GRID_VALUE_CHOICE+":I,A",]
    487487        Types += 7*[wg.GRID_VALUE_FLOAT+':10,4',]
     
    490490            colLabels += ['Mx','My','Mz']
    491491            Types[2] = wg.GRID_VALUE_CHOICE+": ,X,XU,U,M,MX,MXU,MU,F,FX,FXU,FU,FM,FMX,FMU,"
    492             Types += [
    493                 wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_FLOAT+':10,4']
     492            Types += 3*[wg.GRID_VALUE_FLOAT+':10,4',]
    494493        elif generalData['Type'] == 'macromolecular':
    495494            colLabels = ['res no','residue','chain'] + colLabels
     
    622621        def ChangeAtomCell(event):
    623622           
    624             def chkUij(Uij,CSI):
     623            def chkUij(Uij,CSI): #needs to do something!!!
    625624                return Uij
    626625
     
    10311030        cx,ct,cs = drawingData['atomPtrs']
    10321031        atomData = drawingData['Atoms']
    1033         Types = [wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,
    1034             wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,5',    #x,y,z
    1035             wg.GRID_VALUE_STRING,wg.GRID_VALUE_CHOICE+": ,lines,vdW balls,sticks,balls & sticks,ellipsoids,polyhedra",
     1032        Types = [wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,5',]+ \
     1033            [wg.GRID_VALUE_STRING,wg.GRID_VALUE_CHOICE+": ,lines,vdW balls,sticks,balls & sticks,ellipsoids,polyhedra",
    10361034            wg.GRID_VALUE_CHOICE+": ,type,name,number",wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,]
    10371035        styleChoice = [' ','lines','vdW balls','sticks','balls & sticks','ellipsoids','polyhedra']
     
    10401038        if generalData['Type'] == 'macromolecular':
    10411039            colLabels = ['Residue','1-letter','Chain'] + colLabels
    1042             Types = [wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING]+Types
     1040            Types = 3*[wg.GRID_VALUE_STRING,]+Types
    10431041            Types[8] = wg.GRID_VALUE_CHOICE+": ,lines,vdW balls,sticks,balls & sticks,ellipsoids,backbone,ribbons,schematic"
    10441042            styleChoice = [' ','lines','vdW balls','sticks','balls & sticks','ellipsoids','backbone','ribbons','schematic']
     
    19901988            Obj = event.GetEventObject()
    19911989            hist,pid = Indx[Obj.GetId()]
    1992             UseList[hist]['Size'][2][pid] = Obj.GetValue()
     1990            if UseList[hist]['Size'][0] == 'ellipsoidal':
     1991                UseList[hist]['Size'][5][pid] = Obj.GetValue()               
     1992            else:
     1993                UseList[hist]['Size'][2][pid] = Obj.GetValue()
    19931994           
    19941995        def OnSizeVal(event):
    19951996            Obj = event.GetEventObject()
    19961997            hist,pid = Indx[Obj.GetId()]
    1997             try:
    1998                 size = float(Obj.GetValue())
    1999                 if pid == 0 and size <= 0:
    2000                     raise ValueError
    2001                 elif pid == 1 and size <= -UseList[hist]['Size'][1][0]:
    2002                     raise ValueError
    2003                 UseList[hist]['Size'][1][pid] = size
    2004             except ValueError:
    2005                 pass
     1998            if UseList[hist]['Size'][0] == 'ellipsoidal':
     1999                try:
     2000                    size = float(Obj.GetValue())
     2001                    if pid < 3 and size <= 0.01:            #10A lower limit!
     2002                        raise ValueError                   
     2003                except ValueError:
     2004                    pass
     2005            else:
     2006                try:
     2007                    size = float(Obj.GetValue())
     2008                    if size <= 0.01:            #10A lower limit!
     2009                        raise ValueError
     2010                    UseList[hist]['Size'][1][pid] = size
     2011                except ValueError:
     2012                    pass
    20062013            Obj.SetValue("%.1f"%(UseList[hist]['Size'][1][pid]))          #reset in case of error
    20072014           
     
    20462053                    UseList[hist]['Mustrain'][4][pid] = strain
    20472054                else:
    2048                     if pid == 0 and strain < 0:
    2049                         raise ValueError
    2050                     elif pid == 1 and strain < -UseList[hist]['Mustrain'][1][0]:
     2055                    if strain <= 0:
    20512056                        raise ValueError
    20522057                    UseList[hist]['Mustrain'][1][pid] = strain
     
    20542059                pass
    20552060            if UseList[hist]['Mustrain'][0] == 'generalized':
    2056                 Obj.SetValue("%.5f"%(UseList[hist]['Mustrain'][4][pid]))          #reset in case of error
     2061                Obj.SetValue("%.3f"%(UseList[hist]['Mustrain'][4][pid]))          #reset in case of error
    20572062            else:
    2058                 Obj.SetValue("%.5f"%(UseList[hist]['Mustrain'][1][pid]))          #reset in case of error
     2063                Obj.SetValue("%.3f"%(UseList[hist]['Mustrain'][1][pid]))          #reset in case of error
    20592064            G2plt.PlotStrain(self,data)
    20602065           
     
    20732078            Obj.SetValue('%3d %3d %3d'%(h,k,l))
    20742079            G2plt.PlotStrain(self,data)
    2075 
    2076         def OnMDRef(event):
     2080           
     2081        def OnPOType(event):
    20772082            Obj = event.GetEventObject()
    20782083            hist = Indx[Obj.GetId()]
    2079             UseList[hist]['MDtexture'][1] = Obj.GetValue()
    2080            
    2081         def OnMDVal(event):
     2084            if 'March' in POType.GetValue():
     2085                UseList[hist]['Pref.Ori.'][0] = 'MD'
     2086            else:
     2087                UseList[hist]['Pref.Ori.'][0] = 'SH'
     2088            UpdateDData()           
     2089
     2090        def OnPORef(event):
     2091            Obj = event.GetEventObject()
     2092            hist = Indx[Obj.GetId()]
     2093            UseList[hist]['Pref.Ori.'][2] = Obj.GetValue()
     2094           
     2095        def OnPOVal(event):
    20822096            Obj = event.GetEventObject()
    20832097            hist = Indx[Obj.GetId()]
     
    20852099                mdVal = float(Obj.GetValue())
    20862100                if mdVal > 0:
    2087                     UseList[hist]['MDtexture'][0] = mdVal
     2101                    UseList[hist]['Pref.Ori.'][1] = mdVal
    20882102            except ValueError:
    20892103                pass
    2090             Obj.SetValue("%.3f"%(UseList[hist]['MDtexture'][0]))          #reset in case of error
    2091            
    2092         def OnMDAxis(event):
     2104            Obj.SetValue("%.3f"%(UseList[hist]['Pref.Ori.'][1]))          #reset in case of error
     2105           
     2106        def OnPOAxis(event):
    20932107            Obj = event.GetEventObject()
    20942108            hist = Indx[Obj.GetId()]
     
    20972111                hkl = [int(Saxis[i]) for i in range(3)]
    20982112            except (ValueError,IndexError):
    2099                 hkl = UseList[hist]['MDtexture'][2]
     2113                hkl = UseList[hist]['Pref.Ori.'][3]
    21002114            if not np.any(np.array(hkl)):
    2101                 hkl = UseList[hist]['MDtexture'][2]
    2102             UseList[hist]['MDtexture'][2] = hkl
     2115                hkl = UseList[hist]['Pref.Ori.'][3]
     2116            UseList[hist]['Pref.Ori.'][3] = hkl
    21032117            h,k,l = hkl
    21042118            Obj.SetValue('%3d %3d %3d'%(h,k,l))
    21052119           
     2120        def OnPOOrder(event):
     2121            Obj = event.GetEventObject()
     2122            hist = Indx[Obj.GetId()]
     2123            Order = int(Obj.GetValue())
     2124            UseList[hist]['Pref.Ori.'][4] = Order
     2125            UseList[hist]['Pref.Ori.'][5] = SetPOCoef(Order,hist)
     2126            UpdateDData()
     2127
     2128        def SetPOCoef(Order,hist):
     2129            cofNames = G2lat.GenSHCoeff(SGData['SGLaue'],None,Order)     #cylindrical sample symmetry
     2130            newPOCoef = dict(zip(cofNames,np.zeros(len(cofNames))))
     2131            POCoeff = UseList[hist]['Pref.Ori.'][5]
     2132            for cofName in POCoeff:
     2133                if cofName in  cofNames:
     2134                    newPOCoef[cofName] = POCoeff[cofName]
     2135            return newPOCoef
     2136       
    21062137        def OnExtRef(event):
    21072138            Obj = event.GetEventObject()
     
    21442175        shOrder.Bind(wx.EVT_COMBOBOX,OnShOrder)
    21452176        shSizer.Add(shOrder,0,wx.ALIGN_CENTER_VERTICAL)
    2146         shSizer.Add((5,0),0)
    2147         shRef = wx.CheckBox(dataDisplay,label=' Refine texture?')
    2148         shRef.SetValue(textureData['SH Coeff'][0])
    2149         shRef.Bind(wx.EVT_CHECKBOX, OnSHRefine)
    2150         shSizer.Add(shRef,0,wx.ALIGN_CENTER_VERTICAL)
    2151         shShow = wx.CheckBox(dataDisplay,label=' Show coeff.?')
    2152         shShow.SetValue(textureData['SHShow'])
    2153         shShow.Bind(wx.EVT_CHECKBOX, OnSHShow)
    2154         shSizer.Add(shShow,0,wx.ALIGN_CENTER_VERTICAL)
    2155         mainSizer.Add(shSizer,0,0)
    2156         mainSizer.Add((0,5),0)
    2157         PTSizer = wx.FlexGridSizer(2,4,5,5)
    2158         PTSizer.Add(wx.StaticText(dataDisplay,-1,' Texture plot type: '),0,wx.ALIGN_CENTER_VERTICAL)
    2159         choices = ['Axial pole distribution','Pole figure','Inverse pole figure']           
    2160         pfType = wx.ComboBox(dataDisplay,-1,value=str(textureData['PlotType']),choices=choices,
    2161             style=wx.CB_READONLY|wx.CB_DROPDOWN)
    2162         pfType.Bind(wx.EVT_COMBOBOX,OnPfType)
    2163         PTSizer.Add(pfType,0,wx.ALIGN_CENTER_VERTICAL)
    2164         PTSizer.Add(wx.StaticText(dataDisplay,-1,' Projection type: '),0,wx.ALIGN_CENTER_VERTICAL)
    2165         projSel = wx.ComboBox(dataDisplay,-1,value=self.Projection,choices=['equal area','stereographic'],
    2166             style=wx.CB_READONLY|wx.CB_DROPDOWN)
    2167         projSel.Bind(wx.EVT_COMBOBOX,OnProjSel)
    2168         PTSizer.Add(projSel,0,wx.ALIGN_CENTER_VERTICAL)
    2169         if textureData['PlotType'] in ['Pole figure','Axial pole distribution']:
    2170             PTSizer.Add(wx.StaticText(dataDisplay,-1,' Pole figure HKL: '),0,wx.ALIGN_CENTER_VERTICAL)
    2171             PH = textureData['PFhkl']
    2172             pfVal = wx.TextCtrl(dataDisplay,-1,'%d,%d,%d'%(PH[0],PH[1],PH[2]),style=wx.TE_PROCESS_ENTER)
    2173         else:
    2174             PTSizer.Add(wx.StaticText(dataDisplay,-1,' Inverse pole figure XYZ: '),0,wx.ALIGN_CENTER_VERTICAL)
    2175             PX = textureData['PFxyz']
    2176             pfVal = wx.TextCtrl(dataDisplay,-1,'%3.1f,%3.1f,%3.1f'%(PX[0],PX[1],PX[2]),style=wx.TE_PROCESS_ENTER)
    2177         pfVal.Bind(wx.EVT_TEXT_ENTER,OnPFValue)
    2178         pfVal.Bind(wx.EVT_KILL_FOCUS,OnPFValue)
    2179         PTSizer.Add(pfVal,0,wx.ALIGN_CENTER_VERTICAL)
    2180         PTSizer.Add(wx.StaticText(dataDisplay,-1,' Color scheme'),0,wx.ALIGN_CENTER_VERTICAL)
    2181         choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
    2182         choice.sort()
    2183         colorSel = wx.ComboBox(dataDisplay,-1,value=self.ContourColor,choices=choice,
    2184             style=wx.CB_READONLY|wx.CB_DROPDOWN)
    2185         colorSel.Bind(wx.EVT_COMBOBOX,OnColorSel)
    2186         PTSizer.Add(colorSel,0,wx.ALIGN_CENTER_VERTICAL)       
    2187         mainSizer.Add(PTSizer,0,wx.ALIGN_CENTER_VERTICAL)
    2188        
    2189         mainSizer.Add((0,5),0)
    2190         if textureData['SHShow']:
    2191             mainSizer.Add(wx.StaticText(dataDisplay,-1,'Spherical harmonic coefficients: '),0,wx.ALIGN_CENTER_VERTICAL)
     2177        if textureData['Order']:
     2178            shSizer.Add((5,0),0)
     2179            shRef = wx.CheckBox(dataDisplay,label=' Refine texture?')
     2180            shRef.SetValue(textureData['SH Coeff'][0])
     2181            shRef.Bind(wx.EVT_CHECKBOX, OnSHRefine)
     2182            shSizer.Add(shRef,0,wx.ALIGN_CENTER_VERTICAL)
     2183            shShow = wx.CheckBox(dataDisplay,label=' Show coeff.?')
     2184            shShow.SetValue(textureData['SHShow'])
     2185            shShow.Bind(wx.EVT_CHECKBOX, OnSHShow)
     2186            shSizer.Add(shShow,0,wx.ALIGN_CENTER_VERTICAL)
     2187            mainSizer.Add(shSizer,0,0)
    21922188            mainSizer.Add((0,5),0)
    2193             ODFSizer = wx.FlexGridSizer(2,8,2,2)
    2194             ODFIndx = {}
    2195             ODFkeys = textureData['SH Coeff'][1].keys()
    2196             ODFkeys.sort()
    2197             for item in ODFkeys:
    2198                 ODFSizer.Add(wx.StaticText(dataDisplay,-1,item),0,wx.ALIGN_CENTER_VERTICAL)
    2199                 ODFval = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%8.3f'%(textureData['SH Coeff'][1][item]),style=wx.TE_PROCESS_ENTER)
    2200                 ODFIndx[ODFval.GetId()] = item
    2201                 ODFval.Bind(wx.EVT_TEXT_ENTER,OnODFValue)
    2202                 ODFval.Bind(wx.EVT_KILL_FOCUS,OnODFValue)
    2203                 ODFSizer.Add(ODFval,0,wx.ALIGN_CENTER_VERTICAL)
    2204             mainSizer.Add(ODFSizer,0,wx.ALIGN_CENTER_VERTICAL)
     2189            PTSizer = wx.FlexGridSizer(2,4,5,5)
     2190            PTSizer.Add(wx.StaticText(dataDisplay,-1,' Texture plot type: '),0,wx.ALIGN_CENTER_VERTICAL)
     2191            choices = ['Axial pole distribution','Pole figure','Inverse pole figure']           
     2192            pfType = wx.ComboBox(dataDisplay,-1,value=str(textureData['PlotType']),choices=choices,
     2193                style=wx.CB_READONLY|wx.CB_DROPDOWN)
     2194            pfType.Bind(wx.EVT_COMBOBOX,OnPfType)
     2195            PTSizer.Add(pfType,0,wx.ALIGN_CENTER_VERTICAL)
     2196            PTSizer.Add(wx.StaticText(dataDisplay,-1,' Projection type: '),0,wx.ALIGN_CENTER_VERTICAL)
     2197            projSel = wx.ComboBox(dataDisplay,-1,value=self.Projection,choices=['equal area','stereographic'],
     2198                style=wx.CB_READONLY|wx.CB_DROPDOWN)
     2199            projSel.Bind(wx.EVT_COMBOBOX,OnProjSel)
     2200            PTSizer.Add(projSel,0,wx.ALIGN_CENTER_VERTICAL)
     2201            if textureData['PlotType'] in ['Pole figure','Axial pole distribution']:
     2202                PTSizer.Add(wx.StaticText(dataDisplay,-1,' Pole figure HKL: '),0,wx.ALIGN_CENTER_VERTICAL)
     2203                PH = textureData['PFhkl']
     2204                pfVal = wx.TextCtrl(dataDisplay,-1,'%d,%d,%d'%(PH[0],PH[1],PH[2]),style=wx.TE_PROCESS_ENTER)
     2205            else:
     2206                PTSizer.Add(wx.StaticText(dataDisplay,-1,' Inverse pole figure XYZ: '),0,wx.ALIGN_CENTER_VERTICAL)
     2207                PX = textureData['PFxyz']
     2208                pfVal = wx.TextCtrl(dataDisplay,-1,'%3.1f,%3.1f,%3.1f'%(PX[0],PX[1],PX[2]),style=wx.TE_PROCESS_ENTER)
     2209            pfVal.Bind(wx.EVT_TEXT_ENTER,OnPFValue)
     2210            pfVal.Bind(wx.EVT_KILL_FOCUS,OnPFValue)
     2211            PTSizer.Add(pfVal,0,wx.ALIGN_CENTER_VERTICAL)
     2212            PTSizer.Add(wx.StaticText(dataDisplay,-1,' Color scheme'),0,wx.ALIGN_CENTER_VERTICAL)
     2213            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
     2214            choice.sort()
     2215            colorSel = wx.ComboBox(dataDisplay,-1,value=self.ContourColor,choices=choice,
     2216                style=wx.CB_READONLY|wx.CB_DROPDOWN)
     2217            colorSel.Bind(wx.EVT_COMBOBOX,OnColorSel)
     2218            PTSizer.Add(colorSel,0,wx.ALIGN_CENTER_VERTICAL)       
     2219            mainSizer.Add(PTSizer,0,wx.ALIGN_CENTER_VERTICAL)
     2220           
    22052221            mainSizer.Add((0,5),0)
    2206         mainSizer.Add((0,5),0)
    2207         mainSizer.Add(wx.StaticText(dataDisplay,-1,'Sample orientation angles: '),0,wx.ALIGN_CENTER_VERTICAL)
    2208         mainSizer.Add((0,5),0)
    2209         angSizer = wx.BoxSizer(wx.HORIZONTAL)
    2210         angIndx = {}
    2211         valIndx = {}
    2212         for item in ['Sample omega','Sample chi','Sample phi']:
    2213             angRef = wx.CheckBox(dataDisplay,label=item+': ')
    2214             angRef.SetValue(textureData[item][0])
    2215             angIndx[angRef.GetId()] = item
    2216             angRef.Bind(wx.EVT_CHECKBOX, OnAngRef)
    2217             angSizer.Add(angRef,0,wx.ALIGN_CENTER_VERTICAL)
    2218             angVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%8.2f'%(textureData[item][1]),style=wx.TE_PROCESS_ENTER)
    2219             valIndx[angVal.GetId()] = item
    2220             angVal.Bind(wx.EVT_TEXT_ENTER,OnAngValue)
    2221             angVal.Bind(wx.EVT_KILL_FOCUS,OnAngValue)
    2222             angSizer.Add(angVal,0,wx.ALIGN_CENTER_VERTICAL)
    2223             angSizer.Add((5,0),0)
    2224         mainSizer.Add(angSizer,0,wx.ALIGN_CENTER_VERTICAL)
     2222            if textureData['SHShow']:
     2223                mainSizer.Add(wx.StaticText(dataDisplay,-1,'Spherical harmonic coefficients: '),0,wx.ALIGN_CENTER_VERTICAL)
     2224                mainSizer.Add((0,5),0)
     2225                ODFSizer = wx.FlexGridSizer(2,8,2,2)
     2226                ODFIndx = {}
     2227                ODFkeys = textureData['SH Coeff'][1].keys()
     2228                ODFkeys.sort()
     2229                for item in ODFkeys:
     2230                    ODFSizer.Add(wx.StaticText(dataDisplay,-1,item),0,wx.ALIGN_CENTER_VERTICAL)
     2231                    ODFval = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%8.3f'%(textureData['SH Coeff'][1][item]),style=wx.TE_PROCESS_ENTER)
     2232                    ODFIndx[ODFval.GetId()] = item
     2233                    ODFval.Bind(wx.EVT_TEXT_ENTER,OnODFValue)
     2234                    ODFval.Bind(wx.EVT_KILL_FOCUS,OnODFValue)
     2235                    ODFSizer.Add(ODFval,0,wx.ALIGN_CENTER_VERTICAL)
     2236                mainSizer.Add(ODFSizer,0,wx.ALIGN_CENTER_VERTICAL)
     2237                mainSizer.Add((0,5),0)
     2238            mainSizer.Add((0,5),0)
     2239            mainSizer.Add(wx.StaticText(dataDisplay,-1,'Sample orientation angles: '),0,wx.ALIGN_CENTER_VERTICAL)
     2240            mainSizer.Add((0,5),0)
     2241            angSizer = wx.BoxSizer(wx.HORIZONTAL)
     2242            angIndx = {}
     2243            valIndx = {}
     2244            for item in ['Sample omega','Sample chi','Sample phi']:
     2245                angRef = wx.CheckBox(dataDisplay,label=item+': ')
     2246                angRef.SetValue(textureData[item][0])
     2247                angIndx[angRef.GetId()] = item
     2248                angRef.Bind(wx.EVT_CHECKBOX, OnAngRef)
     2249                angSizer.Add(angRef,0,wx.ALIGN_CENTER_VERTICAL)
     2250                angVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%8.2f'%(textureData[item][1]),style=wx.TE_PROCESS_ENTER)
     2251                valIndx[angVal.GetId()] = item
     2252                angVal.Bind(wx.EVT_TEXT_ENTER,OnAngValue)
     2253                angVal.Bind(wx.EVT_KILL_FOCUS,OnAngValue)
     2254                angSizer.Add(angVal,0,wx.ALIGN_CENTER_VERTICAL)
     2255                angSizer.Add((5,0),0)
     2256            mainSizer.Add(angSizer,0,wx.ALIGN_CENTER_VERTICAL)
     2257        else:  #finish the texture output when order = 0
     2258            mainSizer.Add(shSizer,0,0)
     2259            mainSizer.Add((0,5),0)           
    22252260        mainSizer.Add(wx.StaticText(dataDisplay,-1,'Histogram data for '+PhaseName+':'),0,wx.ALIGN_CENTER_VERTICAL)
    22262261        for item in keyList:
     
    22352270            if UseList[item]['Show']:
    22362271                scaleSizer = wx.BoxSizer(wx.HORIZONTAL)
    2237                 scaleRef = wx.CheckBox(dataDisplay,label=' Scale factor: ')
     2272                scaleRef = wx.CheckBox(dataDisplay,label=' Phase fraction: ')
    22382273                scaleRef.SetValue(UseList[item]['Scale'][1])
    22392274                Indx[scaleRef.GetId()] = item
     
    22522287                mainSizer.Add((0,5),0)
    22532288                sizeSizer = wx.BoxSizer(wx.HORIZONTAL)
    2254                 choices = ['isotropic','uniaxial',]
     2289                sizeSizer.Add(wx.StaticText(dataDisplay,-1,' Size model: '),0,wx.ALIGN_CENTER_VERTICAL)
     2290                choices = ['isotropic','uniaxial','ellipsoidal']
    22552291                sizeType = wx.ComboBox(dataDisplay,wx.ID_ANY,value=UseList[item]['Size'][0],choices=choices,
    22562292                    style=wx.CB_READONLY|wx.CB_DROPDOWN)
     
    23002336                    sizeSizer.Add((5,0),0)                   
    23012337                    mainSizer.Add(sizeSizer)
     2338                elif UseList[item]['Size'][0] == 'ellipsoidal':
     2339                    mainSizer.Add(sizeSizer)
     2340                    mainSizer.Add((0,5),0)
     2341                    sizeSizer = wx.BoxSizer(wx.HORIZONTAL)
     2342                    parms = zip(['S11','S22','S33','S12','S13','S23'],UseList[item]['Size'][4],
     2343                        UseList[item]['Size'][5],range(6))
     2344                    for Pa,val,ref,id in parms:
     2345                        sizeRef = wx.CheckBox(dataDisplay,label=Pa)
     2346                        sizeRef.SetValue(ref)
     2347                        Indx[sizeRef.GetId()] = [item,id]
     2348                        sizeRef.Bind(wx.EVT_CHECKBOX, OnSizeRef)
     2349                        sizeSizer.Add(sizeRef,0,wx.ALIGN_CENTER_VERTICAL)
     2350                        sizeVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%.1f'%(val),style=wx.TE_PROCESS_ENTER)
     2351                        Indx[sizeVal.GetId()] = [item,id]
     2352                        sizeVal.Bind(wx.EVT_TEXT_ENTER,OnSizeVal)
     2353                        sizeVal.Bind(wx.EVT_KILL_FOCUS,OnSizeVal)
     2354                        sizeSizer.Add(sizeVal,0,wx.ALIGN_CENTER_VERTICAL)
     2355                        sizeSizer.Add((5,0),0)
     2356                    sizeSizer.Add((5,0),0)                   
     2357                    mainSizer.Add(sizeSizer)
    23022358               
    23032359                strainSizer = wx.BoxSizer(wx.HORIZONTAL)
     2360                strainSizer.Add(wx.StaticText(dataDisplay,-1,' Mustrain model: '),0,wx.ALIGN_CENTER_VERTICAL)
    23042361                choices = ['isotropic','uniaxial','generalized',]
    23052362                strainType = wx.ComboBox(dataDisplay,wx.ID_ANY,value=UseList[item]['Mustrain'][0],choices=choices,
     
    23732430                        strainSizer.Add(strainVal,0,wx.ALIGN_CENTER_VERTICAL)
    23742431                    mainSizer.Add(strainSizer)
    2375                 #MD texture  'MDtexture':[1.0,False,[0,0,1]]
    2376                 mdSizer = wx.BoxSizer(wx.HORIZONTAL)
    2377                 mdRef = wx.CheckBox(dataDisplay,label=' March-Dollase texture ratio: ')
    2378                 mdRef.SetValue(UseList[item]['MDtexture'][1])
    2379                 Indx[mdRef.GetId()] = item
    2380                 mdRef.Bind(wx.EVT_CHECKBOX, OnMDRef)
    2381                 mdSizer.Add(mdRef,0,wx.ALIGN_CENTER_VERTICAL)
    2382                 mdVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,
    2383                     '%.3f'%(UseList[item]['MDtexture'][0]),style=wx.TE_PROCESS_ENTER)
    2384                 Indx[mdVal.GetId()] = item
    2385                 mdVal.Bind(wx.EVT_TEXT_ENTER,OnMDVal)
    2386                 mdVal.Bind(wx.EVT_KILL_FOCUS,OnMDVal)
    2387                 mdSizer.Add(mdVal,0,wx.ALIGN_CENTER_VERTICAL)
    2388                 mdSizer.Add(wx.StaticText(dataDisplay,-1,' Unique axis, H K L: '),0,wx.ALIGN_CENTER_VERTICAL)
    2389                 h,k,l = UseList[item]['MDtexture'][2]
    2390                 mdAxis = wx.TextCtrl(dataDisplay,-1,'%3d %3d %3d'%(h,k,l),style=wx.TE_PROCESS_ENTER)
    2391                 Indx[mdAxis.GetId()] = item
    2392                 mdAxis.Bind(wx.EVT_TEXT_ENTER,OnMDAxis)
    2393                 mdAxis.Bind(wx.EVT_KILL_FOCUS,OnMDAxis)
    2394                 mdSizer.Add(mdAxis,0,wx.ALIGN_CENTER_VERTICAL)
    2395                 mainSizer.Add(mdSizer)
    2396                 mainSizer.Add((0,5),0)
    2397                
     2432                #texture  'Pref. Ori.':['MD',1.0,False,[0,0,1],0,[]] last two for 'SH' are SHorder & coeff
     2433#                poSizer = wx.BoxSizer(wx.HORIZONTAL)
     2434                poSizer = wx.FlexGridSizer(1,6,5,5)
     2435                POData = UseList[item]['Pref.Ori.']
     2436                choice = ['March-Dollase','Spherical harmonics']
     2437                POtype = choice[['MD','SH'].index(POData[0])]
     2438                poSizer.Add(wx.StaticText(dataDisplay,-1,' Preferred orientation model '),0,wx.ALIGN_CENTER_VERTICAL)
     2439                POType = wx.ComboBox(dataDisplay,wx.ID_ANY,value=POtype,choices=choice,
     2440                    style=wx.CB_READONLY|wx.CB_DROPDOWN)
     2441                Indx[POType.GetId()] = item
     2442                POType.Bind(wx.EVT_COMBOBOX, OnPOType)
     2443                poSizer.Add(POType)
     2444                if POData[0] == 'MD':
     2445                    mainSizer.Add(poSizer)
     2446                    poSizer = wx.BoxSizer(wx.HORIZONTAL)
     2447                    poRef = wx.CheckBox(dataDisplay,label=' March-Dollase ratio: ')
     2448                    poRef.SetValue(POData[2])
     2449                    Indx[poRef.GetId()] = item
     2450                    poRef.Bind(wx.EVT_CHECKBOX,OnPORef)
     2451                    poSizer.Add(poRef,0,wx.ALIGN_CENTER_VERTICAL)
     2452                    poVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,
     2453                        '%.3f'%(POData[1]),style=wx.TE_PROCESS_ENTER)
     2454                    Indx[poVal.GetId()] = item
     2455                    poVal.Bind(wx.EVT_TEXT_ENTER,OnPOVal)
     2456                    poVal.Bind(wx.EVT_KILL_FOCUS,OnPOVal)
     2457                    poSizer.Add(poVal,0,wx.ALIGN_CENTER_VERTICAL)
     2458                    poSizer.Add(wx.StaticText(dataDisplay,-1,' Unique axis, H K L: '),0,wx.ALIGN_CENTER_VERTICAL)
     2459                    h,k,l =POData[3]
     2460                    poAxis = wx.TextCtrl(dataDisplay,-1,'%3d %3d %3d'%(h,k,l),style=wx.TE_PROCESS_ENTER)
     2461                    Indx[poAxis.GetId()] = item
     2462                    poAxis.Bind(wx.EVT_TEXT_ENTER,OnPOAxis)
     2463                    poAxis.Bind(wx.EVT_KILL_FOCUS,OnPOAxis)
     2464                    poSizer.Add(poAxis,0,wx.ALIGN_CENTER_VERTICAL)
     2465                    mainSizer.Add(poSizer)
     2466                else:           #'SH'
     2467                    poSizer.Add(wx.StaticText(dataDisplay,-1,' Harmonic order: '),0,wx.ALIGN_CENTER_VERTICAL)
     2468                    poOrder = wx.ComboBox(dataDisplay,wx.ID_ANY,value=str(POData[4]),choices=[str(2*i) for i in range(18)],
     2469                        style=wx.CB_READONLY|wx.CB_DROPDOWN)
     2470                    Indx[poOrder.GetId()] = item
     2471                    poOrder.Bind(wx.EVT_COMBOBOX,OnPOOrder)
     2472                    poSizer.Add(poOrder,0,wx.ALIGN_CENTER_VERTICAL)
     2473                    poRef = wx.CheckBox(dataDisplay,label=' Refine? ')
     2474                    poRef.SetValue(POData[2])
     2475                    Indx[poRef.GetId()] = item
     2476                    poRef.Bind(wx.EVT_CHECKBOX,OnPORef)
     2477                    poSizer.Add(poRef,0,wx.ALIGN_CENTER_VERTICAL)
     2478                    mainSizer.Add(poSizer)
     2479                    if POData[4]:
     2480                        mainSizer.Add(wx.StaticText(dataDisplay,-1,'Spherical harmonic coefficients: '),0,wx.ALIGN_CENTER_VERTICAL)
     2481                        mainSizer.Add((0,5),0)
     2482                        ODFSizer = wx.FlexGridSizer(2,8,2,2)
     2483                        ODFIndx = {}
     2484                        ODFkeys = POData[5].keys()
     2485                        ODFkeys.sort()
     2486                        for odf in ODFkeys:
     2487                            ODFSizer.Add(wx.StaticText(dataDisplay,-1,odf),0,wx.ALIGN_CENTER_VERTICAL)
     2488                            ODFval = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%8.3f'%(POData[5][odf]),style=wx.TE_PROCESS_ENTER)
     2489                            ODFIndx[ODFval.GetId()] = odf
     2490    #                        ODFval.Bind(wx.EVT_TEXT_ENTER,OnODFValue)
     2491    #                        ODFval.Bind(wx.EVT_KILL_FOCUS,OnODFValue)
     2492                            ODFSizer.Add(ODFval,0,wx.ALIGN_CENTER_VERTICAL)
     2493                        mainSizer.Add(ODFSizer,0,wx.ALIGN_CENTER_VERTICAL)
     2494                        mainSizer.Add((0,5),0)
     2495                mainSizer.Add((0,5),0)               
    23982496                #Extinction  'Extinction':[0.0,False]
    23992497                extSizer = wx.BoxSizer(wx.HORIZONTAL)
     
    24692567                    for i in result:
    24702568                        histoName = TextList[i]
     2569                        pId = G2gd.GetPatternTreeItemId(self,self.root,histoName)
    24712570                        UseList[histoName] = {'Histogram':histoName,'Show':False,
    2472                             'Scale':[1.0,False],'MDtexture':[1.0,False,[0,0,1]],
    2473                             'Size':['isotropic',[10000.,0,],[False,False],[0,0,1]],
    2474                             'Mustrain':['isotropic',[1.0,0.0],[False,False],[0,0,1],
     2571                            'Scale':[1.0,False],'Pref.Ori.':['MD',1.0,False,[0,0,1],0,{}],
     2572                            'Size':['isotropic',[10.,10.,],[False,False],[0,0,1],6*[0.0,],6*[False,]],
     2573                            'Mustrain':['isotropic',[1.0,1.0],[False,False],[0,0,1],
    24752574                                NShkl*[0.01,],NShkl*[False,]],                           
    24762575                            'Extinction':[0.0,False]}
     2576                        refList = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,pId,'Reflection Lists'))
     2577                        refList[generalData['Name']] = []                       
    24772578                    data['Histograms'] = UseList
    24782579                    UpdateDData()
     
    25012602
    25022603    def FillPawleyReflectionsGrid():
    2503         if data['Histograms']:
    2504             self.dataFrame.PawleyMenu.FindItemById(G2gd.wxID_PAWLEYLOAD).Enable(True)
    2505             self.dataFrame.PawleyMenu.FindItemById(G2gd.wxID_PAWLEYIMPORT).Enable(True)
    2506         else:
    2507             self.dataFrame.PawleyMenu.FindItemById(G2gd.wxID_PAWLEYLOAD).Enable(False)
    2508             self.dataFrame.PawleyMenu.FindItemById(G2gd.wxID_PAWLEYIMPORT).Enable(False)
    25092604                       
    25102605        def KeyEditPawleyGrid(event):
    2511             colList = PawleyRefl.GetSelectedCols()
     2606            colList = self.PawleyRefl.GetSelectedCols()
    25122607            PawleyPeaks = data['Pawley ref']
    25132608            if event.GetKeyCode() == wx.WXK_RETURN:
     
    25182613                event.Skip(True)
    25192614            elif colList:
    2520                 PawleyRefl.ClearSelection()
     2615                self.PawleyRefl.ClearSelection()
    25212616                key = event.GetKeyCode()
    25222617                for col in colList:
     
    25322627            rowLabels = []
    25332628            for i in range(len(PawleyPeaks)): rowLabels.append(str(i+1))
    2534             colLabels = ['h','k','l','mul','2-theta','sigma','refine','Iobs','Icalc']
    2535             Types = [wg.GRID_VALUE_LONG,wg.GRID_VALUE_LONG,wg.GRID_VALUE_LONG,wg.GRID_VALUE_LONG,
    2536                 wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
    2537                 wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,5']
     2629            colLabels = ['h','k','l','mul','d','refine','Fsq(hkl)','sig(Fsq)']
     2630            Types = 4*[wg.GRID_VALUE_LONG,]+[wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,]+ \
     2631                2*[wg.GRID_VALUE_FLOAT+':10,2',]
    25382632            PawleyTable = G2gd.Table(PawleyPeaks,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    2539             PawleyRefl.SetTable(PawleyTable, True)
    2540             PawleyRefl.Bind(wx.EVT_KEY_DOWN, KeyEditPawleyGrid)                 
    2541             PawleyRefl.SetMargins(0,0)
    2542             PawleyRefl.AutoSizeColumns(False)
     2633            self.PawleyRefl.SetTable(PawleyTable, True)
     2634            self.PawleyRefl.Bind(wx.EVT_KEY_DOWN, KeyEditPawleyGrid)                 
     2635            self.PawleyRefl.SetMargins(0,0)
     2636            self.PawleyRefl.AutoSizeColumns(False)
    25432637            self.dataFrame.setSizePosLeft([500,300])
    25442638                   
    25452639    def OnPawleyLoad(event):
    2546         sig = lambda Th,U,V,W: 1.17741*math.sqrt(U*tand(Th)**2+V*tand(Th)+W)/100.
    2547         gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))/100.
    2548         gamFW = lambda s,g: math.sqrt(s**2+(0.4654996*g)**2)+.5345004*g  #Ubaldo Bafile - private communication
    2549         choice = data['Histograms'].keys()
    2550         dlg = wx.SingleChoiceDialog(self,'Select','Powder histogram',choice)
    2551         if dlg.ShowModal() == wx.ID_OK:
    2552             histogram = choice[dlg.GetSelection()]
    2553             Id  = G2gd.GetPatternTreeItemId(self, self.root, histogram)
    2554             Iparms = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Instrument Parameters'))
    2555             try:                                            #try single wavelength
    2556                 lam = Iparms[1][Iparms[3].index('Lam')]
    2557             except ValueError:                              #Ka1 & Ka2 present
    2558                 lam = Iparms[1][Iparms[3].index('Lam1')]
    2559             GU = Iparms[1][Iparms[3].index('U')]               
    2560             GV = Iparms[1][Iparms[3].index('V')]               
    2561             GW = Iparms[1][Iparms[3].index('W')]               
    2562             LX = Iparms[1][Iparms[3].index('X')]               
    2563             LY = Iparms[1][Iparms[3].index('Y')]
    2564         dlg.Destroy()
    25652640        generalData = data['General']
    25662641        dmin = generalData['Pawley dmin']
     
    25682643        A = G2lat.cell2A(cell)
    25692644        SGData = generalData['SGData']
    2570         Laue = SGData['SGLaue']
    2571         SGLatt = SGData['SGLatt']
    2572         SGUniq = SGData['SGUniq']       
    2573         HKLd = G2lat.GenHLaue(dmin,Laue,SGLatt,SGUniq,A)
     2645        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
    25742646        PawleyPeaks = []
    25752647        wx.BeginBusyCursor()
    25762648        try:
    25772649            for h,k,l,d in HKLd:
    2578                 ext,mul = G2spc.GenHKL([h,k,l],SGData)[:2]
     2650                ext,mul = G2spc.GenHKLf([h,k,l],SGData)[:2]
    25792651                if not ext:
    2580                     th = asind(lam/(2.0*d))
    2581                     H = gamFW(sig(th,GU,GV,GW),gam(th,LX,LY))/2.35482
    2582                     PawleyPeaks.append([h,k,l,mul,2*th,H,False,0,0])
     2652                    PawleyPeaks.append([h,k,l,mul,d,False,100.0,1.0])
    25832653        finally:
    25842654            wx.EndBusyCursor()
    25852655        data['Pawley ref'] = PawleyPeaks
    25862656        FillPawleyReflectionsGrid()
    2587                
    2588            
    2589     def OnPawleyImport(event):
    2590         dlg = wx.FileDialog(self, 'Choose file with Pawley reflections', '.', '',
    2591             'GSAS Pawley files (*.RFL)|*.RFL',wx.OPEN)
    2592         if self.dirname:
    2593             dlg.SetDirectory(self.dirname)
    2594         try:
    2595             if dlg.ShowModal() == wx.ID_OK:
    2596                 PawleyFile = dlg.GetPath()
    2597                 self.dirname = dlg.GetDirectory()
    2598                 choice = data['Histograms'].keys()
    2599                 dlg2 = wx.SingleChoiceDialog(self,'Select','Powder histogram',choice)
    2600                 if dlg2.ShowModal() == wx.ID_OK:
    2601                     histogram = choice[dlg2.GetSelection()]
    2602                     Id  = G2gd.GetPatternTreeItemId(self, self.root, histogram)
    2603                     Iparms = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Instrument Parameters'))
    2604                 dlg2.Destroy()
    2605                
    2606                 PawleyPeaks = G2IO.GetPawleyPeaks(PawleyFile)
    2607                 data['Pawley ref'] = PawleyPeaks
    2608                 FillPawleyReflectionsGrid()
    2609         finally:
    2610             dlg.Destroy()
    2611            
     2657                           
    26122658    def OnPawleyDelete(event):
    26132659        dlg = wx.MessageDialog(self,'Do you really want to delete Pawley reflections?','Delete',
     
    26682714            self.dataFrame.SetMenuBar(self.dataFrame.PawleyMenu)
    26692715            self.dataFrame.Bind(wx.EVT_MENU, OnPawleyLoad, id=G2gd.wxID_PAWLEYLOAD)
    2670             self.dataFrame.Bind(wx.EVT_MENU, OnPawleyImport, id=G2gd.wxID_PAWLEYIMPORT)
    26712716            self.dataFrame.Bind(wx.EVT_MENU, OnPawleyDelete, id=G2gd.wxID_PAWLEYDELETE)           
    26722717            FillPawleyReflectionsGrid()           
     
    26842729        DData = wx.ScrolledWindow(self.dataDisplay)
    26852730        self.dataDisplay.AddPage(DData,'Data')
    2686         PawleyRefl = G2gd.GSGrid(self.dataDisplay)
    2687         self.dataDisplay.AddPage(PawleyRefl,'Pawley reflections')
     2731        self.PawleyRefl = G2gd.GSGrid(self.dataDisplay)
     2732        self.dataDisplay.AddPage(self.PawleyRefl,'Pawley reflections')
    26882733    else:
    26892734        DData = wx.ScrolledWindow(self.dataDisplay)
  • trunk/GSASIIplot.py

    r335 r342  
    456456                if self.itemPicked:
    457457                    Page.canvas.SetToolTipString('%9.3f'%(xpos))
    458                 if self.PickId and self.PatternTree.GetItemText(self.PickId) in ['Index Peak List','Unit Cells List']:
     458                if self.PickId:
    459459                    found = []
    460                     if len(HKL):
    461                         view = Page.toolbar._views.forward()[0][:2]
    462                         wid = view[1]-view[0]
    463                         found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
    464                     if len(found):
    465                         h,k,l = found[0][:3]
    466                         Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
    467                     else:
    468                         Page.canvas.SetToolTipString('')
     460                    if self.PatternTree.GetItemText(self.PickId) in ['Index Peak List','Unit Cells List','Reflection Lists']:
     461                        if len(HKL):
     462                            view = Page.toolbar._views.forward()[0][:2]
     463                            wid = view[1]-view[0]
     464                            found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
     465                        if len(found):
     466                            h,k,l = found[0][:3]
     467                            Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
     468                        else:
     469                            Page.canvas.SetToolTipString('')
    469470
    470471            except TypeError:
     
    573574    Ymax = 1.0
    574575    lenX = 0
    575     HKL = np.array(self.HKL)
     576    if self.PatternTree.GetItemText(PickId) in ['Reflection Lists']:
     577        Phases = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId,'Reflection Lists'))
     578        HKL = []
     579        for peak in Phases[self.RefList]:
     580            HKL.append(peak[:6])
     581        HKL = np.array(HKL)
     582    else:
     583        HKL = np.array(self.HKL)
    576584    for Pattern in PlotList:
    577585        xye = Pattern[1]
     
    659667                else:
    660668                    Plot.plot(X,Y,colors[N%6],picker=False)
    661     if PickId and self.PatternTree.GetItemText(PickId) in ['Index Peak List','Unit Cells List']:
     669    if PickId:
    662670        Values,Names = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Instrument Parameters'))[1::2]
    663671        Parms = dict(zip(Names,Values))
     
    666674        except KeyError:
    667675            wave = Parms['Lam1']
    668         peaks = np.array((self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Index Peak List'))))
    669         for peak in peaks:
    670             if self.qPlot:
    671                 Plot.axvline(4*np.pi*sind(peak[0]/2.0)/wave,color='b')
    672             else:
    673                 Plot.axvline(peak[0],color='b')
    674         for hkl in self.HKL:
    675             if self.qPlot:
    676                 Plot.axvline(4*np.pi*sind(hkl[5]/2.0)/wave,color='r',dashes=(5,5))
    677             else:
    678                 Plot.axvline(hkl[5],color='r',dashes=(5,5))
     676        if self.PatternTree.GetItemText(PickId) in ['Index Peak List','Unit Cells List']:
     677            peaks = np.array((self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Index Peak List'))))
     678            for peak in peaks:
     679                if self.qPlot:
     680                    Plot.axvline(4*np.pi*sind(peak[0]/2.0)/wave,color='b')
     681                else:
     682                    Plot.axvline(peak[0],color='b')
     683            for hkl in self.HKL:
     684                if self.qPlot:
     685                    Plot.axvline(4*np.pi*sind(hkl[5]/2.0)/wave,color='r',dashes=(5,5))
     686                else:
     687                    Plot.axvline(hkl[5],color='r',dashes=(5,5))
     688        elif self.PatternTree.GetItemText(PickId) in ['Reflection Lists']:
     689            Phases = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId,'Reflection Lists'))
     690            for pId,phase in enumerate(Phases):
     691                pId += 1
     692                peaks = Phases[phase]
     693                for peak in peaks:
     694                    if self.qPlot:
     695                        Plot.plot(2*np.pi/peak[4],offset*N-pId*Ymax*.005,colors[pId%6]+'|',mew=1,ms=8,picker=2.)
     696                    else:
     697                        Plot.plot(peak[5],offset*N-pId*Ymax*.005,colors[pId%6]+'|',mew=1,ms=8,picker=2.)
     698           
    679699    if self.Contour:
    680700        acolor = mpl.cm.get_cmap(self.ContourColor)
     
    11571177            elif muStrain[0] == 'uniaxial':
    11581178                def uniaxMustrain(xyz,muiso,muaniso,axes):
    1159                     cp = abs(np.dot(xyz,axes))
    1160                     S = muiso+muaniso*cp
    1161                     return S*xyz*math.pi/0.018      #centidegrees to radians!
     1179                    Z = np.array(axes)
     1180                    cp = abs(np.dot(xyz,Z))
     1181                    sp = np.sqrt(1.-cp**2)
     1182                    R = muiso*muaniso/np.sqrt((muiso*cp)**2+(muaniso*sp)**2)
     1183#                    S = muiso+muaniso*cp           #old GSAS - wrong math!!
     1184                    return R*xyz*math.pi/0.018      #centidegrees to radians!
    11621185                muiso,muaniso = muStrain[1][:2]
    11631186                axes = np.inner(A,np.array(muStrain[3]))
     
    11761199                    uvw = np.inner(A.T,xyz)
    11771200                    Strm = np.array(G2spc.MustrainCoeff(uvw,SGData))
    1178                     sum = np.sqrt(np.sum(np.multiply(Shkl,Strm)))*math.pi/0.018      #centidegrees to radians!
     1201                    sum = np.sum(np.multiply(Shkl,Strm))
     1202                    sum = np.where(sum > 0.01,sum,0.01)
     1203                    sum = np.sqrt(sum)*math.pi/0.018      #centidegrees to radians!
    11791204                    return sum*xyz
    11801205                Shkl = np.array(muStrain[4])
     
    11891214                   
    11901215            if np.any(X) and np.any(Y) and np.any(Z):
    1191                 Plot.plot_surface(X,Y,Z,rstride=1,cstride=1,color='g')
     1216                Plot.plot_surface(X,Y,Z,rstride=1,cstride=1,color='g',linewidth=1)
    11921217               
    11931218            Plot.set_xlabel('X')
  • trunk/GSASIIpwd.py

    r335 r342  
    9292    return plist
    9393
    94 def ValuesOut(parmdict, varylist):
    95     '''Use before call to leastsq to setup list of values for the parameters
    96     in parmdict, as selected by key in varylist'''
    97     return [parmdict[key] for key in varylist]
    98    
    99 def ValuesIn(parmdict, varylist, values):
    100     ''' Use after call to leastsq to update the parameter dictionary with
    101     values corresponding to keys in varylist'''
    102     parmdict.update(zip(varylist,values))
    103    
    10494def Transmission(Geometry,Abs,Diam):
    10595#Calculate sample transmission
     
    482472fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
    483473               
    484 def getBackground(parmDict,bakType,xdata):
     474def getBackground(pfx,parmDict,bakType,xdata):
     475    yb = np.zeros_like(xdata)
    485476    if bakType == 'chebyschev':
    486         yb = np.zeros_like(xdata)
    487477        iBak = 0
    488478        while True:
    489             key = 'Back'+str(iBak)
     479            key = pfx+'Back:'+str(iBak)
    490480            try:
    491481                yb += parmDict[key]*(xdata-xdata[0])**iBak
    492482                iBak += 1
    493483            except KeyError:
    494                 return yb
    495                
     484                break
     485    return yb
     486
     487def calcPeakFFT(x,fxns,widths,pos,args):
     488    dx = x[1]-x[0]
     489    z = np.empty([len(fxns),len(x)])
     490    for i,(fxn,width,arg) in enumerate(zip(fxns,widths,args)):
     491        z[i] = fxn.pdf(x,loc=pos,scale=width,*arg)*dx
     492    Z = fft(z)
     493    if len(fxns)%2:
     494        return ifft(Z.prod(axis=0)).real
     495    else:
     496        return fftshift(ifft(Z.prod(axis=0))).real
     497                   
     498def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
     499   
     500    DX = xdata[1]-xdata[0]
     501    widths,fmin,fmax = getWidths(pos,sig,gam,shl)
     502    x = np.linspace(pos-fmin,pos+fmin,1024)
     503    if pos > 90:
     504        fmin,fmax = [fmax,fmin]         
     505    dx = x[1]-x[0]
     506    arg = [pos,shl/10.,dx,]
     507    Df = fcjde.pdf(x,*arg,loc=pos,scale=1.0)
     508    if len(np.nonzero(Df)[0])>5:
     509        fxns = [st.norm,st.cauchy,fcjde]
     510        args = [[],[],arg]
     511    else:
     512        fxns = [st.norm,st.cauchy]
     513        args = [[],[]]
     514    Df = calcPeakFFT(x,fxns,widths,pos,args)
     515    Df /= np.sum(Df)
     516    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
     517    return intens*Df(xdata)*DX/dx       #*10 to get close to old fxn???
     518                                   
    496519def getPeakProfile(parmDict,xdata,varyList,bakType):
    497520   
    498     def calcPeakFFT(x,fxns,widths,pos,args):
    499         dx = x[1]-x[0]
    500         z = np.empty([len(fxns),len(x)])
    501         for i,(fxn,width,arg) in enumerate(zip(fxns,widths,args)):
    502             z[i] = fxn.pdf(x,loc=pos,scale=width,*arg)*dx
    503         Z = fft(z)
    504         if len(fxns)%2:
    505             return ifft(Z.prod(axis=0)).real
    506         else:
    507             return fftshift(ifft(Z.prod(axis=0))).real
    508                        
    509     def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
    510        
    511         DX = xdata[1]-xdata[0]
    512         widths,fmin,fmax = getWidths(pos,sig,gam,shl)
    513         x = np.linspace(pos-fmin,pos+fmin,1024)
    514         if pos > 90:
    515             fmin,fmax = [fmax,fmin]         
    516         dx = x[1]-x[0]
    517         arg = [pos,shl/10.,dx,]
    518         Df = fcjde.pdf(x,*arg,loc=pos,scale=1.0)
    519         if len(np.nonzero(Df)[0])>5:
    520             fxns = [st.norm,st.cauchy,fcjde]
    521             args = [[],[],arg]
    522         else:
    523             fxns = [st.norm,st.cauchy]
    524             args = [[],[]]
    525         Df = calcPeakFFT(x,fxns,widths,pos,args)
    526         Df /= np.sum(Df)
    527         Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
    528         return intens*Df(xdata)*DX/dx       #*10 to get close to old fxn???
    529                    
    530     yb = getBackground(parmDict,bakType,xdata)
     521    yb = getBackground('',parmDict,bakType,xdata)
    531522    yc = np.zeros_like(yb)
    532523    U = parmDict['U']
     
    591582    return widths,fmin,fmax
    592583               
    593 def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data):
     584def Dict2Values(parmdict, varylist):
     585    '''Use before call to leastsq to setup list of values for the parameters
     586    in parmdict, as selected by key in varylist'''
     587    return [parmdict[key] for key in varylist]
     588   
     589def Values2Dict(parmdict, varylist, values):
     590    ''' Use after call to leastsq to update the parameter dictionary with
     591    values corresponding to keys in varylist'''
     592    parmdict.update(zip(varylist,values))
     593   
     594def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False):
    594595   
    595596    def SetBackgroundParms(Background):
    596597        bakType,bakFlag = Background[:2]
    597598        backVals = Background[3:]
    598         backNames = ['Back'+str(i) for i in range(len(backVals))]
     599        backNames = ['Back:'+str(i) for i in range(len(backVals))]
    599600        if bakFlag: #returns backNames as varyList = backNames
    600601            return bakType,dict(zip(backNames,backVals)),backNames
     
    606607        while True:
    607608            try:
    608                 bakName = 'Back'+str(iBak)
     609                bakName = 'Back:'+str(iBak)
    609610                Background[iBak+3] = parmList[bakName]
    610611                iBak += 1
     
    620621            for i,back in enumerate(Background[3:]):
    621622                ptstr += ptfmt % (back)
    622                 sigstr += ptfmt % (sigDict['Back'+str(i)])
     623                sigstr += ptfmt % (sigDict['Back:'+str(i)])
    623624            print ptstr
    624625            print sigstr
     
    675676        print sigstr
    676677
    677     def setPeaksParms(Peaks):
     678    def SetPeaksParms(Peaks):
    678679        peakNames = []
    679680        peakVary = []
     
    734735        return M
    735736       
     737    Ftol = 0.01
     738    if oneCycle:
     739        Ftol = 1.0
    736740    x,y,w,yc,yb,yd = data               #these are numpy arrays!
    737741    xBeg = np.searchsorted(x,Limits[0])
     
    739743    bakType,bakDict,bakVary = SetBackgroundParms(Background)
    740744    dataType,insDict,insVary = SetInstParms(Inst)
    741     peakDict,peakVary = setPeaksParms(Peaks)
     745    peakDict,peakVary = SetPeaksParms(Peaks)
    742746    parmDict = {}
    743747    parmDict.update(bakDict)
     
    747751    while True:
    748752        begin = time.time()
    749         values =  np.array(ValuesOut(parmDict, varyList))
     753        values =  np.array(Dict2Values(parmDict, varyList))
    750754        if FitPgm == 'LSQ':
    751755            dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0,
     
    755759            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
    756760            try:
    757                 result = so.leastsq(errPeakProfile,values,
    758                     args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg),full_output=True)
     761                result = so.leastsq(errPeakProfile,values,full_output=True,ftol=Ftol,
     762                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
    759763            finally:
    760764                dlg.Destroy()
     
    762766            chisq = np.sum(result[2]['fvec']**2)
    763767            ncyc = int(result[2]['nfev']/len(varyList))
    764             ValuesIn(parmDict, varyList, result[0])
     768            Values2Dict(parmDict, varyList, result[0])
    765769            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
    766770            GOF = chisq/(xFin-xBeg-len(varyList))
     
    786790       
    787791    sigDict = dict(zip(varyList,sig))
    788     yb[xBeg:xFin] = getBackground(parmDict,bakType,x[xBeg:xFin])
     792    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
    789793    yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType)
    790794    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
     
    913917    msg = 'test '
    914918    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
    915     gplot.plot(xdata,getBackground(parmDict0,bakType,xdata))   
     919    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
    916920    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
    917921    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
    918     fplot.plot(xdata,getBackground(parmDict1,bakType,xdata))   
     922    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
    919923    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
     924   
     925def test1():
     926    if NeedTestData: TestData()
    920927    time0 = time.time()
    921928    for i in range(100):
    922929        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
    923     print '100+6*Ka1-2 peaks=3600 peaks',time.time()-time0
     930    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
     931   
    924932   
    925933
  • trunk/GSASIIpwdGUI.py

    r335 r342  
    8181        OnPeakFit('LSQ')
    8282       
     83    def OnOneCycle(event):
     84        OnPeakFit('LSQ',oneCycle=True)
     85       
    8386    def OnBGFSPeakFit(event):
    8487        OnPeakFit('BFGS')
    8588               
    86     def OnPeakFit(FitPgm):
     89    def OnPeakFit(FitPgm,oneCycle=False):
    8790        SaveState()
    8891        print 'Peak Fitting:'
     
    99102        wx.BeginBusyCursor()
    100103        try:
    101             G2pwd.DoPeakFit(FitPgm,peaks,background,limits,inst,data)
     104            G2pwd.DoPeakFit(FitPgm,peaks,background,limits,inst,data,oneCycle)
    102105        finally:
    103106            wx.EndBusyCursor()   
     
    237240    self.Bind(wx.EVT_MENU, OnUnDo, id=G2gd.wxID_UNDO)
    238241    self.Bind(wx.EVT_MENU, OnLSQPeakFit, id=G2gd.wxID_LSQPEAKFIT)
     242    self.Bind(wx.EVT_MENU, OnOneCycle, id=G2gd.wxID_LSQONECYCLE)
    239243#    self.Bind(wx.EVT_MENU, OnBGFSPeakFit, id=G2gd.wxID_BFGSPEAKFIT)
    240244    self.Bind(wx.EVT_MENU, OnResetSigGam, id=G2gd.wxID_RESETSIGGAM)
     
    242246    if data:
    243247        self.dataFrame.PeakFit.Enable(True)
     248        self.dataFrame.PFOneCycle.Enable(True)
    244249    self.PickTable = []
    245250    rowLabels = []
     
    346351    colLabels = ['Tmin','Tmax']
    347352    rowLabels = ['original','changed']
    348     Types = [wg.GRID_VALUE_FLOAT+':10,3',wg.GRID_VALUE_FLOAT+':10,3']
     353    Types = 2*[wg.GRID_VALUE_FLOAT+':10,3',]
    349354    self.LimitsTable = G2gd.Table(data,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    350355    self.dataFrame.SetLabel('Limits')
     
    598603        Status = self.dataFrame.CreateStatusBar()   
    599604    self.dataDisplay = wx.Panel(self.dataFrame)
    600    
     605
     606#patch
     607    if not 'Gonio. radius' in data:
     608        data['Gonio. radius'] = 200.0
     609#patch end
     610   
     611    parms = [['Gonio. radius',' Goniometer radius(mm): ','%.2f',]]
    601612    if data['Type'] == 'Debye-Scherrer':
    602         parms = [['DisplaceX',' Sample X displacement: ','%.4f',],
    603             ['DisplaceY',' Sample Y displacement: ','%.4f',],
    604             ['Absorption',' Sample absorption: ','%.4f',],]
     613        parms += [['DisplaceX',' Sample X displacement(\xb5m): ','%.2f',],
     614            ['DisplaceY',' Sample Y displacement(\xb5m): ','%.2f',],
     615            ['Absorption',' Sample absorption(\xb5r): ','%.4f',],]
    605616    elif data['Type'] == 'Bragg-Brentano':
    606         parms = [['Shift',' Sample displacement: ','%.4f',],
    607             ['Transparency',' Sample transparency: ','%.4f'],]
    608     parms.append(['Temperature',' Sample temperature: ','%.2f'])
    609     parms.append(['Pressure',' Sample pressure: ','%.3f'])
    610     parms.append(['Humidity',' Sample humidity: ','%.1f'])
    611     parms.append(['Voltage',' Sample voltage: ','%.3f'])
    612     parms.append(['Force',' Applied load: ','%.3f'])
     617        parms += [['Shift',' Sample displacement(\xb5m): ','%.2f',],
     618            ['Transparency',' Sample transparency(1/\xb5eff,\xc5): ','%.4f'],]
     619    parms.append(['Temperature',' Sample temperature(K): ','%.2f'])
     620    parms.append(['Pressure',' Sample pressure(MPa): ','%.3f'])
     621    parms.append(['Humidity',' Sample humidity(%): ','%.1f'])
     622    parms.append(['Voltage',' Sample voltage(V): ','%.3f'])
     623    parms.append(['Force',' Applied load(MN): ','%.3f'])
    613624    objList = {}
    614625
     
    657668               
    658669    mainSizer = wx.BoxSizer(wx.VERTICAL)
    659     mainSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' Sample parameters: '),0,wx.ALIGN_CENTER_VERTICAL)
     670    mainSizer.Add(wx.StaticText(self.dataDisplay,label=' Sample parameters: '),0,wx.ALIGN_CENTER_VERTICAL)
    660671    mainSizer.Add((5,5),0)
    661     scaleSizer = wx.BoxSizer(wx.HORIZONTAL)
     672    parmSizer = wx.FlexGridSizer(9,2,5,0)
    662673    scaleRef = wx.CheckBox(self.dataDisplay,label=' Histogram scale factor: ')
    663674    scaleRef.SetValue(data['Scale'][1])
    664675    scaleRef.Bind(wx.EVT_CHECKBOX, OnScaleRef)
    665     scaleSizer.Add(scaleRef,0,wx.ALIGN_CENTER_VERTICAL)
     676    parmSizer.Add(scaleRef,0,wx.ALIGN_CENTER_VERTICAL)
    666677    scaleVal = wx.TextCtrl(self.dataDisplay,wx.ID_ANY,
    667678        '%.4f'%(data['Scale'][0]),style=wx.TE_PROCESS_ENTER)
    668679    scaleVal.Bind(wx.EVT_TEXT_ENTER,OnScaleVal)
    669680    scaleVal.Bind(wx.EVT_KILL_FOCUS,OnScaleVal)
    670     scaleSizer.Add(scaleVal,0,wx.ALIGN_CENTER_VERTICAL)
    671     mainSizer.Add(scaleSizer)
    672     mainSizer.Add((0,5),0)
     681    parmSizer.Add(scaleVal,0,wx.ALIGN_CENTER_VERTICAL)
    673682    typeSizer = wx.BoxSizer(wx.HORIZONTAL)
    674683    choices = ['Debye-Scherrer','Bragg-Brentano',]
     
    676685        style=wx.CB_READONLY|wx.CB_DROPDOWN)
    677686    histoType.Bind(wx.EVT_COMBOBOX, OnHistoType)
    678     typeSizer.Add(histoType)
    679     mainSizer.Add(typeSizer)
    680     mainSizer.Add((0,5),0)
     687    parmSizer.Add(histoType)
     688    parmSizer.Add((5,5),0)
    681689   
    682690    for parm in parms:
    683         parmSizer = wx.BoxSizer(wx.HORIZONTAL)
    684691        if 'list' in str(type(data[parm[0]])):
    685692            parmRef = wx.CheckBox(self.dataDisplay,label=parm[1])
     
    698705        parmVal.Bind(wx.EVT_TEXT_ENTER,OnParmVal)
    699706        parmVal.Bind(wx.EVT_KILL_FOCUS,OnParmVal)
    700         parmSizer.Add(parmVal,0,wx.ALIGN_CENTER_VERTICAL)
    701         mainSizer.Add(parmSizer)
    702         mainSizer.Add((0,5),0)   
     707        parmSizer.Add(parmVal,1,wx.EXPAND)
     708    mainSizer.Add(parmSizer)
     709    mainSizer.Add((0,5),0)   
    703710   
    704711    mainSizer.Layout()   
     
    927934            G2plt.PlotPatterns(self)
    928935       
    929        
     936    def CopyUnitCell(event):
     937        controls,bravais,cells,dmin = self.PatternTree.GetItemPyData(UnitCellsId)
     938        for Cell in cells:
     939            if Cell[-1]:
     940                break
     941        cell = Cell[2:9]
     942        controls[4] = 1
     943        controls[5] = bravaisSymb[cell[0]]
     944        controls[6:12] = cell[1:8]
     945        controls[12] = G2lat.calc_V(G2lat.cell2A(controls[6:12]))
     946        self.PatternTree.SetItemPyData(UnitCellsId,[controls,bravais,cells,dmin])
     947        UpdateUnitCellsGrid(self,data)
     948       
     949        self.dataFrame.RefineCell.Enable(True)
     950               
    930951    def RefineCell(event):
    931952        def cellPrint(ibrav,A):
     
    955976        cell = controls[6:12]
    956977        A = G2lat.cell2A(cell)
    957         print controls[5]
    958978        ibrav = bravaisSymb.index(controls[5])
    959979        dmin = G2indx.getDmin(peaks)-0.005
     980        self.HKL = G2lat.GenHBravais(dmin,ibrav,A)
     981        G2indx.IndexPeaks(peaks,self.HKL)
    960982        if controls[0]:
    961983            Lhkl,M20,X20,Aref,Zero = G2indx.refinePeaksZ(peaks,inst[1],ibrav,A,controls[1])           
     
    10171039            UpdateUnitCellsGrid(self,data)
    10181040               
    1019     def CopyUnitCell(event):
    1020         controls,bravais,cells,dmin = self.PatternTree.GetItemPyData(UnitCellsId)
    1021         for Cell in cells:
    1022             if Cell[-1]:
    1023                 break
    1024         cell = Cell[2:9]
    1025         controls[4] = 1
    1026         controls[5] = bravaisSymb[cell[0]]
    1027         controls[6:12] = cell[1:8]
    1028         controls[12] = G2lat.calc_V(G2lat.cell2A(controls[6:12]))
    1029         self.PatternTree.SetItemPyData(UnitCellsId,[controls,bravais,cells,dmin])
    1030         UpdateUnitCellsGrid(self,data)
    1031         self.dataFrame.RefineCell.Enable(True)
    1032        
    1033     def MakeNewPhase(event):
    1034         if not G2gd.GetPatternTreeItemId(self,self.root,'Phases'):
    1035             sub = self.PatternTree.AppendItem(parent=self.root,text='Phases')
    1036         else:
    1037             sub = G2gd.GetPatternTreeItemId(self,self.root,'Phases')
    1038         PhaseName = ''
    1039         dlg = wx.TextEntryDialog(None,'Enter a name for this phase','Phase Name Entry','New phase',
    1040             style=wx.OK)
    1041         if dlg.ShowModal() == wx.ID_OK:
    1042             PhaseName = dlg.GetValue()
    1043         dlg.Destroy()
    1044         cells = self.PatternTree.GetItemPyData(UnitCellsId)[2]
    1045         for Cell in cells:
    1046             if Cell[-1]:
    1047                 break
    1048         cell = Cell[2:10]       
    1049         sub = self.PatternTree.AppendItem(parent=sub,text=PhaseName)
    1050         E,SGData = G2spc.SpcGroup(spaceGroups[cell[0]])
    1051         self.PatternTree.SetItemPyData(sub, \
    1052             {'General':{'Name':'phase name','Type':'nuclear','SGData':SGData,
    1053             'Cell':[False,]+cell[1:],
    1054             'Pawley dmin':0.25},'Atoms':[],'Drawing':{},'Histograms':{}})
    1055         Status.SetStatusText('Change space group if needed')
    1056            
    10571041    def RefreshUnitCellsGrid(event):
    10581042        cells,dmin = self.PatternTree.GetItemPyData(UnitCellsId)[2:]
     
    10761060                    G2plt.PlotPatterns(self)
    10771061       
     1062    def MakeNewPhase(event):
     1063        if not G2gd.GetPatternTreeItemId(self,self.root,'Phases'):
     1064            sub = self.PatternTree.AppendItem(parent=self.root,text='Phases')
     1065        else:
     1066            sub = G2gd.GetPatternTreeItemId(self,self.root,'Phases')
     1067        PhaseName = ''
     1068        dlg = wx.TextEntryDialog(None,'Enter a name for this phase','Phase Name Entry','New phase',
     1069            style=wx.OK)
     1070        try:
     1071            if dlg.ShowModal() == wx.ID_OK:
     1072                PhaseName = dlg.GetValue()
     1073                cells = self.PatternTree.GetItemPyData(UnitCellsId)[2]
     1074                for Cell in cells:
     1075                    if Cell[-1]:
     1076                        break
     1077                cell = Cell[2:10]       
     1078                sub = self.PatternTree.AppendItem(parent=sub,text=PhaseName)
     1079                E,SGData = G2spc.SpcGroup(spaceGroups[cell[0]])
     1080                self.PatternTree.SetItemPyData(sub, \
     1081                    {'General':{'Name':PhaseName,'Type':'nuclear','SGData':SGData,
     1082                    'Cell':[False,]+cell[1:],
     1083                    'Pawley dmin':1.00},'Atoms':[],'Drawing':{},'Histograms':{}})
     1084                Status.SetStatusText('Change space group if needed')
     1085        finally:
     1086            dlg.Destroy()
     1087           
    10781088    if self.dataDisplay:
    10791089        self.dataFrame.Clear()
     
    11701180    zeroVar.Bind(wx.EVT_CHECKBOX,OnZeroVar)
    11711181    littleSizer.Add(zeroVar,0,wx.ALIGN_CENTER_VERTICAL)
    1172     hklShow = wx.CheckBox(self.dataDisplay,label="  Show starting hkl positions")
     1182    hklShow = wx.CheckBox(self.dataDisplay,label="  Show hkl positions")
    11731183    hklShow.Bind(wx.EVT_CHECKBOX,OnHklShow)
    11741184    littleSizer.Add(hklShow,0,wx.ALIGN_CENTER_VERTICAL)
     
    12171227        rowLabels = []
    12181228        colLabels = ['M20','X20','use','Bravais','a','b','c','alpha','beta','gamma','Volume']
    1219         Types = [wg.GRID_VALUE_FLOAT+':10,2',wg.GRID_VALUE_NUMBER,wg.GRID_VALUE_BOOL,wg.GRID_VALUE_STRING,
    1220             wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_FLOAT+':10,5',
    1221             wg.GRID_VALUE_FLOAT+':10,3',wg.GRID_VALUE_FLOAT+':10,3',wg.GRID_VALUE_FLOAT+':10,3',
    1222             wg.GRID_VALUE_FLOAT+':10,2']
     1229        Types = [wg.GRID_VALUE_FLOAT+':10,2',wg.GRID_VALUE_NUMBER,wg.GRID_VALUE_BOOL,wg.GRID_VALUE_STRING,]+ \
     1230            3*[wg.GRID_VALUE_FLOAT+':10,5',]+3*[wg.GRID_VALUE_FLOAT+':10,3',]+ \
     1231            [wg.GRID_VALUE_FLOAT+':10,2']
    12231232        numRows = len(cells)
    12241233        table = []
     
    12371246        gridDisplay.SetTable(UnitCellsTable, True)
    12381247        self.dataFrame.CopyCell.Enable(True)
    1239 #        gridDisplay.Bind(wg.EVT_GRID_CELL_CHANGE, RefreshUnitCellsGrid)
    12401248        gridDisplay.Bind(wg.EVT_GRID_CELL_LEFT_CLICK,RefreshUnitCellsGrid)
    12411249        gridDisplay.SetMargins(0,0)
     
    12491257                    gridDisplay.SetReadOnly(r,c,isReadOnly=True)
    12501258        gridDisplay.SetSize(bottomSize)
     1259
     1260def UpdateReflectionGrid(self,data):
     1261    if not data:
     1262        print 'No phases, no reflections'
     1263        return
     1264    phases = data.keys()
     1265   
     1266    def OnSelectPhase(event):
     1267        dlg = wx.SingleChoiceDialog(self,'Select','Phase',phases)
     1268        try:
     1269            if dlg.ShowModal() == wx.ID_OK:
     1270                sel = dlg.GetSelection()
     1271                self.RefList = phases[sel]
     1272                UpdateReflectionGrid(self,data)
     1273        finally:
     1274            dlg.Destroy()
     1275        G2plt.PlotPatterns(self)
     1276       
     1277       
     1278    if self.dataDisplay:
     1279        self.dataFrame.Clear()
     1280    self.dataFrame.SetMenuBar(self.dataFrame.ReflMenu)
     1281    if not self.dataFrame.GetStatusBar():
     1282        Status = self.dataFrame.CreateStatusBar()   
     1283    self.dataDisplay = wx.Panel(self.dataFrame)
     1284    self.Bind(wx.EVT_MENU, OnSelectPhase, id=G2gd.wxID_SELECTPHASE)
     1285    self.dataFrame.SelectPhase.Enable(False)
     1286    if len(data) > 1:
     1287        self.dataFrame.SelectPhase.Enable(True)
     1288    rowLabels = []
     1289    refList = []
     1290    for h,k,l,m,d,pos,sig,gam,fo,fc,x in data[self.RefList]:
     1291        refList.append([h,k,l,m,d,pos,sig,gam,fo,fc])
     1292    for i in range(len(refList)): rowLabels.append(str(i+1))
     1293    colLabels = ['H','K','L','mul','d','pos','sig','gam','Fosq','Fcsq',]
     1294    Types = 4*[wg.GRID_VALUE_LONG,]+6*[wg.GRID_VALUE_FLOAT+':10,4',]
     1295    self.PeakTable = G2gd.Table(refList,rowLabels=rowLabels,colLabels=colLabels,types=Types)
     1296    self.dataFrame.SetLabel('Reflection List for '+self.RefList)
     1297    self.dataDisplay = G2gd.GSGrid(parent=self.dataFrame)
     1298    self.dataDisplay.SetTable(self.PeakTable, True)
     1299    self.dataDisplay.SetMargins(0,0)
     1300    self.dataDisplay.AutoSizeColumns(False)
     1301    self.dataFrame.setSizePosLeft([535,350])
    12511302
    12521303def UpdatePDFGrid(self,data):
  • trunk/GSASIIspc.py

    r230 r342  
    299299    else:
    300300        return zip(XYZEquiv,Idup,Cell)
     301
     302def GenHKLf(HKL,SGData,Friedel=False):
     303    '''
     304    Uses old GSAS Fortran routine genhkl.for
     305    input:
     306        HKL - [h,k,l]
     307        SGData - space group data obtained from SpcGroup
     308        Friedel = True to retain Friedel pairs for centrosymmetric case
     309    returns:
     310        iabsnt = True is reflection is forbidden by symmetry
     311        mulp = reflection multiplicity including Fridel pairs
     312        Uniq = numpy array of equivalent hkl in descending order of h,k,l
     313    '''
     314    hklf = HKL+[0,]
     315    Ops = SGData['SGOps']
     316    OpM = np.array([op[0] for op in Ops])
     317    OpT = np.array([op[1] for op in Ops])
     318    Inv = SGData['SGInv']
     319    Cen = np.array([cen for cen in SGData['SGCen']])
     320   
     321    Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
     322    h,k,l,f = Uniq
     323    Uniq=np.array(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq]))
     324    Uniq = np.array(Uniq)
     325   
     326    return iabsnt,2*mulp,Uniq       #include Friedel pairs in powder mulp
     327
    301328       
    302 def GenHKL(HKL,SGData,Friedel=False):
     329def GenHKL(HKL,SGData,Friedel=False):          #gives wrong answers!
    303330    '''
    304331    Generate set of equivalent reflections
     
    312339        Uniq = numpy array of equivalent hkl in descending order of h,k,l
    313340        Phs = numpy array of corresponding phase factors (multiples of 2pi)
     341       
     342    NB: only works on one HKL at a time -
     343        it can be made to do an array of HKLs - not any faster!
    314344    '''
    315345    iabsnt = False
     
    322352    # get operators & expand if centrosymmetric
    323353    Ops = SGData['SGOps']
    324     opM = np.array([op[0].T for op in Ops])
     354    opM = np.array([op[0] for op in Ops])
    325355    opT = np.array([op[1] for op in Ops])
    326356    if Inv:
     
    330360    eqH = np.inner(opM,H)
    331361    HT = np.dot(opT,H)%1.
    332     dup = len(ma.nonzero([ma.allclose(HKL,hkl,atol=0.001) for hkl in eqH])[0])
     362    dup = len(ma.nonzero([ma.allclose(H,hkl,atol=0.001) for hkl in eqH])[0])
    333363    mulp = len(eqH)/dup
    334364    # generate unique reflection set (with or without Friedel pairs)
     
    347377    HPT = np.around(HPT-HPT[0],4)%1.
    348378    iabsnt |= np.any(HPT)
    349     return iabsnt,mulp,Uniq,Phs
     379    return iabsnt,mulp,Uniq
     380                                   
     381#def GenHKLs(H,SGData,Friedel=False):
     382#    '''
     383#    Generate set of equivalent reflections
     384#    input:
     385#        H - np.array of [h,k,l] assumed to exclude lattice centering extinctions
     386#        SGData - space group data obtained from SpcGroup
     387#        Friedel = True to retain Friedel pairs for centrosymmetric case
     388#    returns:
     389#        iabsnt = True is reflection is forbidden by symmetry
     390#        mulp = reflection multiplicity including Fridel pairs
     391#        Uniq = numpy array of equivalent hkl in descending order of h,k,l
     392#        Phs = numpy array of corresponding phase factors (multiples of 2pi)
     393#    this is not any faster than GenHKL above - oh well
     394#    '''
     395#    H = np.array(H,dtype=float)
     396#    Inv = SGData['SGInv']
     397#    Ops = SGData['SGOps']
     398#    opM = np.array([op[0] for op in Ops])
     399#    opT = np.array([op[1] for op in Ops])
     400#    if Inv:
     401#        opM = np.concatenate((opM,-opM))
     402#        opT = np.concatenate((opT,-opT))
     403#    # generate full hkl table and phase factors; find duplicates for multiplicity
     404#    eqH = np.swapaxes(np.swapaxes(np.inner(opM,H),0,2),1,2)
     405#    HeqH = zip(H,eqH)
     406#    dup = np.array([sum([ma.allclose(h,hkl,atol=0.001) for hkl in HKL]) for h,HKL in HeqH])
     407#    mulp = len(eqH[0])/dup
     408#    # generate unique reflection set (with or without Friedel pairs)
     409#    HT = np.dot(H,opT.T)%1.
     410#    HPT = zip(mulp,dup,[[tuple(hp) for hp in HP] for HP in np.dstack((eqH,HT))])
     411#    Uniq = []
     412#    Phs = []
     413#    iabsnt = []
     414#    for mp,dp,hpt in HPT:
     415#        hpt.sort()
     416#        upt = hpt[::dp]
     417#        upt.reverse()
     418#        if Inv and not Friedel:
     419#            upt = upt[:len(upt)/2]           
     420#        uniq = np.split(upt,[3,4],axis=1)
     421#        Phs.append(uniq[1].T[0])
     422#        Uniq.append(uniq[0])
     423#        # determine if extinct
     424#        hpt = np.split(hpt,[3,4],axis=1)[1].T[0]
     425#        hpt = np.reshape(hpt,(mp,-1)).T
     426#        hpt = np.around(hpt-hpt[0],4)%1.
     427#        iabsnt.append(np.any(hpt))
     428#    return iabsnt,mulp,Uniq,Phs
    350429                                   
    351430def GetOprPtrName(key):           
  • trunk/GSASIIstruct.py

    r230 r342  
    88########### SVN repository information ###################
    99import sys
     10import os
    1011import os.path as ospath
    1112import numpy as np
     
    1314import time
    1415import math
     16import wx
    1517import GSASIIpath
    1618import GSASIIElem as G2el
    1719import GSASIIlattice as G2lat
    1820import GSASIIspc as G2spc
     21import GSASIIpwd as G2pwd
     22import GSASIImapvars as G2mv
     23import scipy.optimize as so
     24
     25sind = lambda x: np.sin(x*np.pi/180.)
     26cosd = lambda x: np.cos(x*np.pi/180.)
     27tand = lambda x: np.tan(x*np.pi/180.)
     28asind = lambda x: 180.*np.arcsin(x)/np.pi
     29
    1930
    2031def ShowBanner():
     
    3445    '''
    3546    file = open(GPXfile,'rb')
    36     HistogramNames = []
    3747    while True:
    3848        try:
     
    8292    return PhaseNames
    8393
    84 def GetPhaseData(GPXfile,PhaseName):
    85     ''' Returns the 'General' and 'Atoms' objects for PhaseName from GSASII gpx file
     94def GetAllPhaseData(GPXfile,PhaseName):
     95    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
    8696    input:
    8797        GPXfile = gpx full file name
    8898        PhaseName = phase name
    8999    return:
    90         General = dictionary of general phase info.
    91         Atoms = list of atom parameters
    92         these are returned empty if PhaseName not found
     100        phase dictionary
    93101    '''       
    94102    file = open(GPXfile,'rb')
     
    104112            for datus in data[1:]:
    105113                if datus[0] == PhaseName:
    106                     General = datus[1]['General']
    107                     Atoms = datus[1]['Atoms']                                           
     114                    break
    108115    file.close()
    109     return {'General':General,'Atoms':Atoms}
    110    
    111 def ShowPhaseData(phaseNames,PhaseData):
    112     print ' Phases:'
    113     for name in phaseNames:
    114         General = PhaseData[name]['General']
    115         Atoms = PhaseData[name]['Atoms']
    116         print '\n Phase name: ',General['Name']
    117         SGtext = G2spc.SGPrint(General['SGData'])
    118         for line in SGtext: print line
    119         cell = General['Cell']
    120         print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
    121             ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =','%.3f'%(cell[6]),' volume =','%.3f'%(cell[7])
    122         print ' Refine?',cell[0]
    123         print '\n Atoms:'
    124         line = '   name    type  refine?   x         y         z    '+ \
    125             '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
    126         if General['Type'] == 'magnetic':
    127             line += '   Mx     My     Mz'
    128         elif General['Type'] == 'macromolecular':
    129             line = ' res no  residue  chain '+line
    130         print line
    131         if General['Type'] == 'nuclear':
    132             print 135*'-'
    133             for at in Atoms:
    134                 line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
    135                     '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
    136                 if at[9] == 'I':
    137                     line += '%8.4f'%(at[10])+48*' '
    138                 else:
    139                     line += 8*' '
    140                     for i in range(6):
    141                         line += '%8.4f'%(at[11+i])
    142                 print line
    143 #        elif General['Type'] == 'magnetic':
    144 #        elif General['Type'] == 'macromolecular':
    145        
    146 
     116    return datus[1]
     117   
    147118def GetHistogramNames(GPXfile):
    148119    ''' Returns a list of histogram names found in GSASII gpx file
     
    165136    return HistogramNames
    166137   
     138def GetUsedHistogramsAndPhases(GPXfile):
     139    ''' Returns all histograms that are found in any phase
     140    and any phase that uses a histogram
     141    input:
     142        GPXfile = .gpx full file name
     143    return:
     144        Histograms = dictionary of histograms as {name:data,...}
     145        Phases = dictionary of phases that use histograms
     146    '''
     147    phaseNames = GetPhaseNames(GPXfile)
     148    phaseData = {}
     149    for name in phaseNames:
     150        phaseData[name] =  GetAllPhaseData(GPXfile,name)
     151    Histograms = {}
     152    Phases = {}
     153    pId = 0
     154    hId = 0
     155    for phase in phaseData:
     156        Phase = phaseData[phase]
     157        if Phase['Histograms']:
     158            if phase not in Phases:
     159                Phase['pId'] = pId
     160                pId += 1
     161                Phases[phase] = Phase
     162            for hist in Phase['Histograms']:
     163                if hist not in Histograms:
     164                    if 'PWDR' in hist[:4]:
     165                        Histograms[hist] = GetPWDRdata(GPXfile,hist)
     166                    elif 'HKLF' in hist[:4]:
     167                        Histograms[hist] = GetHKLFdata(GPXfile,hist)
     168                    #future restraint, etc. histograms here           
     169                    Histograms[hist]['hId'] = hId
     170                    hId += 1
     171    return Histograms,Phases
     172   
     173def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases):
     174    ''' Updates gpxfile from all histograms that are found in any phase
     175    and any phase that used a histogram
     176    input:
     177        GPXfile = .gpx full file name
     178        Histograms = dictionary of histograms as {name:data,...}
     179        Phases = dictionary of phases that use histograms
     180    '''
     181                       
     182    def GPXBackup(GPXfile):
     183        import distutils.file_util as dfu
     184        GPXpath,GPXname = ospath.split(GPXfile)
     185        Name = ospath.splitext(GPXname)[0]
     186        files = os.listdir(GPXpath)
     187        last = 0
     188        for name in files:
     189            name = name.split('.')
     190            if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
     191                last = max(last,int(name[1].strip('bak'))+1)
     192        GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
     193        dfu.copy_file(GPXfile,GPXback)
     194        return GPXback
     195       
     196    GPXback = GPXBackup(GPXfile)
     197    print '\n',130*'-'
     198    print 'Read from file:',GPXback
     199    print 'Save to file  :',GPXfile
     200    infile = open(GPXback,'rb')
     201    outfile = open(GPXfile,'wb')
     202    while True:
     203        try:
     204            data = cPickle.load(infile)
     205        except EOFError:
     206            break
     207        datum = data[0]
     208        print 'read: ',datum[0]
     209        if datum[0] == 'Phases':
     210            for iphase in range(len(data)):
     211                if data[iphase][0] in Phases:
     212                    phaseName = data[iphase][0]
     213                    data[iphase][1] = Phases[phaseName]
     214        try:
     215            histogram = Histograms[datum[0]]
     216            print 'found ',datum[0]
     217            data[0][1][1] = histogram['Data']
     218            for datus in data[1:]:
     219                print '    read: ',datus[0]
     220                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
     221                    datus[1] = histogram[datus[0]]
     222        except KeyError:
     223            pass
     224                               
     225        cPickle.dump(data,outfile,1)
     226    infile.close()
     227    outfile.close()
     228    print 'refinement save successful'
     229                   
    167230def GetPWDRdata(GPXfile,PWDRname):
    168231    ''' Returns powder data from GSASII gpx file
     
    171234        PWDRname = powder histogram name as obtained from GetHistogramNames
    172235    return:
    173         PWDRdata = powder data list:
     236        PWDRdata = powder data dictionary with:
     237            Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
    174238       
    175239    '''
    176240    file = open(GPXfile,'rb')
    177     PWDRdata = []
     241    PWDRdata = {}
    178242    while True:
    179243        try:
     
    183247        datum = data[0]
    184248        if datum[0] == PWDRname:
    185             PWDRdata = datum[1:][0]
    186                                                
    187            
     249            PWDRdata['Data'] = datum[1][1]          #powder data arrays
     250            PWDRdata[data[2][0]] = data[2][1]       #Limits
     251            PWDRdata[data[3][0]] = data[3][1]       #Background
     252            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
     253            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
     254            try:
     255                PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
     256            except IndexError:
     257                PWDRdata['Reflection lists'] = {}
    188258    file.close()
    189     return PWDRdata#,Limits,InstrumentParms
     259    return PWDRdata
    190260   
    191261def GetHKLFdata(GPXfile,HKLFname):
     
    226296                FFtable[El] = item
    227297    return FFtable
    228        
     298   
     299def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
     300    if SGLaue in ['-1','2/m','mmm']:
     301        return                      #no Pawley symmetry required constraints
     302    for varyI in pawleyVary:
     303        refI = int(varyI.split(':')[-1])
     304        ih,ik,il = PawleyRef[refI][:3]
     305        for varyJ in pawleyVary[:refI]:
     306            refJ = int(varyJ.split(':')[-1])
     307            jh,jk,jl = PawleyRef[refJ][:3]
     308            if SGLaue in ['4/m','4/mmm']:
     309                isum = ih**2+ik**2
     310                jsum = jh**2+jk**2
     311                if abs(il) == abs(jl) and isum == jsum:
     312                    G2mv.StoreEquivalence(varyJ,(varyI,))
     313            elif SGLaue in ['3R','3mR']:
     314                isum = ih**2+ik**2+il**2
     315                jsum = jh**2+jk**2*jl**2
     316                isum2 = ih*ik+ih*il+ik*il
     317                jsum2 = jh*jk+jh*jl+jk*jl
     318                if isum == jsum and isum2 == jsum2:
     319                    G2mv.StoreEquivalence(varyJ,(varyI,))
     320            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
     321                isum = ih**2+ik**2+ih*ik
     322                jsum = jh**2+jk**2+jh*jk
     323                if abs(il) == abs(jl) and isum == jsum:
     324                    G2mv.StoreEquivalence(varyJ,(varyI,))
     325            elif SGLaue in ['m3','m3m']:
     326                isum = ih**2+ik**2+il**2
     327                jsum = jh**2+jk**2+jl**2
     328                if isum == jsum:
     329                    G2mv.StoreEquivalence(varyJ,(varyI,))
     330       
     331def GetPhaseData(PhaseData):
     332   
     333    def cellVary(pfx,SGData):
     334        if SGData['SGLaue'] in ['-1',]:
     335            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
     336        elif SGData['SGLaue'] in ['2/m',]:
     337            if SGData['SGUniq'] == 'a':
     338                return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
     339            elif SGData['SGUniq'] == 'b':
     340                return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
     341            else:
     342                return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
     343        elif SGData['SGLaue'] in ['mmm',]:
     344            return [pfx+'A0',pfx+'A1',pfx+'A2']
     345        elif SGData['SGLaue'] in ['4/m','4/mmm']:
     346            return [pfx+'A0',pfx+'A2']
     347        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
     348            return [pfx+'A0',pfx+'A2']
     349        elif SGData['SGLaue'] in ['3R', '3mR']:
     350            return [pfx+'A0',pfx+'A3']                       
     351        elif SGData['SGLaue'] in ['m3m','m3']:
     352            return [pfx+'A0']
     353       
     354       
     355    print ' Phases:'
     356    phaseVary = []
     357    phaseDict = {}
     358    phaseConstr = {}
     359    pawleyLookup = {}
     360    for name in PhaseData:
     361        General = PhaseData[name]['General']
     362        pId = PhaseData[name]['pId']
     363        pfx = str(pId)+'::'
     364        Atoms = PhaseData[name]['Atoms']
     365        try:
     366            PawleyRef = PhaseData[name]['Pawley ref']
     367        except KeyError:
     368            PawleyRef = []
     369        print '\n Phase name: ',General['Name']
     370        SGData = General['SGData']
     371        SGtext = G2spc.SGPrint(SGData)
     372        for line in SGtext: print line
     373        cell = General['Cell']
     374        A = G2lat.cell2A(cell[1:7])
     375        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]})
     376        if cell[0]:
     377            phaseVary = cellVary(pfx,SGData)
     378        print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
     379            ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
     380            '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
     381        if Atoms:
     382            print '\n Atoms:'
     383            line = '   name    type  refine?   x         y         z    '+ \
     384                '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
     385            if General['Type'] == 'magnetic':
     386                line += '   Mx     My     Mz'
     387            elif General['Type'] == 'macromolecular':
     388                line = ' res no  residue  chain '+line
     389            print line
     390            if General['Type'] == 'nuclear':
     391                print 135*'-'
     392                for at in Atoms:
     393                    line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
     394                        '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
     395                    if at[9] == 'I':
     396                        line += '%8.4f'%(at[10])+48*' '
     397                    else:
     398                        line += 8*' '
     399                        for i in range(6):
     400                            line += '%8.4f'%(at[11+i])
     401                    print line
     402                    if 'X' in at[2]:
     403                        xId,xCoef = G2spc.GetCSxinel(at[7])
     404                    if 'U' in at[2]:
     405                        uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
     406#        elif General['Type'] == 'magnetic':
     407#        elif General['Type'] == 'macromolecular':
     408#       PWDR: harmonics texture, unit cell, atom parameters
     409        elif PawleyRef:
     410            pawleyVary = []
     411            for i,refl in enumerate(PawleyRef):
     412                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
     413                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
     414                if refl[5]:
     415                    pawleyVary.append(pfx+'PWLref:'+str(i))
     416            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
     417            phaseVary += pawleyVary   
     418               
     419    return phaseVary,phaseDict,pawleyLookup
     420   
     421def SetPhaseData(parmDict,sigDict,Phases):
     422   
     423    def cellFill(pfx,SGData,parmDict,sigDict):
     424        if SGData['SGLaue'] in ['-1',]:
     425            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
     426                parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
     427            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
     428                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
     429        elif SGData['SGLaue'] in ['2/m',]:
     430            if SGData['SGUniq'] == 'a':
     431                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
     432                    parmDict[pfx+'A3'],0,0]
     433                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
     434                    sigDict[pfx+'A3'],0,0]
     435            elif SGData['SGUniq'] == 'b':
     436                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
     437                    0,parmDict[pfx+'A4'],0]
     438                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
     439                    0,sigDict[pfx+'A4'],0]
     440            else:
     441                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
     442                    0,0,parmDict[pfx+'A5']]
     443                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
     444                    0,0,sigDict[pfx+'A5']]
     445        elif SGData['SGLaue'] in ['mmm',]:
     446            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
     447            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
     448        elif SGData['SGLaue'] in ['4/m','4/mmm']:
     449            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
     450            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
     451        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
     452            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
     453                parmDict[pfx+'A0'],0,0]
     454            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
     455        elif SGData['SGLaue'] in ['3R', '3mR']:
     456            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
     457                parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
     458            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
     459        elif SGData['SGLaue'] in ['m3m','m3']:
     460            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
     461            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
     462        return A,sigA
     463           
     464    print '\n Phases:'
     465    for phase in Phases:
     466        print ' Result for phase: ',phase
     467        Phase = Phases[phase]
     468        General = Phase['General']
     469        SGData = General['SGData']
     470        cell = General['Cell']
     471        pId = Phase['pId']
     472        pfx = str(pId)+'::'
     473        if cell[0]:
     474            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
     475            print ' Reciprocal metric tensor: '
     476            ptfmt = "%15.9f"
     477            names = ['A11','A22','A33','A12','A13','A23']
     478            namstr = '  names :'
     479            ptstr =  '  values:'
     480            sigstr = '  esds  :'
     481            for name,a,siga in zip(names,A,sigA):
     482                namstr += '%15s'%(name)
     483                ptstr += ptfmt%(a)
     484                if siga:
     485                    sigstr += ptfmt%(siga)
     486                else:
     487                    sigstr += 15*' '
     488            print namstr
     489            print ptstr
     490            print sigstr
     491            cell[1:7] = G2lat.A2cell(A)
     492            cell[7] = G2lat.calc_V(A)
     493            print ' New unit cell: a =%.5f'%(cell[1]),' b =%.5f'%(cell[2]), \
     494                ' c =%.5f'%(cell[3]),' alpha =%.3f'%(cell[4]),' beta =%.3f'%(cell[5]), \
     495                ' gamma =%.3f'%(cell[6]),' volume =%.3f'%(cell[7])
     496        if 'Pawley' in Phase['General']['Type']:
     497            pawleyRef = Phase['Pawley ref']
     498            for i,refl in enumerate(pawleyRef):
     499                key = pfx+'PWLref:'+str(i)
     500                refl[6] = parmDict[key]
     501                if key in sigDict:
     502                    refl[7] = sigDict[key]
     503                else:
     504                    refl[7] = 0
     505
     506def GetHistogramPhaseData(Phases,Histograms):
     507   
     508    def PrintSize(hapData):
     509        line = '\n Size model    : '+hapData[0]
     510        if hapData[0] in ['isotropic','uniaxial']:
     511            line += ' equatorial:'+'%12.2f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
     512            if hapData[0] == 'uniaxial':
     513                line += ' axial:'+'%12.2f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
     514            print line
     515        else:
     516            print line
     517            Snames = ['S11','S22','S33','S12','S13','S23']
     518            ptlbls = ' names :'
     519            ptstr =  ' values:'
     520            varstr = ' refine:'
     521            for i,name in enumerate(Snames):
     522                ptlbls += '%12s' % (name)
     523                ptstr += '%12.6f' % (hapData[4][i])
     524                varstr += '%12s' % (str(hapData[5][i]))
     525            print ptlbls
     526            print ptstr
     527            print varstr
     528       
     529    def PrintMuStrain(hapData,SGData):
     530        line = '\n Mustrain model: '+hapData[0]
     531        if hapData[0] in ['isotropic','uniaxial']:
     532            line += ' equatorial:'+'%12.4f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
     533            if hapData[0] == 'uniaxial':
     534                line += ' axial:'+'%12.4f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
     535            print line
     536        else:
     537            print line
     538            Snames = G2spc.MustrainNames(SGData)
     539            ptlbls = ' names :'
     540            ptstr =  ' values:'
     541            varstr = ' refine:'
     542            for i,name in enumerate(Snames):
     543                ptlbls += '%12s' % (name)
     544                ptstr += '%12.6f' % (hapData[4][i])
     545                varstr += '%12s' % (str(hapData[5][i]))
     546            print ptlbls
     547            print ptstr
     548            print varstr
     549       
     550   
     551    hapDict = {}
     552    hapVary = []
     553    controlDict = {}
     554    poType = {}
     555    poAxes = {}
     556    spAxes = {}
     557    spType = {}
     558   
     559    for phase in Phases:
     560        HistoPhase = Phases[phase]['Histograms']
     561        SGData = Phases[phase]['General']['SGData']
     562        cell = Phases[phase]['General']['Cell'][1:7]
     563        A = G2lat.cell2A(cell)
     564        pId = Phases[phase]['pId']
     565        for histogram in HistoPhase:
     566            print '\n Phase: ',phase,' in histogram: ',histogram
     567            hapData = HistoPhase[histogram]
     568            Histogram = Histograms[histogram]
     569            hId = Histogram['hId']
     570            limits = Histogram['Limits'][1]
     571            inst = Histogram['Instrument Parameters']
     572            inst = dict(zip(inst[3],inst[1]))
     573            Zero = inst['Zero']
     574            if 'C' in inst['Type']:
     575                try:
     576                    wave = inst['Lam']
     577                except KeyError:
     578                    wave = inst['Lam1']
     579                dmin = wave/(2.0*sind(limits[1]/2.0))
     580            pfx = str(pId)+':'+str(hId)+':'
     581            for item in ['Scale','Extinction']:
     582                hapDict[pfx+item] = hapData[item][0]
     583                if hapData[item][1]:
     584                    hapVary.append(pfx+item)
     585            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
     586            if hapData['Pref.Ori.'][0] == 'MD':
     587                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
     588                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
     589                if hapData['Pref.Ori.'][2]:
     590                    hapVary.append(pfx+'MD')
     591            else:                           #'SH' spherical harmonics
     592                for item in hapData['Pref.Ori.'][5]:
     593                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
     594                    if hapData['Pref.Ori.'][2]:
     595                        hapVary.append(pfx+item)                   
     596            for item in ['Mustrain','Size']:
     597                controlDict[pfx+item+'Type'] = hapData[item][0]
     598                if hapData[item][0] in ['isotropic','uniaxial']:
     599                    hapDict[pfx+item+':0'] = hapData[item][1][0]
     600                    if hapData[item][2][0]:
     601                        hapVary.append(pfx+item+':0')
     602                    if hapData[item][0] == 'uniaxial':
     603                        controlDict[pfx+item+'Axis'] = hapData[item][3]
     604                        hapDict[pfx+item+':1'] = hapData[item][1][1]
     605                        if hapData[item][2][1]:
     606                            hapVary.append(pfx+item+':1')
     607                else:       #generalized for mustrain or ellipsoidal for size
     608                    if item == 'Mustrain':
     609                        names = G2spc.MustrainNames(SGData)
     610                        pwrs = []
     611                        for name in names:
     612                            h,k,l = name[1:]
     613                            pwrs.append([int(h),int(k),int(l)])
     614                        controlDict[pfx+'MuPwrs'] = pwrs
     615                    for i in range(len(hapData[item][4])):
     616                        sfx = ':'+str(i)
     617                        hapDict[pfx+item+sfx] = hapData[item][4][i]
     618                        if hapData[item][5][i]:
     619                            hapVary.append(pfx+item+sfx)
     620                           
     621            print '\n Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
     622            print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
     623            if hapData['Pref.Ori.'][0] == 'MD':
     624                Ax = hapData['Pref.Ori.'][3]
     625                print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
     626                    ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])               
     627            PrintSize(hapData['Size'])
     628            PrintMuStrain(hapData['Mustrain'],SGData)
     629            HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
     630            refList = []
     631            for h,k,l,d in HKLd:
     632                ext,mul,Uniq = G2spc.GenHKLf([h,k,l],SGData)
     633                if ext:
     634                    continue
     635                if 'C' in inst['Type']:
     636                    pos = 2.0*asind(wave/(2.0*d))
     637                    if limits[0] < pos < limits[1]:
     638                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,Uniq])
     639                else:
     640                    raise ValueError
     641            Histogram['Reflection Lists'][phase] = refList
     642    return hapVary,hapDict,controlDict
     643   
     644def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms):
     645   
     646    def PrintSizeAndSig(hapData,sizeSig):
     647        line = '\n Size model: '+hapData[0]
     648        if hapData[0] in ['isotropic','uniaxial']:
     649            line += ' equatorial:%12.2f'%(hapData[1][0])
     650            if sizeSig[0][0]:
     651                line += ', sig: %8.2f'%(sizeSig[0][0])
     652            if hapData[0] == 'uniaxial':
     653                line += ' axial:%12.2f'%(hapData[1][1])
     654                if sizeSig[0][1]:
     655                    line += ', sig: %8.2f'%(sizeSig[0][1])
     656            print line
     657        else:
     658            print line
     659            Snames = ['S11','S22','S33','S12','S13','S23']
     660            ptlbls = ' name  :'
     661            ptstr =  ' value :'
     662            sigstr = ' sig   :'
     663            for i,name in enumerate(Snames):
     664                ptlbls += '%12s' % (name)
     665                ptstr += '%12.6f' % (hapData[4][i])
     666                if sizeSig[1][i]:
     667                    sigstr += '%12.6f' % (sizeSig[1][i])
     668                else:
     669                    sigstr += 12*' '
     670            print ptlbls
     671            print ptstr
     672            print sigstr
     673       
     674    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
     675        line = '\n Mustrain model: '+hapData[0]
     676        if hapData[0] in ['isotropic','uniaxial']:
     677            line += ' equatorial:%12.4f'%(hapData[1][0])
     678            if mustrainSig[0][0]:
     679                line += ',sig: %12.4f'%(mustrainSig[0][0])
     680            if hapData[0] == 'uniaxial':
     681                line += ' axial:%12.4f'%(hapData[1][1])
     682                if mustrainSig[0][1]:
     683                     line += ',sig: %12.4f'%(mustrainSig[0][1])
     684            print line
     685        else:
     686            print line
     687            Snames = G2spc.MustrainNames(SGData)
     688            ptlbls = ' name  :'
     689            ptstr =  ' value :'
     690            sigstr = ' sig   :'
     691            for i,name in enumerate(Snames):
     692                ptlbls += '%12s' % (name)
     693                ptstr += '%12.6f' % (hapData[4][i])
     694                sigstr += '%12.6f' % (mustrainSig[1][i])
     695            print ptlbls
     696            print ptstr
     697            print sigstr
     698       
     699    for phase in Phases:
     700        HistoPhase = Phases[phase]['Histograms']
     701        SGData = Phases[phase]['General']['SGData']
     702        pId = Phases[phase]['pId']
     703        for histogram in HistoPhase:
     704            print '\n Phase: ',phase,' in histogram: ',histogram
     705            hapData = HistoPhase[histogram]
     706            Histogram = Histograms[histogram]
     707            hId = Histogram['hId']
     708            pfx = str(pId)+':'+str(hId)+':'
     709           
     710# Add this stuff for Rietveld refinement!!!!           
     711#            for item in ['Scale','Extinction']:
     712#                hapDict[pfx+item] = hapData[item][0]
     713#                if hapData[item][1]:
     714#                    hapVary.append(pfx+item)
     715#            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
     716#            if hapData['Pref.Ori.'][0] == 'MD':
     717#                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
     718#                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
     719#                if hapData['Pref.Ori.'][2]:
     720#                    hapVary.append(pfx+'MD')
     721#            else:                           #'SH' spherical harmonics
     722#                for item in hapData['Pref.Ori.'][5]:
     723#                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
     724#                    if hapData['Pref.Ori.'][2]:
     725#                        hapVary.append(pfx+item)
     726#            print '\n Phase fraction  : %10.4f'%(hapData['Scale'][0])
     727#            print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0])
     728#            if hapData['Pref.Ori.'][0] == 'MD':
     729#                Ax = hapData['Pref.Ori.'][3]
     730#                print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
     731#                    ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
     732               
     733            SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
     734                'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]]}                 
     735            for item in ['Mustrain','Size']:
     736                if hapData[item][0] in ['isotropic','uniaxial']:
     737                    hapData[item][1][0] = parmDict[pfx+item+':0']
     738                    if hapData[item][2][0]:
     739                        SizeMuStrSig[item][0][0] = sigDict[pfx+item+':0']
     740                    if hapData[item][0] == 'uniaxial':
     741                        hapData[item][1][1] = parmDict[pfx+item+':1']
     742                        if hapData[item][2][1]:
     743                            SizeMuStrSig[item][0][1] = sigDict[pfx+item+':1']
     744                else:       #generalized for mustrain or ellipsoidal for size
     745                    for i in range(len(hapData[item][4])):
     746                        sfx = ':'+str(i)
     747                        hapData[item][4][i] = parmDict[pfx+item+sfx]
     748                        if hapData[item][5][i]:
     749                            SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
     750                           
     751            PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
     752            PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
     753   
     754def GetHistogramData(Histograms):
     755   
     756    def GetBackgroundParms(hId,Background):
     757        bakType,bakFlag = Background[:2]
     758        backVals = Background[3:]
     759        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
     760        if bakFlag:                                 #returns backNames as varyList = backNames
     761            return bakType,dict(zip(backNames,backVals)),backNames
     762        else:                                       #no background varied; varyList = []
     763            return bakType,dict(zip(backNames,backVals)),[]
     764       
     765    def GetInstParms(hId,Inst):
     766        insVals,insFlags,insNames = Inst[1:4]
     767        dataType = insVals[0]
     768        instDict = {}
     769        insVary = []
     770        pfx = ':'+str(hId)+':'
     771        for i,flag in enumerate(insFlags):
     772            insName = pfx+insNames[i]
     773            instDict[insName] = insVals[i]
     774            if flag:
     775                insVary.append(insName)
     776        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.1)
     777        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.1)
     778        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.001)
     779        return dataType,instDict,insVary
     780       
     781    def GetSampleParms(hId,Sample):
     782        sampVary = []
     783        pfx = ':'+str(hId)+':'       
     784        sampDict = {pfx+'Gonio. radius':Sample['Gonio. radius']}
     785        Type = Sample['Type']
     786        if 'Bragg' in Type:             #Bragg-Brentano
     787            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
     788                sampDict[pfx+item] = Sample[item][0]
     789                if Sample[item][1]:
     790                    sampVary.append(pfx+item)
     791        elif 'Debye' in Type:        #Debye-Scherrer
     792            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
     793                sampDict[pfx+item] = Sample[item][0]
     794                if Sample[item][1]:
     795                    sampVary.append(pfx+item)
     796        return Type,sampDict,sampVary
     797       
     798    def PrintBackground(Background):
     799        print '\n Background function: ',Background[0],' Refine?',bool(Background[1])
     800        line = ' Coefficients: '
     801        for back in Background[3:]:
     802            line += '%10.3f'%(back)
     803        print line
     804       
     805    def PrintInstParms(Inst):
     806        print '\n Instrument Parameters:'
     807        ptlbls = ' name  :'
     808        ptstr =  ' value :'
     809        varstr = ' refine:'
     810        instNames = Inst[3][1:]
     811        for i,name in enumerate(instNames):
     812            ptlbls += '%12s' % (name)
     813            ptstr += '%12.6f' % (Inst[1][i+1])
     814            if name in ['Lam1','Lam2','Azimuth']:
     815                varstr += 12*' '
     816            else:
     817                varstr += '%12s' % (str(bool(Inst[2][i+1])))
     818        print ptlbls
     819        print ptstr
     820        print varstr
     821       
     822    def PrintSampleParms(Sample):
     823        print '\n Sample Parameters:'
     824        ptlbls = ' name  :'
     825        ptstr =  ' value :'
     826        varstr = ' refine:'
     827        if 'Bragg' in Sample['Type']:
     828            for item in ['Scale','Shift','Transparency']:
     829                ptlbls += '%14s'%(item)
     830                ptstr += '%14.4f'%(Sample[item][0])
     831                varstr += '%14s'%(str(bool(Sample[item][1])))
     832           
     833        elif 'Debye' in Type:        #Debye-Scherrer
     834            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
     835                ptlbls += '%14s'%(item)
     836                ptstr += '%14.4f'%(Sample[item][0])
     837                varstr += '%14s'%(str(bool(Sample[item][1])))
     838
     839        print ptlbls
     840        print ptstr
     841        print varstr
     842       
     843
     844    histDict = {}
     845    histVary = []
     846    controlDict = {}
     847    for histogram in Histograms:
     848        Histogram = Histograms[histogram]
     849        hId = Histogram['hId']
     850        pfx = ':'+str(hId)+':'
     851        controlDict[pfx+'Limits'] = Histogram['Limits'][1]
     852       
     853        Background = Histogram['Background'][0]
     854        Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
     855        controlDict[pfx+'bakType'] = Type
     856        histDict.update(bakDict)
     857        histVary += bakVary
     858       
     859        Inst = Histogram['Instrument Parameters']
     860        Type,instDict,insVary = GetInstParms(hId,Inst)
     861        controlDict[pfx+'histType'] = Type
     862        histDict.update(instDict)
     863        histVary += insVary
     864       
     865        Sample = Histogram['Sample Parameters']
     866        Type,sampDict,sampVary = GetSampleParms(hId,Sample)
     867        controlDict[pfx+'instType'] = Type
     868        histDict.update(sampDict)
     869        histVary += sampVary
     870
     871        print '\n Histogram: ',histogram,' histogram Id: ',hId
     872        print 135*'-'
     873        Units = {'C':' deg','T':' msec'}
     874        units = Units[controlDict[pfx+'histType'][2]]
     875        Limits = controlDict[pfx+'Limits']
     876        print ' Instrument type: ',Sample['Type']
     877        print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
     878        PrintSampleParms(Sample)
     879        PrintInstParms(Inst)
     880        PrintBackground(Background)
     881       
     882    return histVary,histDict,controlDict
     883   
     884def SetHistogramData(parmDict,sigDict,Histograms):
     885   
     886    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
     887        backVals = Background[3:]
     888        backSig = [0 for i in range(len(Background[3:]))]
     889        for i in range(len(backVals)):
     890            backVals[i] = parmDict[pfx+'Back:'+str(i)]
     891            if Background[1]:
     892                backSig[i] = sigDict[pfx+'Back:'+str(i)]
     893        return backSig
     894       
     895    def SetInstParms(pfx,Inst,parmDict,sigDict):
     896        insVals,insFlags,insNames = Inst[1:4]
     897        instSig = [0 for i in range(len(insVals))]
     898        for i,flag in enumerate(insFlags):
     899            insName = pfx+insNames[i]
     900            insVals[i] = parmDict[insName]
     901            if flag:
     902                instSig[i] = sigDict[insName]
     903        return instSig
     904       
     905    def SetSampleParms(pfx,Sample,parmDict,sigDict):
     906        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
     907            sampSig = [0 for i in range(3)]
     908            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
     909                Sample[item][0] = parmDict[pfx+item]
     910                if Sample[item][1]:
     911                    sampSig[i] = sigDict[pfx+item]
     912        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
     913            sampSig = [0 for i in range(4)]
     914            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
     915                Sample[item][0] = parmDict[pfx+item]
     916                if Sample[item][1]:
     917                    sampSig[i] = sigDict[pfx+item]
     918        return sampSig
     919       
     920    def PrintBackgroundSig(Background,backSig):
     921        print '\n Background function: ',Background[0]
     922        valstr = ' value : '
     923        sigstr = ' sig   : '
     924        for i,back in enumerate(Background[3:]):
     925            valstr += '%10.4f'%(back)
     926            if Background[1]:
     927                sigstr += '%10.4f'%(backSig[i])
     928            else:
     929                sigstr += 10*' '
     930        print valstr
     931        print sigstr
     932       
     933    def PrintInstParmsSig(Inst,instSig):
     934        print '\n Instrument Parameters:'
     935        ptlbls = ' names :'
     936        ptstr =  ' value :'
     937        sigstr = ' sig   :'
     938        instNames = Inst[3][1:]
     939        for i,name in enumerate(instNames):
     940            ptlbls += '%12s' % (name)
     941            ptstr += '%12.6f' % (Inst[1][i+1])
     942            if instSig[i+1]:
     943                sigstr += '%12.6f' % (instSig[i+1])
     944            else:
     945                sigstr += 12*' '
     946        print ptlbls
     947        print ptstr
     948        print sigstr
     949       
     950    def PrintSampleParmsSig(Sample,sampleSig):
     951        print '\n Sample Parameters:'
     952        ptlbls = ' names :'
     953        ptstr =  ' values:'
     954        sigstr = ' sig   :'
     955        if 'Bragg' in Sample['Type']:
     956            for i,item in enumerate(['Scale','Shift','Transparency']):
     957                ptlbls += '%14s'%(item)
     958                ptstr += '%14.4f'%(Sample[item][0])
     959                sigstr += '%14.4f'%(sampleSig[i])
     960           
     961        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
     962            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
     963                ptlbls += '%14s'%(item)
     964                ptstr += '%14.4f'%(Sample[item][0])
     965                sigstr += '%14.4f'%(sampleSig[i])
     966
     967        print ptlbls
     968        print ptstr
     969        print sigstr
     970       
     971    for histogram in Histograms:
     972        if 'PWDR' in histogram:
     973            Histogram = Histograms[histogram]
     974            hId = Histogram['hId']
     975            pfx = ':'+str(hId)+':'
     976           
     977            Background = Histogram['Background'][0]
     978            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
     979           
     980            Inst = Histogram['Instrument Parameters']
     981            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
     982       
     983            Sample = Histogram['Sample Parameters']
     984            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
     985
     986            print '\n Histogram: ',histogram,' histogram Id: ',hId
     987            print 135*'-'
     988            print ' Instrument type: ',Sample['Type']
     989            PrintSampleParmsSig(Sample,sampSig)
     990            PrintInstParmsSig(Inst,instSig)
     991            PrintBackgroundSig(Background,backSig)
     992
     993
    229994def SetupSFcalc(General,Atoms):
    230995    ''' setup data for use in StructureFactor; mostly rearranging arrays
     
    3031068        print 'Absent'
    3041069       
     1070def Dict2Values(parmdict, varylist):
     1071    '''Use before call to leastsq to setup list of values for the parameters
     1072    in parmdict, as selected by key in varylist'''
     1073    return [parmdict[key] for key in varylist]
     1074   
     1075def Values2Dict(parmdict, varylist, values):
     1076    ''' Use after call to leastsq to update the parameter dictionary with
     1077    values corresponding to keys in varylist'''
     1078    parmdict.update(zip(varylist,values))
     1079   
     1080def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
     1081   
     1082    def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
     1083        if calcControls[phfx+'SizeType'] == 'isotropic':
     1084            gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*cosd(refl[5]/2.))
     1085        elif calcControls[phfx+'SizeType'] == 'uniaxial':
     1086            H = np.array(refl[:3])
     1087            P = np.array(calcControls[phfx+'SizeAxis'])
     1088            cosP,sinP = G2lat.CosSinAngle(H,V,P)
     1089            gam = (1.8*wave/np.pi)*parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1']
     1090            gam /= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2)*cosd(refl[5]/2.)
     1091        else:           #ellipsoidal crystallites
     1092            H = np.array(refl[:3])
     1093            gam += 1.8*wave/(np.pi*cosd(refl[5])*np.inner(H,np.inner(sizeEllipse,H)))           
     1094        if calcControls[phfx+'MustrainType'] == 'isotropic':
     1095            gam += 0.018*parmDict[phfx+'Mustrain:0']*tand(refl[5]/2.)/np.pi
     1096        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
     1097            H = np.array(refl[:3])
     1098            P = np.array(calcControls[phfx+'MustrainAxis'])
     1099            cosP,sinP = G2lat.CosSinAngle(H,V,P)
     1100            Si = parmDict[phfx+'Mustrain:0']
     1101            Sa = parmDict[phfx+'Mustrain:1']
     1102            gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
     1103        else:       #generalized - P.W. Stephens model
     1104            pwrs = calcControls[phfx+'MuPwrs']
     1105            sum = 0
     1106            for i,pwr in enumerate(pwrs):
     1107                sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
     1108            gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum
     1109           
     1110        return gam 
     1111       
     1112       
     1113    hId = Histogram['hId']
     1114    hfx = ':%d:'%(hId)
     1115    bakType = calcControls[hfx+'bakType']
     1116    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
     1117    yc = np.zeros_like(yb)
     1118       
     1119    if 'C' in calcControls[hfx+'histType']:   
     1120        U = parmDict[hfx+'U']
     1121        V = parmDict[hfx+'V']
     1122        W = parmDict[hfx+'W']
     1123        X = parmDict[hfx+'X']
     1124        Y = parmDict[hfx+'Y']
     1125        shl = max(parmDict[hfx+'SH/L'],0.001)
     1126        Zero = parmDict[hfx+'Zero']
     1127        pola = parmDict[hfx+'Polariz.']
     1128        Ka2 = False
     1129        if hfx+'Lam1' in parmDict.keys():
     1130            wave = parmDict[hfx+'Lam1']
     1131            Ka2 = True
     1132            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
     1133            kRatio = parmDict[hfx+'I(L2)/I(L1)']
     1134        else:
     1135            wave = parmDict[hfx+'Lam']
     1136    else:
     1137        print 'TOF Undefined at present'
     1138        raise ValueError
     1139    for phase in Histogram['Reflection Lists']:
     1140        refList = Histogram['Reflection Lists'][phase]
     1141        Phase = Phases[phase]
     1142        pId = Phase['pId']
     1143        pfx = '%d::'%(pId)
     1144        phfx = '%d:%d:'%(pId,hId)
     1145        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
     1146        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
     1147        sizeEllipse = []
     1148        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
     1149            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
     1150        for refl in refList:
     1151            if 'C' in calcControls[hfx+'histType']:
     1152                h,k,l,m = refl[:4]
     1153                dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
     1154                d = np.sqrt(dsq)
     1155                pos = 2.0*asind(wave/(2.0*d))+Zero
     1156                const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
     1157                if 'Bragg' in calcControls[hfx+'instType']:
     1158                    pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
     1159                        1.e-7*parmDict[hfx+'Transparency']*sind(pos))            #trans(=1/mueff) in Angstroms
     1160                else:               #Debye-Scherrer - simple but maybe not right
     1161                    pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
     1162                refl[5] = pos
     1163                tanPos = tand(pos/2.0)
     1164                refl[6] = U*tanPos**2+V*tanPos+W        #save peak sigma
     1165                refl[7] = X/cosd(pos/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma
     1166                if 'Pawley' in Phase['General']['Type']:
     1167                    try:
     1168                        refl[8] = parmDict[hfx+'Scale']*m*parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]
     1169                    except KeyError:
     1170                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
     1171                        continue
     1172                else:
     1173                    raise ValueError       #wants strctrfacr calc here
     1174                Wd,fmin,fmax = G2pwd.getWidths(pos,refl[6],refl[7],shl)
     1175                if pos > 90:
     1176                    fmin,fmax = [fmax,fmin]         
     1177                iBeg = np.searchsorted(x,pos-fmin)
     1178                iFin = np.searchsorted(x,pos+fmax)
     1179                if not iBeg+iFin:       #peak below low limit - skip peak
     1180                    continue
     1181                elif not iBeg-iFin:     #peak above high limit - done
     1182                    return yc,yb
     1183                yc[iBeg:iFin] += G2pwd.getFCJVoigt(pos,refl[8],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
     1184                if Ka2:
     1185                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
     1186                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
     1187                    if pos > 90:
     1188                        fmin,fmax = [fmax,fmin]         
     1189                    iBeg = np.searchsorted(x,pos2-fmin)
     1190                    iFin = np.searchsorted(x,pos2+fmax)
     1191                    yc[iBeg:iFin] += G2pwd.getFCJVoigt(pos2,refl[8]*kRatio,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
     1192            else:
     1193                raise ValueError
     1194    return yc,yb   
     1195           
     1196
     1197   
    3051198def Refine(GPXfile):
     1199   
     1200    def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
     1201        parmdict.update(zip(varylist,values))
     1202        G2mv.Dict2Map(parmDict)
     1203        Histograms,Phases = HistoPhases
     1204        M = np.empty(0)
     1205        sumwYo = 0
     1206        Nobs = 0
     1207        for histogram in Histograms:
     1208            if 'PWDR' in histogram[:4]:
     1209                Histogram = Histograms[histogram]
     1210                hId = Histogram['hId']
     1211                hfx = ':%d:'%(hId)
     1212                Limits = calcControls[hfx+'Limits']
     1213                x,y,w,yc,yb,yd = Histogram['Data']
     1214                yc *= 0.0                           #zero full calcd profiles
     1215                yb *= 0.0
     1216                yd *= 0.0
     1217                xB = np.searchsorted(x,Limits[0])
     1218                xF = np.searchsorted(x,Limits[1])
     1219                Nobs += xF-xB
     1220                Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
     1221                sumwYo += Histogram['sumwYo']
     1222                yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
     1223                    varylist,Histogram,Phases,calcControls,pawleyLookup)
     1224                yc[xB:xF] *= parmDict[hfx+'Scale']
     1225                yc[xB:xF] += yb[xB:xF]
     1226                yd[xB:xF] = y[xB:xF]-yc[xB:xF]
     1227                Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
     1228                M = np.concatenate((M,np.sqrt(w[xB:xF])*(yd[xB:xF])))
     1229#                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])
     1230        Histograms['sumwYo'] = sumwYo
     1231        Histograms['Nobs'] = Nobs
     1232        Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
     1233        if dlg:
     1234            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
     1235            if not GoOn:
     1236                return -M           #abort!!
     1237        return M
     1238   
    3061239    ShowBanner()
    307     Controls = GetControls(GPXfile)
    308     ShowControls(Controls)
    309    
    310     phaseNames = GetPhaseNames(GPXfile)
    311     phaseData = {}
    312     for name in phaseNames:
    313         phaseData[name] =  GetPhaseData(GPXfile,name)
    314     ShowPhaseData(phaseNames,phaseData)
    315        
    316     histograms = GetHistogramNames(GPXfile)
    317     for hist in histograms:
    318         if 'PWDR' in hist[:4]:
    319             print hist
    320             PWDRdata = GetPWDRdata(GPXfile,hist)
    321             print PWDRdata
     1240    varyList = []
     1241    parmDict = {}
     1242    calcControls = {}   
     1243#    Controls = GetControls(GPXfile)
     1244#    ShowControls(Controls)           
     1245    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
     1246    if not Phases:
     1247        print ' *** ERROR - you have no phases to refine! ***'
     1248        print ' *** Refine aborted ***'
     1249        raise Exception
     1250    if not Histograms:
     1251        print ' *** ERROR - you have no data to refine with! ***'
     1252        print ' *** Refine aborted ***'
     1253        raise Exception
     1254    phaseVary,phaseDict,pawleyLookup = GetPhaseData(Phases)
     1255    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
     1256    calcControls.update(controlDict)
     1257    histVary,histDict,controlDict = GetHistogramData(Histograms)
     1258    calcControls.update(controlDict)
     1259    varyList = phaseVary+histVary+hapVary
     1260    parmDict.update(phaseDict)
     1261    parmDict.update(hapDict)
     1262    parmDict.update(histDict)
     1263    constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
     1264    groups,parmlist = G2mv.GroupConstraints(constrDict)
     1265    G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList)
     1266    G2mv.Map2Dict(parmDict,varyList)
     1267
     1268    while True:
     1269        begin = time.time()
     1270        values =  np.array(Dict2Values(parmDict, varyList))
     1271        dlg = wx.ProgressDialog('Residual','Powder profile Rwp =',101.0,
     1272            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
     1273        screenSize = wx.ClientDisplayRect()
     1274        Size = dlg.GetSize()
     1275        dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
     1276        try:
     1277            result = so.leastsq(errRefine,values,full_output=True,  #factor=1.,epsfcn=0.00001,ftol=0.0001,
     1278                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     1279            print len(result)
     1280        finally:
     1281            dlg.Destroy()
     1282        runtime = time.time()-begin   
     1283        chisq = np.sum(result[2]['fvec']**2)
     1284        ncyc = int(result[2]['nfev']/len(varyList))
     1285        Values2Dict(parmDict, varyList, result[0])
     1286        G2mv.Dict2Map(parmDict)
     1287        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
     1288        GOF = chisq/(Histograms['Nobs']-len(varyList))
     1289        print '\n Refinement results:'
     1290        print 135*'-'
     1291        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
     1292        print 'Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
     1293        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
     1294        try:
     1295            sig = np.sqrt(np.diag(result[1])*GOF)
     1296            if np.any(np.isnan(sig)):
     1297                print '*** Least squares aborted - some invalid esds possible ***'
     1298            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
     1299#            for item in table: print item,table[item]               #useful debug - are things shifting?
     1300            break                   #refinement succeeded - finish up!
     1301        except ValueError:          #result[1] is None on singular matrix
     1302            print '**** Refinement failed - singular matrix ****'
     1303            Ipvt = result[2]['ipvt']
     1304            for i,ipvt in enumerate(Ipvt):
     1305                if not np.sum(result[2]['fjac'],axis=1)[i]:
     1306                    print 'Removing parameter: ',varyList[ipvt-1]
     1307                    del(varyList[ipvt-1])
     1308                    break
     1309
     1310    sigDict = dict(zip(varyList,sig))
     1311    SetPhaseData(parmDict,sigDict,Phases)
     1312    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
     1313    SetHistogramData(parmDict,sigDict,Histograms)
     1314    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases)
    3221315
    3231316def main():
Note: See TracChangeset for help on using the changeset viewer.