Changeset 345


Ignore:
Timestamp:
Aug 18, 2011 8:54:16 PM (12 years ago)
Author:
vondreele
Message:

put new banner in GSASII.BAT
fix for vertical tilt & offset images
speedup of peak calcs
allow change from lam to lam1 & lam2
Pawley refinement fixes

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASII.BAT

    r328 r345  
    1 @echo ==========================================================
    2 @echo                 This file starts GSAS-II
    3 @echo = Run this from the directory where GSAS-II is installed =
    4 @echo ==========================================================
     1@echo ========================================================================
     2@echo                General Structure Analysis System-II
     3@echo      by Robert B. Von Dreele, Argonne National Laboratory(C), 2010
     4@echo  This product includes software developed by the UChicago Argonne, LLC,
     5@echo             as Operator of Argonne National Laboratory.'
     6@echo ========================================================================
    57@
    6 @REM Get this script's directory; make sure that the path ends 
     8@REM Get this script's directory; make sure that the path ends
    79@REM    with a single backslash
    810@set gsasloc=%~dp0\*

  • trunk/GSASII.py

    r342 r345  
    13511351        try:
    13521352            if dlg.ShowModal() == wx.ID_OK:
     1353                Id = 0
    13531354                self.PatternTree.DeleteChildren(self.root)
    13541355                if self.HKL: self.HKL = []
    … …  
    13561357                    self.G2plotNB.clear()
    13571358                G2IO.ProjFileOpen(self)
     1359                item, cookie = self.PatternTree.GetFirstChild(self.root)
     1360                while item and not Id:
     1361                    name = self.PatternTree.GetItemText(item)
     1362                    if name[:4] in ['PWDR','HKLF']:
     1363                        Id = item
     1364                    item, cookie = self.PatternTree.GetNextChild(self.root, cookie)               
     1365                if Id:
     1366                    self.PatternTree.SelectItem(Id)
    13581367        finally:
    13591368            dlg.Destroy()

  • trunk/GSASIIimage.py

    r298 r345  
    180180    ValSig = zip(names,fmt,vals,sig)
    181181    CalibPrint(ValSig)
    182     if len(rings) > 1:
     182    try:
    183183        print 'Trial refinement of wavelength - not used for calibration'
    184184        p0 = result[0]
    … …  
    189189        ValSig = zip(names,fmt,resultW[0],sig*np.sqrt(np.diag(resultW[1])))
    190190        CalibPrint(ValSig)
    191         return result[0],resultW[0][-1]
    192     else:
     191        return result[0],resultW[0][-1]       
     192    except ValueError:
     193        print 'Bad wavelength refinement - no result'
    193194        return result[0],wave
    194195           
    … …  
    458459        radius = dist*ttth
    459460        elcent,phi,radii = ellipse
    460         radii[1] = dist*stth*ctth*cosB/(cosB**2-stth**2)
    461         radii[0] = math.sqrt(radii[1]*radius*cosB)
    462         zdis,cosB = calcZdisCosB(radius,tth,radii)
    463         zsinp = zdis*sind(phi)
    464         zcosp = zdis*cosd(phi)
    465         cent = data['center']
    466         elcent = [cent[0]+zsinp,cent[1]-zcosp]
     461        if i:
     462            radii[1] = dist*stth*ctth*cosB/(cosB**2-stth**2)
     463            radii[0] = math.sqrt(radii[1]*radius*cosB)
     464            zdis,cosB = calcZdisCosB(radius,tth,radii)
     465            zsinp = zdis*sind(phi)
     466            zcosp = zdis*cosd(phi)
     467            cent = data['center']
     468            elcent = [cent[0]-zsinp,cent[1]+zcosp]
     469            ellipse = (elcent,phi,radii)
    467470        ratio = radii[1]/radii[0]
    468471        Ring,delt = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
    … …  
    549552    data['tilt'] = result[4]
    550553    N = len(data['ellipses'])
    551     data['ellipses'] = []           #clear away individual ellipse fits
     554#    data['ellipses'] = []           #clear away individual ellipse fits
    552555    for H in HKL[:N]:
    553556        ellipse = GetEllipse(H[3],data)

  • trunk/GSASIIphsGUI.py

    r342 r345  
    320320            except ValueError:
    321321                pass
    322             pawlVal.SetValue("%.2f"%(generalData['Pawley dmin']))          #reset in case of error           
     322            pawlVal.SetValue("%.3f"%(generalData['Pawley dmin']))          #reset in case of error           
    323323                                   
    324324        cellGUIlist = [[['m3','m3m'],4,zip([" Unit cell: a = "," Vol = "],["%.5f","%.3f"],[True,False],[0,0])],
    … …  
    455455            pawlSizer = wx.BoxSizer(wx.HORIZONTAL)
    456456            pawlSizer.Add(wx.StaticText(dataDisplay,label=' Pawley dmin: '),0,wx.ALIGN_CENTER_VERTICAL)
    457             pawlVal = wx.TextCtrl(dataDisplay,value='%.2f'%(generalData['Pawley dmin']),style=wx.TE_PROCESS_ENTER)
     457            pawlVal = wx.TextCtrl(dataDisplay,value='%.3f'%(generalData['Pawley dmin']),style=wx.TE_PROCESS_ENTER)
    458458            pawlVal.Bind(wx.EVT_TEXT_ENTER,OnPawleyVal)       
    459459            pawlVal.Bind(wx.EVT_KILL_FOCUS,OnPawleyVal)

  • trunk/GSASIIpwd.py

    r343 r345  
    410410np.seterr(divide='ignore')
    411411
    412 # Voigt function = convolution (norm,cauchy)
    413 
    414 class voigt_gen(st.rv_continuous):
    415     '''
    416     Voigt function
    417     Parameters
    418     -----------------------------------------
    419     '''
    420     def _pdf(self,x,s,g):
    421         z = np.empty([2,len(x)])
    422         z[0] = st.norm.pdf(x,scale=s)
    423         z[1] = st.cauchy.pdf(x,scale=g)
    424         Z = fft(z)
    425         return fftshift(ifft(Z.prod(axis=0))).real
    426 #    def _cdf(self, x):
    427 #    def _ppf(self, q):
    428 #    def _sf(self, x):
    429 #    def _isf(self, q):
    430 #    def _stats(self):
    431 #    def _entropy(self):
    432 
    433 voigt = voigt_gen(name='voigt',shapes='ts,g')
    434        
    435 # Finger-Cox_Jephcoat D(2phi,2th) function for S/L = H/L
     412# Normal distribution
     413
     414# loc = mu, scale = std
     415_norm_pdf_C = 1./math.sqrt(2*math.pi)
     416class norm_gen(st.rv_continuous):
     417       
     418    def pdf(self,x,*args,**kwds):
     419        loc,scale=kwds['loc'],kwds['scale']
     420        x = (x-loc)/scale
     421        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
     422       
     423norm = norm_gen(name='norm',longname='A normal',extradoc="""
     424
     425Normal distribution
     426
     427The location (loc) keyword specifies the mean.
     428The scale (scale) keyword specifies the standard deviation.
     429
     430normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
     431""")
     432
     433## Cauchy
     434
     435# median = loc
     436
     437class cauchy_gen(st.rv_continuous):
     438
     439    def pdf(self,x,*args,**kwds):
     440        loc,scale=kwds['loc'],kwds['scale']
     441        x = (x-loc)/scale
     442        return 1.0/np.pi/(1.0+x*x) / scale
     443       
     444cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
     445
     446Cauchy distribution
     447
     448cauchy.pdf(x) = 1/(pi*(1+x**2))
     449
     450This is the t distribution with one degree of freedom.
     451""")
     452   
     453   
     454#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
     455
    436456
    437457class fcjde_gen(st.rv_continuous):
    … …  
    457477    def _pdf(self,x,t,s,dx):
    458478        T = dx*x+t
    459         ax = npcosd(T)**2
     479        ax2 = abs(npcosd(T))
     480        ax = ax2**2
    460481        bx = npcosd(t)**2
    461482        bx = np.where(ax>bx,bx,ax)
    462         fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/np.sqrt(ax),0.0)
     483        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
    463484        fx = np.where(fx > 0.,fx,0.0)
    464485        return fx
    465 #    def _cdf(self, x):
    466 #    def _ppf(self, q):
    467 #    def _sf(self, x):
    468 #    def _isf(self, q):
    469 #    def _stats(self):
    470 #    def _entropy(self):
     486             
     487    def pdf(self,x,*args,**kwds):
     488        loc=kwds['loc']
     489        return self._pdf(x-loc,*args)
    471490       
    472491fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
    … …  
    485504    return yb
    486505
    487 def 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                    
    498 def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
    499    
     506def getWidths(pos,sig,gam,shl):
     507    widths = [np.sqrt(sig)/100.,gam/200.]
     508    fwhm = 2.355*widths[0]+2.*widths[1]
     509    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
     510    fmax = 15.0*fwhm
     511    if pos > 90:
     512        fmin,fmax = [fmax,fmin]         
     513    return widths,fmin,fmax
     514               
     515def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
    500516    DX = xdata[1]-xdata[0]
    501517    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]         
     518    x = np.linspace(pos-fmin,pos+fmin,256)
    505519    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]
     520    Norm = norm.pdf(x,loc=pos,scale=widths[0])
     521    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
     522    arg = [pos,shl/57.2958,dx,]
     523    FCJ = fcjde.pdf(x,*arg,loc=pos)
     524    if len(np.nonzero(FCJ)[0])>5:
     525        z = np.column_stack([Norm,Cauchy,FCJ]).T
     526        Z = fft(z)
     527        Df = ifft(Z.prod(axis=0)).real
    511528    else:
    512         fxns = [st.norm,st.cauchy]
    513         args = [[],[]]
    514     Df = calcPeakFFT(x,fxns,widths,pos,args)
     529        z = np.column_stack([Norm,Cauchy]).T
     530        Z = fft(z)
     531        Df = fftshift(ifft(Z.prod(axis=0))).real
    515532    Df /= np.sum(Df)
    516533    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
    517534    return intens*Df(xdata)*DX/dx
    518                                    
     535
    519536def getPeakProfile(parmDict,xdata,varyList,bakType):
    520537   
    521538    yb = getBackground('',parmDict,bakType,xdata)
    522539    yc = np.zeros_like(yb)
     540    dx = xdata[1]-xdata[0]
    523541    U = parmDict['U']
    524542    V = parmDict['V']
    … …  
    526544    X = parmDict['X']
    527545    Y = parmDict['Y']
    528     shl = max(parmDict['SH/L'],0.001)
     546    shl = max(parmDict['SH/L'],0.0005)
    529547    Ka2 = False
    530548    if 'Lam1' in parmDict.keys():
    531549        Ka2 = True
    532         lamRatio = parmDict['Lam2']/parmDict['Lam1']
     550        lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
    533551        kRatio = parmDict['I(L2)/I(L1)']
    534552    iPeak = 0
    … …  
    548566            else:
    549567                gam = X/cosd(pos/2.0)+Y*tand(pos/2.0)
    550             gam = max(gam,0.01)             #avoid neg gamma
     568            gam = max(gam,0.001)             #avoid neg gamma
    551569            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
    552             if pos > 90:
    553                 fmin,fmax = [fmax,fmin]         
    554570            iBeg = np.searchsorted(xdata,pos-fmin)
    555             iFin = np.searchsorted(xdata,pos+fmax)
     571            lenX = len(xdata)
     572            if not iBeg:
     573                iFin = np.searchsorted(xdata,pos+fmin)
     574            elif iBeg == lenX:
     575                iFin = iBeg
     576            else:
     577                iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
    556578            if not iBeg+iFin:       #peak below low limit
    557579                iPeak += 1
    … …  
    561583            yc[iBeg:iFin] += getFCJVoigt(pos,intens,sig,gam,shl,xdata[iBeg:iFin])
    562584            if Ka2:
    563                 pos2 = 2.0*asind(lamRatio*sind(pos/2.0))
    564                 Wd,fmin,fmax = getWidths(pos2,sig,gam,shl)
    565                 if pos > 90:
    566                     fmin,fmax = [fmax,fmin]         
    567                 iBeg = np.searchsorted(xdata,pos2-fmin)
    568                 iFin = np.searchsorted(xdata,pos2+fmax)
    569                 yc[iBeg:iFin] += getFCJVoigt(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin])
     585                pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
     586                kdelt = int((pos2-pos)/dx)               
     587                iBeg = min(lenX,iBeg+kdelt)
     588                iFin = min(lenX,iFin+kdelt)
     589                if iBeg-iFin:
     590                    yc[iBeg:iFin] += getFCJVoigt(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin])
    570591            iPeak += 1
    571592        except KeyError:        #no more peaks to process
    572593            return yb+yc
    573        
    574 def getWidths(pos,sig,gam,shl):
    575     if shl:
    576         widths = [np.sqrt(sig)/100.,gam/200.,.1]
    577     else:
    578         widths = [np.sqrt(sig)/100.,gam/200.]
    579     fwhm = 2.355*widths[0]+2.*widths[1]
    580     fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
    581     fmax = 15.0*fwhm
    582     return widths,fmin,fmax
    583                
     594
    584595def Dict2Values(parmdict, varylist):
    585596    '''Use before call to leastsq to setup list of values for the parameters
    … …  
    635646                insVary.append(insNames[i])
    636647        instDict = dict(zip(insNames,insVals))
    637         instDict['X'] = max(instDict['X'],0.1)
    638         instDict['Y'] = max(instDict['Y'],0.1)
    639         instDict['SH/L'] = max(instDict['SH/L'],0.001)
     648        instDict['X'] = max(instDict['X'],0.01)
     649        instDict['Y'] = max(instDict['Y'],0.01)
     650        instDict['SH/L'] = max(instDict['SH/L'],0.0001)
    640651        return dataType,instDict,insVary
    641652       
    … …  
    735746        return M
    736747       
    737     Ftol = 0.01
     748    Ftol = 0.0001
    738749    if oneCycle:
    739750        Ftol = 1.0
    … …  
    759770            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
    760771            try:
    761                 result = so.leastsq(errPeakProfile,values,full_output=True,ftol=Ftol,
     772                result = so.leastsq(errPeakProfile,values,full_output=True,             #ftol=Ftol,
    762773                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
    763774            finally:

  • trunk/GSASIIpwdGUI.py

    r344 r345  
    8484        OnPeakFit('LSQ',oneCycle=True)
    8585       
    86     def OnBGFSPeakFit(event):
    87         OnPeakFit('BFGS')
    88                
     86    def OnClearPeaks(event):
     87        dlg = wx.MessageDialog(self,'Delete all peaks?','Clear peak list',wx.OK|wx.CANCEL)
     88        try:
     89            if dlg.ShowModal() == wx.ID_OK:
     90                peaks = []
     91        finally:
     92            dlg.Destroy()
     93        UpdatePeakGrid(self,peaks)
     94        G2plt.PlotPatterns(self)
     95       
    8996    def OnPeakFit(FitPgm,oneCycle=False):
    9097        SaveState()
    … …  
    138145        for key in T: X.append(D[key])
    139146        data = X       
    140         G2plt.PlotPatterns(self)
    141147       
    142148    def setBackgroundColors():
    … …  
    148154                   else:
    149155                       self.dataDisplay.SetCellBackgroundColour(r,c,wx.WHITE)
    150                        
    151     def RefineSelect(event):
    152         data = self.PatternTree.GetItemPyData(self.PickId)
    153         r,c =  event.GetRow(),event.GetCol()
    154         if r < 0 and self.dataDisplay.GetColLabelValue(c) == 'refine':
    155             self.dataDisplay.SelectCol(c,False)
    156            
    157     def OnRefine(event):
    158         r,c =  event.GetRow(),event.GetCol()
    159         if self.dataDisplay.GetColLabelValue(c) == 'refine':
    160             if self.PeakTable.GetValue(r,c):
    161                 data[r][c] = False
    162             else:
    163                 data[r][c] = True
    164             print r,c,data[r][c]
    165        
    166                        
    167     def RowSelect(event):
    168         r,c =  event.GetRow(),event.GetCol()
    169         if r < 0 and c < 0:
    170             if self.dataDisplay.IsSelection():
    171                 self.dataDisplay.ClearSelection()
    172         elif c < 0:                   #only row clicks
    173             if event.ControlDown():                   
    174                 if r in self.dataDisplay.GetSelectedRows():
    175                     self.dataDisplay.DeselectRow(r)
    176                 else:
    177                     self.dataDisplay.SelectRow(r,True)
    178             elif event.ShiftDown():
    179                 for row in range(r+1):
    180                     self.dataDisplay.SelectRow(row,True)
    181             else:
    182                 self.dataDisplay.ClearSelection()
    183                 self.dataDisplay.SelectRow(r,True)               
    184        
    185                            
     156                                                 
    186157    def KeyEditPeakGrid(event):
    187158        rowList = self.dataDisplay.GetSelectedRows()
    … …  
    241212    self.Bind(wx.EVT_MENU, OnLSQPeakFit, id=G2gd.wxID_LSQPEAKFIT)
    242213    self.Bind(wx.EVT_MENU, OnOneCycle, id=G2gd.wxID_LSQONECYCLE)
    243 #    self.Bind(wx.EVT_MENU, OnBGFSPeakFit, id=G2gd.wxID_BFGSPEAKFIT)
     214    self.Bind(wx.EVT_MENU, OnClearPeaks, id=G2gd.wxID_CLEARPEAKS)
    244215    self.Bind(wx.EVT_MENU, OnResetSigGam, id=G2gd.wxID_RESETSIGGAM)
    245216    self.dataFrame.PeakFit.Enable(False)
    … …  
    271242    self.dataDisplay.Bind(wg.EVT_GRID_CELL_CHANGE, RefreshPeakGrid)
    272243    self.dataDisplay.Bind(wx.EVT_KEY_DOWN, KeyEditPeakGrid)
    273     self.dataDisplay.Bind(wg.EVT_GRID_LABEL_LEFT_CLICK, RowSelect)                 
    274     self.dataDisplay.Bind(wg.EVT_GRID_LABEL_LEFT_DCLICK, RefineSelect)
    275     self.dataDisplay.Bind(wg.EVT_GRID_CELL_RIGHT_CLICK,OnRefine)
    276244    self.dataDisplay.SetMargins(0,0)
    277245    self.dataDisplay.AutoSizeColumns(False)
    … …  
    413381        UpdateInstrumentGrid(self,data)
    414382       
     383    def OnWaveChange(event):
     384        if 'Lam' in insVal:           
     385            data[0] = data[0][:1]+tuple(waves['CuKa'])+(.5,)+data[0][2:]
     386            data[1] = data[1][:1]+waves['CuKa']+[.5,]+data[1][2:]
     387            data[2] = data[2][:1]+[0,0,0,]+data[2][2:]
     388            data[3] = data[3][:1]+['Lam1','Lam2','I(L2)/I(L1)',]+data[3][2:]           
     389        else:
     390            data[0] = data[0][:2]+data[0][4:]
     391            data[1] = data[1][:2]+data[1][4:]
     392            data[2] = data[2][:2]+data[2][4:]
     393            data[3] = data[3][:1]+['Lam',]+data[3][4:]           
     394        UpdateInstrumentGrid(self,data)
     395               
    415396    def OnNewType(event):
    416397        insVal['Type'] = typePick.GetValue()
    … …  
    424405        data = updateData(insVal,insRef)
    425406        UpdateInstrumentGrid(self,data)
    426          
    427        
     407                 
    428408    def OnRatValue(event):
    429409        try:
    … …  
    461441        try:
    462442            value = float(Obj.GetValue())
    463             if value < 0:
    464                 raise ValueError
    465443        except ValueError:
    466444            value = insVal[item]
    … …  
    487465        if not self.dataFrame.GetStatusBar():
    488466            Status = self.dataFrame.CreateStatusBar()
    489         self.Bind(wx.EVT_MENU, OnReset, id=G2gd.wxID_INSTPRMRESET)
     467        self.Bind(wx.EVT_MENU, OnReset,id=G2gd.wxID_INSTPRMRESET)
     468        self.Bind(wx.EVT_MENU,OnWaveChange,id=G2gd.wxID_CHANGEWAVETYPE)
    490469        typePick = wx.ComboBox(self.dataDisplay,value=insVal['Type'],
    491470            choices=['PXC','PNC','PNT'],style=wx.CB_READONLY|wx.CB_DROPDOWN)
    … …  
    753732            key = event.GetKeyCode()
    754733            for col in colList:
    755                 if self.IndexPeaksTable.GetTypeName(0,col) == wg.GRID_VALUE_BOOL:
     734                if self.IndexPeaksTable.GetColLabelValue(col) in ['use','refine']:
    756735                    if key == 89: #'Y'
    757736                        for row in range(self.IndexPeaksTable.GetNumberRows()): data[row][col]=True
    … …  
    12911270    for i in range(len(refList)): rowLabels.append(str(i+1))
    12921271    colLabels = ['H','K','L','mul','d','pos','sig','gam','Fosq','Fcsq',]
    1293     Types = 4*[wg.GRID_VALUE_LONG,]+6*[wg.GRID_VALUE_FLOAT+':10,4',]
     1272    Types = 4*[wg.GRID_VALUE_LONG,]+4*[wg.GRID_VALUE_FLOAT+':10,4',]+2*[wg.GRID_VALUE_FLOAT+':10,2',]
    12941273    self.PeakTable = G2gd.Table(refList,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    12951274    self.dataFrame.SetLabel('Reflection List for '+self.RefList)

  • trunk/GSASIIstruct.py

    r342 r345  
    300300    if SGLaue in ['-1','2/m','mmm']:
    301301        return                      #no Pawley symmetry required constraints
    302     for varyI in pawleyVary:
     302    for i,varyI in enumerate(pawleyVary):
    303303        refI = int(varyI.split(':')[-1])
    304304        ih,ik,il = PawleyRef[refI][:3]
    305         for varyJ in pawleyVary[:refI]:
     305        for varyJ in pawleyVary[0:i]:
    306306            refJ = int(varyJ.split(':')[-1])
    307307            jh,jk,jl = PawleyRef[refJ][:3]
    … …  
    491491            cell[1:7] = G2lat.A2cell(A)
    492492            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])
     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])
    496496        if 'Pawley' in Phase['General']['Type']:
    497497            pawleyRef = Phase['Pawley ref']
    … …  
    708708            pfx = str(pId)+':'+str(hId)+':'
    709709           
    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])
     710            PhFrExtPOSig = {}
     711            for item in ['Scale','Extinction']:
     712                hapData[item][0] = parmDict[pfx+item]
     713                if hapData[item][1]:
     714                    PhFrExtPOSig[item] = sigDict[pfx+item]           
     715            if hapData['Pref.Ori.'][0] == 'MD':
     716                hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
     717                if hapData['Pref.Ori.'][2]:
     718                    PhFrExtPOSig[item] = sigDict[pfx+item]
     719            else:                           #'SH' spherical harmonics
     720                for item in hapData['Pref.Ori.'][5]:
     721                    hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
     722                    if hapData['Pref.Ori.'][2]:
     723                        PhFrExtPOSig[item] = sigDict[pfx+item]
     724#            print '\n Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
     725#            print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
    728726#            if hapData['Pref.Ori.'][0] == 'MD':
    729727#                Ax = hapData['Pref.Ori.'][3]
    … …  
    774772            if flag:
    775773                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)
     774        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.01)
     775        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.01)
     776        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
    779777        return dataType,instDict,insVary
    780778       
    … …  
    11061104            for i,pwr in enumerate(pwrs):
    11071105                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        
     1106            gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum           
     1107        return gam
     1108       
     1109    def GetIntensityCorr(refl,phfx,hfx,calcControls,parmDict):
     1110        Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
     1111        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
     1112       
     1113        return Icorr
     1114       
     1115    def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
     1116        h,k,l = refl[:3]
     1117        dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
     1118        d = np.sqrt(dsq)
     1119        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
     1120        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
     1121        if 'Bragg' in calcControls[hfx+'instType']:
     1122            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
     1123                1.e-7*parmDict[hfx+'Transparency']*sind(pos))            #trans(=1/mueff) in Angstroms
     1124        else:               #Debye-Scherrer - simple but maybe not right
     1125            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
     1126        return pos
     1127   
     1128    def GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse):
     1129        U = parmDict[hfx+'U']
     1130        V = parmDict[hfx+'V']
     1131        W = parmDict[hfx+'W']
     1132        X = parmDict[hfx+'X']
     1133        Y = parmDict[hfx+'Y']
     1134        tanPos = tand(refl[5]/2.0)
     1135        sig = U*tanPos**2+V*tanPos+W        #save peak sigma
     1136        gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma
     1137        return sig,gam
     1138               
    11131139    hId = Histogram['hId']
    11141140    hfx = ':%d:'%(hId)
    … …  
    11181144       
    11191145    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.']
     1146        shl = max(parmDict[hfx+'SH/L'],0.0005)
    11281147        Ka2 = False
    11291148        if hfx+'Lam1' in parmDict.keys():
    … …  
    11501169        for refl in refList:
    11511170            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
     1171                h,k,l = refl[:3]
     1172                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
     1173                refl[6:8] = GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse)    #peak sig & gam
     1174                Icorr = GetIntensityCorr(refl,phfx,hfx,calcControls,parmDict)
    11661175                if 'Pawley' in Phase['General']['Type']:
    11671176                    try:
    1168                         refl[8] = parmDict[hfx+'Scale']*m*parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]
     1177                        refl[8] = parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]
    11691178                    except KeyError:
    11701179                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
    … …  
    11721181                else:
    11731182                    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)
     1183                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
     1184                iBeg = np.searchsorted(x,refl[5]-fmin)
     1185                iFin = np.searchsorted(x,refl[5]+fmax)
    11791186                if not iBeg+iFin:       #peak below low limit - skip peak
    11801187                    continue
    11811188                elif not iBeg-iFin:     #peak above high limit - done
    11821189                    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
     1190                yc[iBeg:iFin] += G2pwd.getFCJVoigt(refl[5],Icorr*refl[8],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
    11841191                if Ka2:
    1185                     pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
     1192                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
    11861193                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
    1187                     if pos > 90:
    1188                         fmin,fmax = [fmax,fmin]         
    11891194                    iBeg = np.searchsorted(x,pos2-fmin)
    11901195                    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
     1196                    yc[iBeg:iFin] += G2pwd.getFCJVoigt(pos2,Icorr*refl[8]*kRatio,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
    11921197            else:
    11931198                raise ValueError
    … …  
    12271232                Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
    12281233                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])
    12301234        Histograms['sumwYo'] = sumwYo
    12311235        Histograms['Nobs'] = Nobs
    … …  
    12771281            result = so.leastsq(errRefine,values,full_output=True,  #factor=1.,epsfcn=0.00001,ftol=0.0001,
    12781282                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
    1279             print len(result)
    12801283        finally:
    12811284            dlg.Destroy()
Note: See TracChangeset for help on using the changeset viewer.