Changeset 537


Ignore:
Timestamp:
Apr 13, 2012 11:34:34 AM (11 years ago)
Author:
vondreele
Message:

some refactoring/rearrangement in GSASIIgrid.py
fixes to powderFxyeSave
move Fourier & map search math to GSASIImath.py
implement plot of map peaks & their move to atom list
begin charge flipping GUI

Location:
trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASII.py

    r536 r537  
    446446        self.Image = 0
    447447        self.oldImagefile = ''
     448        self.ImageZ = []
    448449        self.Integrate = 0
    449450        self.Pwdr = False

  • trunk/GSASIIIO.py

    r526 r537  
    490490        All files (*.*)|*.*',wx.OPEN|wx.CHANGE_DIR)
    491491        try:
    492             dlg.SetFilename(ospath.split(imagefile)[1])
     492            dlg.SetFilename(''+ospath.split(imagefile)[1])
    493493            if dlg.ShowModal() == wx.ID_OK:
    494494                imagefile = dlg.GetPath()
    … …  
    916916    head,tail = ospath.split(powderfile)
    917917    name,ext = tail.split('.')
    918     wx.BeginBusyCursor()
    919918    for i,export in enumerate(exports):
    920919        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
    … …  
    943942        file = open(filename,'w')
    944943        print 'save powder pattern to file: ',filename
    945         try:
    946             x,y,w,yc,yb,yd = G2frame.PatternTree.GetItemPyData(PickId)[1]
    947             file.write(powderfile+'\n')
    948             file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
    949                 100.*x[0],100.*(x[1]-x[0])))
    950             s = list(np.sqrt(1./np.array(w)))       
    951             XYW = zip(x,y,s)
    952             for X,Y,S in XYW:
    953                 file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,max(S,1.0)))
    954             file.close()
    955         finally:
    956             wx.EndBusyCursor()
    957         print 'powder pattern file written'
     944        x,y,w,yc,yb,yd = G2frame.PatternTree.GetItemPyData(PickId)[1]
     945        file.write(powderfile+'\n')
     946        file.write('Instrument parameter file:'+ospath.split(prmname)[1]+'\n')
     947        file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
     948            100.*x[0],100.*(x[1]-x[0])))
     949        s = list(np.sqrt(1./np.array(w)))       
     950        XYW = zip(x,y,s)
     951        for X,Y,S in XYW:
     952            file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,max(S,1.0)))
     953        file.close()
     954        print 'powder pattern file '+filename+' written'
    958955       
    959956def powderXyeSave(G2frame,exports,powderfile):
    … …  
    966963        file.write('#%s\n'%(export))
    967964        print 'save powder pattern to file: ',filename
    968         wx.BeginBusyCursor()
    969         try:
    970             x,y,w,yc,yb,yd = G2frame.PatternTree.GetItemPyData(PickId)[1]
    971             s = list(np.sqrt(1./np.array(w)))       
    972             XYW = zip(x,y,s)
    973             for X,Y,W in XYW:
    974                 file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
    975             file.close()
    976         finally:
    977             wx.EndBusyCursor()
    978         print 'powder pattern file written'
     965        x,y,w,yc,yb,yd = G2frame.PatternTree.GetItemPyData(PickId)[1]
     966        s = list(np.sqrt(1./np.array(w)))       
     967        XYW = zip(x,y,s)
     968        for X,Y,W in XYW:
     969            file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
     970        file.close()
     971        print 'powder pattern file '+filename+' written'
    979972       
    980973def PDFSave(G2frame,exports):   

  • trunk/GSASIIgrid.py

    r535 r537  
    3838htmlFirstUse = True
    3939
    40 [ wxID_FOURCALC,wxID_FOURSEARCH,
    41 ] = [wx.NewId() for item in range(2)]
     40[ wxID_FOURCALC,wxID_FOURSEARCH, wxID_PEAKSMOVE, wxID_PEAKSCLEAR,
     41] = [wx.NewId() for item in range(4)]
    4242
    4343[ wxID_PWDRADD, wxID_HKLFADD, wxID_DATADELETE,
    … …  
    448448            text='Search map')
    449449       
     450# Phase / Data tab
     451        self.DataMenu = wx.MenuBar()
     452        self.DataEdit = wx.Menu(title='')
     453        self.DataMenu.Append(menu=self.DataEdit, title='Edit')
     454        self.DataMenu.Append(menu=MyHelp(self,helpType='Data'),title='&Help')
     455        self.DataEdit.Append(id=wxID_PWDRADD, kind=wx.ITEM_NORMAL,text='Add powder histograms',
     456            help='Select new powder histograms to be used for this phase')
     457        self.DataEdit.Append(id=wxID_HKLFADD, kind=wx.ITEM_NORMAL,text='Add single crystal histograms',
     458            help='Select new single crystal histograms to be used for this phase')
     459        self.DataEdit.Append(id=wxID_DATADELETE, kind=wx.ITEM_NORMAL,text='Delete histograms',
     460            help='Delete histograms from use for this phase')
     461           
    450462# Phase / Atoms tab
    451463        self.AtomsMenu = wx.MenuBar()
    … …  
    476488            help='Compute distances & angles for selected atoms')   
    477489                 
    478 # Phase / Data tab
    479         self.DataMenu = wx.MenuBar()
    480         self.DataEdit = wx.Menu(title='')
    481         self.DataMenu.Append(menu=self.DataEdit, title='Edit')
    482         self.DataMenu.Append(menu=MyHelp(self,helpType='Data'),title='&Help')
    483         self.DataEdit.Append(id=wxID_PWDRADD, kind=wx.ITEM_NORMAL,text='Add powder histograms',
    484             help='Select new powder histograms to be used for this phase')
    485         self.DataEdit.Append(id=wxID_HKLFADD, kind=wx.ITEM_NORMAL,text='Add single crystal histograms',
    486             help='Select new single crystal histograms to be used for this phase')
    487         self.DataEdit.Append(id=wxID_DATADELETE, kind=wx.ITEM_NORMAL,text='Delete histograms',
    488             help='Delete histograms from use for this phase')
    489            
    490 # Phase / Texture tab
    491         self.TextureMenu = wx.MenuBar()
    492         self.TextureEdit = wx.Menu(title='')
    493         self.TextureMenu.Append(menu=self.TextureEdit, title='Texture')
    494         self.TextureMenu.Append(menu=MyHelp(self,helpType='Texture'),title='&Help')
    495         self.TextureEdit.Append(id=wxID_REFINETEXTURE, kind=wx.ITEM_NORMAL,text='Refine texture',
    496             help='Refine the texture coefficients from sequential Pawley results')
    497         self.TextureEdit.Append(id=wxID_CLEARTEXTURE, kind=wx.ITEM_NORMAL,text='Clear texture',
    498             help='Clear the texture coefficients' )
    499            
    500 # Phase / Pawley tab
    501         self.PawleyMenu = wx.MenuBar()
    502         self.PawleyEdit = wx.Menu(title='')
    503         self.PawleyMenu.Append(menu=self.PawleyEdit,title='Operations')
    504         self.PawleyMenu.Append(menu=MyHelp(self,helpType='Pawley'),title='&Help')
    505         self.PawleyEdit.Append(id=wxID_PAWLEYLOAD, kind=wx.ITEM_NORMAL,text='Pawley create',
    506             help='Initialize Pawley reflection list')
    507         self.PawleyEdit.Append(id=wxID_PAWLEYESTIMATE, kind=wx.ITEM_NORMAL,text='Pawley estimate',
    508             help='Estimate initial Pawley intensities')
    509         self.PawleyEdit.Append(id=wxID_PAWLEYDELETE, kind=wx.ITEM_NORMAL,text='Pawley delete',
    510             help='Delete Pawley reflection list')
    511 
    512490# Phase / Draw Options tab
    513491        self.DataDrawOptions = wx.MenuBar()
    … …  
    545523        self.DrawAtomCompute.Append(id=wxID_DRAWPLANE, kind=wx.ITEM_NORMAL,text='Best plane',
    546524            help='Compute best plane for 4+ selected atoms')   
     525           
     526# Phase / Texture tab
     527        self.TextureMenu = wx.MenuBar()
     528        self.TextureEdit = wx.Menu(title='')
     529        self.TextureMenu.Append(menu=self.TextureEdit, title='Texture')
     530        self.TextureMenu.Append(menu=MyHelp(self,helpType='Texture'),title='&Help')
     531        self.TextureEdit.Append(id=wxID_REFINETEXTURE, kind=wx.ITEM_NORMAL,text='Refine texture',
     532            help='Refine the texture coefficients from sequential Pawley results')
     533        self.TextureEdit.Append(id=wxID_CLEARTEXTURE, kind=wx.ITEM_NORMAL,text='Clear texture',
     534            help='Clear the texture coefficients' )
     535           
     536# Phase / Pawley tab
     537        self.PawleyMenu = wx.MenuBar()
     538        self.PawleyEdit = wx.Menu(title='')
     539        self.PawleyMenu.Append(menu=self.PawleyEdit,title='Operations')
     540        self.PawleyMenu.Append(menu=MyHelp(self,helpType='Pawley'),title='&Help')
     541        self.PawleyEdit.Append(id=wxID_PAWLEYLOAD, kind=wx.ITEM_NORMAL,text='Pawley create',
     542            help='Initialize Pawley reflection list')
     543        self.PawleyEdit.Append(id=wxID_PAWLEYESTIMATE, kind=wx.ITEM_NORMAL,text='Pawley estimate',
     544            help='Estimate initial Pawley intensities')
     545        self.PawleyEdit.Append(id=wxID_PAWLEYDELETE, kind=wx.ITEM_NORMAL,text='Pawley delete',
     546            help='Delete Pawley reflection list')
     547           
     548# Phase / Map peaks tab
     549        self.MapPeaksMenu = wx.MenuBar()
     550        self.MapPeaksEdit = wx.Menu(title='')
     551        self.MapPeaksMenu.Append(menu=self.MapPeaksEdit, title='Map peaks')
     552        self.MapPeaksMenu.Append(menu=MyHelp(self,helpType='Map peaks'),title='&Help')
     553        self.MapPeaksEdit.Append(id=wxID_PEAKSMOVE, kind=wx.ITEM_NORMAL,text='Move peaks',
     554            help='Move selected peaks to atom list')
     555        self.MapPeaksEdit.Append(id=wxID_PEAKSCLEAR, kind=wx.ITEM_NORMAL,text='Clear peaks',
     556            help='Clear the map peak list')
    547557           
    548558# end of GSAS-II menu definitions

  • trunk/GSASIIimage.py

    r536 r537  
    673673    print 'Begin image integration'
    674674    LUtth = data['IOtth']
    675     LRazm = np.array(data['LRazimuth'])
     675    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
    676676    numAzms = data['outAzimuths']
    677677    numChans = data['outChannels']

  • trunk/GSASIImath.py

    r530 r537  
    1212import numpy as np
    1313import numpy.linalg as nl
     14import numpy.ma as ma
    1415import cPickle
    1516import time
    1617import math
     18import copy
    1719import GSASIIpath
     20import GSASIIlattice as G2lat
    1821import GSASIIspc as G2spc
    1922import scipy.optimize as so
    … …  
    514517            return text
    515518
    516    
     519def FourierMap(data,reflData):
     520   
     521    import scipy.fftpack as fft
     522    generalData = data['General']
     523    if not generalData['Map']['MapType']:
     524        print '**** ERROR - Fourier map not defined'
     525        return
     526    mapData = generalData['Map']
     527    dmin = mapData['Resolution']
     528    SGData = generalData['SGData']
     529    cell = generalData['Cell'][1:8]       
     530    A = G2lat.cell2A(cell[:6])
     531    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
     532    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
     533#    Fhkl[0,0,0] = generalData['F000X']
     534    time0 = time.time()
     535    for ref in reflData:
     536        if ref[4] >= dmin:
     537            for i,hkl in enumerate(ref[11]):
     538                hkl = np.asarray(hkl,dtype='i')
     539                Fosq,Fcsq,ph = ref[8:11]
     540                dp = 360.*ref[12][i]
     541                a = cosd(ph+dp)
     542                b = sind(ph+dp)
     543                phasep = complex(a,b)
     544                phasem = complex(a,-b)
     545                if 'Fobs' in mapData['MapType']:
     546                    F = np.sqrt(Fosq)
     547                    h,k,l = hkl+Hmax
     548                    Fhkl[h,k,l] = F*phasep
     549                    h,k,l = -hkl+Hmax
     550                    Fhkl[h,k,l] = F*phasem
     551                elif 'Fcalc' in mapData['MapType']:
     552                    F = np.sqrt(Fcsq)
     553                    h,k,l = hkl+Hmax
     554                    Fhkl[h,k,l] = F*phasep
     555                    h,k,l = -hkl+Hmax
     556                    Fhkl[h,k,l] = F*phasem
     557                elif 'delt-F' in mapData['MapType']:
     558                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
     559                    h,k,l = hkl+Hmax
     560                    Fhkl[h,k,l] = dF*phasep
     561                    h,k,l = -hkl+Hmax
     562                    Fhkl[h,k,l] = dF*phasem
     563                elif '2*Fo-Fc' in mapData['MapType']:
     564                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
     565                    h,k,l = hkl+Hmax
     566                    Fhkl[h,k,l] = F*phasep
     567                    h,k,l = -hkl+Hmax
     568                    Fhkl[h,k,l] = F*phasem
     569                elif 'Patterson' in mapData['MapType']:
     570                    h,k,l = hkl+Hmax
     571                    Fhkl[h,k,l] = complex(Fosq,0.)
     572                    h,k,l = -hkl+Hmax
     573                    Fhkl[h,k,l] = complex(Fosq,0.)
     574    Fhkl = fft.fftshift(Fhkl)
     575    rho = fft.fftn(Fhkl)/cell[6]
     576    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
     577    mapData['rho'] = np.real(rho)
     578    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
     579    return mapData
     580   
     581def SearchMap(data,keepDup=False):
     582   
     583    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
     584   
     585    def noDuplicate(xyz,peaks,SGData):                  #be careful where this is used - it's slow
     586        equivs = G2spc.GenAtom(xyz,SGData)
     587        xyzs = [equiv[0] for equiv in equivs]
     588        for x in xyzs:
     589            if True in [np.allclose(x,peak,atol=0.002) for peak in peaks]:
     590                return False
     591        return True
     592           
     593    def findRoll(rhoMask,mapHalf):
     594        indx = np.array(ma.nonzero(rhoMask)).T
     595        rhoList = np.array([rho[i,j,k] for i,j,k in indx])
     596        rhoMax = np.max(rhoList)
     597        return mapHalf-indx[np.argmax(rhoList)]
     598       
     599    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
     600        Mag,x0,y0,z0,sig = parms
     601        if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
     602            return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(x0-rX)*(y0-rY)+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
     603        else:
     604            return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
     605       
     606    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
     607        Mag,x0,y0,z0,sig = parms
     608        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
     609        return M
     610       
     611    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
     612        Mag,x0,y0,z0,sig = parms
     613        dMdv = np.zeros(([5,]+list(rX.shape)))
     614        delt = .01
     615        for i in range(5):
     616            parms[i] -= delt
     617            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
     618            parms[i] += 2.*delt
     619            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
     620            parms[i] -= delt
     621            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
     622        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
     623        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
     624        dMdv = np.reshape(dMdv,(5,rX.size))
     625        Hess = np.inner(dMdv,dMdv)
     626       
     627        return Vec,Hess
     628       
     629    generalData = data['General']
     630    phaseName = generalData['Name']
     631    SGData = generalData['SGData']
     632    cell = generalData['Cell'][1:8]       
     633    A = G2lat.cell2A(cell[:6])
     634    drawingData = data['Drawing']
     635    peaks = []
     636    mags = []
     637    try:
     638        mapData = generalData['Map']
     639        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
     640        rho = copy.copy(mapData['rho'])     #don't mess up original
     641        mapHalf = np.array(rho.shape)/2
     642        res = mapData['Resolution']
     643        incre = 1./np.array(rho.shape)
     644        step = max(1.0,1./res)+1
     645        steps = np.array(3*[step,])
     646    except KeyError:
     647        print '**** ERROR - Fourier map not defined'
     648        return peaks,mags
     649    while True:
     650        rhoMask = ma.array(rho,mask=(rho<contLevel))
     651        if not ma.count(rhoMask):
     652            break
     653        rMI = findRoll(rhoMask,mapHalf)
     654        rho = np.roll(np.roll(np.roll(rho,rMI[0],axis=0),rMI[1],axis=1),rMI[2],axis=2)
     655        rMM = mapHalf-steps
     656        rMP = mapHalf+steps+1
     657        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
     658        peakInt = np.sum(rhoPeak)*res**3
     659        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
     660        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
     661        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
     662            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.0001,maxcyc=100)
     663        x1 = result[0]
     664        if np.any(x1 < 0):
     665            break
     666        peak = (np.array(x1[1:4])-rMI)*incre
     667        if not len(peaks):
     668            peaks.append(peak)
     669            mags.append(x1[0])
     670        else:
     671            if keepDup:
     672                peaks.append(peak)
     673                mags.append(x1[0])
     674            elif noDuplicate(peak,peaks,SGData) and x1[0] > 0.:
     675                peaks.append(peak)
     676                mags.append(x1[0])
     677        rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] = peakFunc(result[0],rX,rY,rZ,rhoPeak,res,SGData['SGLaue'])
     678        rho = np.roll(np.roll(np.roll(rho,-rMI[2],axis=2),-rMI[1],axis=1),-rMI[0],axis=0)
     679    return np.array(peaks),np.array([mags,]).T

  • trunk/GSASIIphsGUI.py

    r530 r537  
    255255        if 'Map' not in generalData:
    256256            generalData['Map'] = {'MapType':'','RefList':'','Resolution':1.0,
    257                 'rhoMax':100.,'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.}
     257                'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.}
     258        if 'Flip' not in generalData:
     259            generalData['Flip'] = {'RefList':'','Resolution':1.0,'Norm element':'C',
     260                'k-factor':1.1,'rhoSig':0.0,'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.} #???
     261           
    258262#        if 'SH Texture' not in generalData:
    259263#            generalData['SH Texture'] = data['SH Texture']
    … …  
    674678            return mapSizer
    675679               
     680        def FlipSizer():
     681           
     682            def OnMapType(event):
     683                Map['MapType'] = mapType.GetValue()
     684               
     685            def OnRefList(event):
     686                Flip['RefList'] = refList.GetValue()
     687               
     688            def OnResVal(event):
     689                try:
     690                    res = float(mapRes.GetValue())
     691                    if 0.25 <= res <= 20.:
     692                        Flip['Resolution'] = res
     693                except ValueError:
     694                    pass
     695                flipRes.SetValue("%.2f"%(Flip['Resolution']))          #reset in case of error
     696           
     697            def OnCutOff(event):
     698                try:
     699                    res = float(cutOff.GetValue())
     700                    if 1.0 <= res <= 100.:
     701                        Map['cutOff'] = res
     702                except ValueError:
     703                    pass
     704                cutOff.SetValue("%.1f"%(Map['cutOff']))          #reset in case of error
     705           
     706            refList = data['Histograms'].keys()
     707            flipSizer = wx.BoxSizer(wx.VERTICAL)
     708            lineSizer = wx.BoxSizer(wx.HORIZONTAL)
     709            lineSizer.Add(wx.StaticText(dataDisplay,label=' Charge flip controls: Reflection set from: '),0,wx.ALIGN_CENTER_VERTICAL)
     710            refList = wx.ComboBox(dataDisplay,-1,value=Map['RefList'],choices=refList,
     711                style=wx.CB_READONLY|wx.CB_DROPDOWN)
     712            refList.Bind(wx.EVT_COMBOBOX,OnRefList)
     713            lineSizer.Add(refList,0,wx.ALIGN_CENTER_VERTICAL)
     714            flipSizer.Add(lineSizer,0,wx.ALIGN_CENTER_VERTICAL)
     715            line2Sizer = wx.BoxSizer(wx.HORIZONTAL)
     716            line2Sizer.Add(wx.StaticText(dataDisplay,label=' Resolution: '),0,wx.ALIGN_CENTER_VERTICAL)
     717            mapRes =  wx.TextCtrl(dataDisplay,value='%.2f'%(Map['Resolution']),style=wx.TE_PROCESS_ENTER)
     718            mapRes.Bind(wx.EVT_TEXT_ENTER,OnResVal)       
     719            mapRes.Bind(wx.EVT_KILL_FOCUS,OnResVal)
     720            line2Sizer.Add(mapRes,0,wx.ALIGN_CENTER_VERTICAL)
     721            line2Sizer.Add(wx.StaticText(dataDisplay,label=' Peak cutoff %: '),0,wx.ALIGN_CENTER_VERTICAL)
     722            cutOff =  wx.TextCtrl(dataDisplay,value='%.1f'%(Map['cutOff']),style=wx.TE_PROCESS_ENTER)
     723            cutOff.Bind(wx.EVT_TEXT_ENTER,OnCutOff)       
     724            cutOff.Bind(wx.EVT_KILL_FOCUS,OnCutOff)
     725            line2Sizer.Add(cutOff,0,wx.ALIGN_CENTER_VERTICAL)
     726            flipSizer.Add(line2Sizer,0,wx.ALIGN_CENTER_VERTICAL)
     727            return flipSizer
     728               
    676729        General.DestroyChildren()
    677730        dataDisplay = wx.Panel(General)
    … …  
    694747        mainSizer.Add((5,5),0)
    695748        mainSizer.Add(MapSizer())
     749
     750        mainSizer.Add((5,5),0)
     751        mainSizer.Add(FlipSizer())
    696752
    697753        dataDisplay.SetSizer(mainSizer)
    … …  
    33723428
    33733429################################################################################
    3374 ##### End of main routines
     3430##### Fourier routines
    33753431################################################################################
     3432
     3433    def FillMapPeaksGrid():
     3434                       
     3435        def RowSelect(event):
     3436            r,c =  event.GetRow(),event.GetCol()
     3437            if r < 0 and c < 0:
     3438                if MapPeaks.IsSelection():
     3439                    MapPeaks.ClearSelection()
     3440                else:
     3441                    for row in range(MapPeaks.GetNumberRows()):
     3442                        MapPeaks.SelectRow(row,True)
     3443                   
     3444            elif c < 0:                   #only row clicks
     3445                if event.ControlDown():                   
     3446                    if r in MapPeaks.GetSelectedRows():
     3447                        MapPeaks.DeselectRow(r)
     3448                    else:
     3449                        MapPeaks.SelectRow(r,True)
     3450                elif event.ShiftDown():
     3451                    for row in range(r+1):
     3452                        MapPeaks.SelectRow(row,True)
     3453                else:
     3454                    MapPeaks.ClearSelection()
     3455                    MapPeaks.SelectRow(r,True)               
     3456            G2plt.PlotStructure(G2frame,data)                   
     3457           
     3458        G2frame.dataFrame.setSizePosLeft([450,300])
     3459        if 'Map Peaks' in data:
     3460            mapPeaks = data['Map Peaks']                       
     3461            rowLabels = []
     3462            for i in range(len(mapPeaks)): rowLabels.append(str(i))
     3463            colLabels = ['mag','x','y','z']
     3464            Types = 4*[wg.GRID_VALUE_FLOAT+':10,4',]
     3465            MapPeaksTable = G2gd.Table(mapPeaks,rowLabels=rowLabels,colLabels=colLabels,types=Types)
     3466            MapPeaks.SetTable(MapPeaksTable, True)
     3467            MapPeaks.Bind(wg.EVT_GRID_LABEL_LEFT_CLICK, RowSelect)
     3468            for r in range(MapPeaks.GetNumberRows()):
     3469                for c in range(MapPeaks.GetNumberCols()):
     3470                    MapPeaks.SetCellStyle(r,c,VERY_LIGHT_GREY,True)
     3471            MapPeaks.SetMargins(0,0)
     3472            MapPeaks.AutoSizeColumns(False)
     3473                   
     3474    def OnPeaksMove(event):
     3475        if 'Map Peaks' in data:
     3476            mapPeaks = data['Map Peaks']                       
     3477            Ind = MapPeaks.GetSelectedRows()
     3478            for ind in Ind:
     3479                print mapPeaks[ind]
     3480                x,y,z = mapPeaks[ind][1:]
     3481                AtomAdd(x,y,z)
     3482   
     3483    def OnPeaksClear(event):
     3484        data['Map Peaks'] = []
     3485        FillMapPeaksGrid()
     3486        G2plt.PlotStructure(G2frame,data)
    33763487   
    33773488    def OnFourierMaps(event):
    3378         import scipy.fftpack as fft
    33793489        generalData = data['General']
    3380         if not generalData['Map']['MapType']:
    3381             print '**** ERROR - Fourier map not defined'
    3382             return
    33833490        mapData = generalData['Map']
    3384         dmin = mapData['Resolution']
     3491        reflName = mapData['RefList']
    33853492        phaseName = generalData['Name']
    3386         SGData = generalData['SGData']
    3387         cell = generalData['Cell'][1:8]       
    3388         A = G2lat.cell2A(cell[:6])
    3389         Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
    3390         Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
    3391         reflName = mapData['RefList']
    33923493        if 'PWDR' in reflName:
    33933494            PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, reflName)
    … …  
    33973498            print 'single crystal reflections'
    33983499            return
    3399 #        Fhkl[0,0,0] = generalData['F000X']
    3400         time0 = time.time()
    3401         for ref in reflData:
    3402             if ref[4] >= dmin:
    3403                 for i,hkl in enumerate(ref[11]):
    3404                     hkl = np.asarray(hkl,dtype='i')
    3405                     Fosq,Fcsq,ph = ref[8:11]
    3406                     dp = 360.*ref[12][i]
    3407                     a = cosd(ph+dp)
    3408                     b = sind(ph+dp)
    3409                     phasep = complex(a,b)
    3410                     phasem = complex(a,-b)
    3411                     if 'Fobs' in mapData['MapType']:
    3412                         F = np.sqrt(Fosq)
    3413                         h,k,l = hkl+Hmax
    3414                         Fhkl[h,k,l] = F*phasep
    3415                         h,k,l = -hkl+Hmax
    3416                         Fhkl[h,k,l] = F*phasem
    3417                     elif 'Fcalc' in mapData['MapType']:
    3418                         F = np.sqrt(Fcsq)
    3419                         h,k,l = hkl+Hmax
    3420                         Fhkl[h,k,l] = F*phasep
    3421                         h,k,l = -hkl+Hmax
    3422                         Fhkl[h,k,l] = F*phasem
    3423                     elif 'delt-F' in mapData['MapType']:
    3424                         dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
    3425                         h,k,l = hkl+Hmax
    3426                         Fhkl[h,k,l] = dF*phasep
    3427                         h,k,l = -hkl+Hmax
    3428                         Fhkl[h,k,l] = dF*phasem
    3429                     elif '2*Fo-Fc' in mapData['MapType']:
    3430                         F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
    3431                         h,k,l = hkl+Hmax
    3432                         Fhkl[h,k,l] = F*phasep
    3433                         h,k,l = -hkl+Hmax
    3434                         Fhkl[h,k,l] = F*phasem
    3435                     elif 'Patterson' in mapData['MapType']:
    3436                         h,k,l = hkl+Hmax
    3437                         Fhkl[h,k,l] = complex(Fosq,0.)
    3438                         h,k,l = -hkl+Hmax
    3439                         Fhkl[h,k,l] = complex(Fosq,0.)
    3440         Fhkl = fft.fftshift(Fhkl)
    3441         rho = fft.fftn(Fhkl)/cell[6]
    3442         print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
    3443         mapData['rho'] = np.real(rho)
    3444         mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
     3500        mapData.update(G2mth.FourierMap(data,reflData))
     3501        mapSig = np.std(mapData['rho'])
    34453502        data['Drawing']['contourLevel'] = 1.
    34463503        data['Drawing']['mapSize'] = 10.
    3447         print mapData['MapType']+' computed: rhomax = %.3f rhomin = %.3f'%(np.max(mapData['rho']),np.min(mapData['rho']))
     3504        print mapData['MapType']+' computed: rhomax = %.3f rhomin = %.3f sigma = %.3f'%(np.max(mapData['rho']),np.min(mapData['rho']),mapSig)
    34483505        G2plt.PlotStructure(G2frame,data)
    34493506       
    … …  
    34653522    def OnSearchMaps(event):
    34663523       
    3467         norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
    3468        
    3469         def noDuplicate(xyz,peaks,SGData):                  #be careful where this is used - it's slow
    3470             equivs = G2spc.GenAtom(xyz,SGData)
    3471             xyzs = [equiv[0] for equiv in equivs]
    3472             for x in xyzs:
    3473                 if True in [np.allclose(x,peak,atol=0.002) for peak in peaks]:
    3474                     return False
    3475             return True
    3476                
    3477         def findRoll(rhoMask,mapHalf):
    3478             indx = np.array(ma.nonzero(rhoMask)).T
    3479             rhoList = np.array([rho[i,j,k] for i,j,k in indx])
    3480             rhoMax = np.max(rhoList)
    3481             return mapHalf-indx[np.argmax(rhoList)]
    3482            
    3483         def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
    3484             Mag,x0,y0,z0,sig = parms
    3485             if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
    3486                 return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(x0-rX)*(y0-rY)+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
    3487             else:
    3488                 return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
    3489            
    3490         def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
    3491             Mag,x0,y0,z0,sig = parms
    3492             M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
    3493             return M
    3494            
    3495         def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
    3496             Mag,x0,y0,z0,sig = parms
    3497             dMdv = np.zeros(([5,]+list(rX.shape)))
    3498             delt = .01
    3499             for i in range(5):
    3500                 parms[i] -= delt
    3501                 rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
    3502                 parms[i] += 2.*delt
    3503                 rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
    3504                 parms[i] -= delt
    3505                 dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
    3506             rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
    3507             Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
    3508             dMdv = np.reshape(dMdv,(5,rX.size))
    3509             Hess = np.inner(dMdv,dMdv)
    3510            
    3511             return Vec,Hess
    3512            
    3513         generalData = data['General']
    3514         phaseName = generalData['Name']
    3515         SGData = generalData['SGData']
    3516         cell = generalData['Cell'][1:8]       
    3517         A = G2lat.cell2A(cell[:6])
    3518         drawingData = data['Drawing']
    3519         try:
    3520             mapData = generalData['Map']
    3521             contLevel = mapData['cutOff']*mapData['rhoMax']/100.
    3522             rho = copy.copy(mapData['rho'])     #don't mess up original
    3523             mapHalf = np.array(rho.shape)/2
    3524             res = mapData['Resolution']
    3525             incre = 1./np.array(rho.shape)
    3526             step = max(1.0,1./res)+1
    3527             steps = np.array(3*[step,])
    3528         except KeyError:
    3529             print '**** ERROR - Fourier map not defined'
    3530             return
    35313524        peaks = []
    35323525        mags = []
    … …  
    35353528        wx.BeginBusyCursor()
    35363529        try:
    3537             while True:
    3538                 rhoMask = ma.array(rho,mask=(rho<contLevel))
    3539                 if not ma.count(rhoMask):
    3540                     break
    3541                 rMI = findRoll(rhoMask,mapHalf)
    3542                 rho = np.roll(np.roll(np.roll(rho,rMI[0],axis=0),rMI[1],axis=1),rMI[2],axis=2)
    3543                 rMM = mapHalf-steps
    3544                 rMP = mapHalf+steps+1
    3545                 rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
    3546                 peakInt = np.sum(rhoPeak)*res**3
    3547                 rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
    3548                 x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
    3549                 result = G2mth.HessianLSQ(peakFunc,x0,Hess=peakHess,
    3550                     args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.0001,maxcyc=100)
    3551                 x1 = result[0]
    3552                 if np.any(x1 < 0):
    3553                     break
    3554                 peak = (np.array(x1[1:4])-rMI)*incre
    3555                 if not len(peaks):
    3556                     peaks.append(peak)
    3557                     mags.append(x1[0])
    3558                 else:
    3559                     if noDuplicate(peak,peaks,SGData) and x1[0] > 0.:
    3560                         peaks.append(peak)
    3561                         mags.append(x1[0])
    3562                 rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] = peakFunc(result[0],rX,rY,rZ,rhoPeak,res,SGData['SGLaue'])
    3563                 rho = np.roll(np.roll(np.roll(rho,-rMI[2],axis=2),-rMI[1],axis=1),-rMI[0],axis=0)
     3530            peaks,mags = G2mth.SearchMap(data,keepDup=True)
    35643531        finally:
    35653532            wx.EndBusyCursor()
    3566         sortIdx = np.argsort(mags)       
    3567         print ' Map search peaks found:'
    3568         print '  No.    Mag.      x        y        z'
    3569         for j,i in enumerate(sortIdx[::-1]):
    3570             print ' %3d %8.3f %8.5f %8.5f %8.5f'%(j,mags[i],peaks[i][0],peaks[i][1],peaks[i][2])
    3571         print ' Map search finished, time = %.2fs'%(time.time()-time0)
    3572        
     3533        sortIdx = np.argsort(mags.flatten())
     3534        if len(peaks):
     3535            data['Map Peaks'] = np.concatenate((mags,peaks),axis=1)
     3536           
     3537            print ' Map search peaks found:'
     3538            print '  No.    Mag.      x        y        z'
     3539            for j,i in enumerate(sortIdx[::-1]):
     3540                print ' %3d %8.3f %8.5f %8.5f %8.5f'%(j,mags[i],peaks[i][0],peaks[i][1],peaks[i][2])
     3541            print ' Map search finished, time = %.2fs'%(time.time()-time0)
     3542        Page = G2frame.dataDisplay.FindPage('Map peaks')
     3543        G2frame.dataDisplay.ChangeSelection(Page)
     3544        G2frame.dataFrame.SetMenuBar(G2frame.dataFrame.MapPeaksMenu)
     3545        G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksMove, id=G2gd.wxID_PEAKSMOVE)
     3546        G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksClear, id=G2gd.wxID_PEAKSCLEAR)
     3547        FillMapPeaksGrid()
     3548        G2plt.PlotStructure(G2frame,data)
     3549               
    35733550    def OnTextureRefine(event):
    35743551        print 'refine texture?'
    … …  
    36393616            UpdateTexture()                       
    36403617            G2plt.PlotTexture(G2frame,data,Start=True)
     3618        elif text == 'Map peaks':
     3619            G2frame.dataFrame.SetMenuBar(G2frame.dataFrame.MapPeaksMenu)
     3620            G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksMove, id=G2gd.wxID_PEAKSMOVE)
     3621            G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksClear, id=G2gd.wxID_PEAKSCLEAR)
     3622            FillMapPeaksGrid()
     3623            G2plt.PlotStructure(G2frame,data)
     3624           
    36413625        else:
    36423626            G2frame.dataFrame.SetMenuBar(G2frame.dataFrame.BlankMenu)
    … …  
    36703654        Texture = wx.ScrolledWindow(G2frame.dataDisplay)
    36713655        G2frame.dataDisplay.AddPage(Texture,'Texture')
    3672 
     3656        MapPeaks = G2gd.GSGrid(G2frame.dataDisplay)
     3657        G2frame.dataDisplay.AddPage(MapPeaks,'Map peaks')
     3658           
    36733659    G2frame.dataDisplay.Bind(wx.EVT_NOTEBOOK_PAGE_CHANGED, OnPageChanged)
    36743660    G2frame.dataDisplay.SetSelection(oldPage)

  • trunk/GSASIIplot.py

    r536 r537  
    17071707        Page.canvas.SetToolTipString('')
    17081708        sizexy = Data['size']
    1709         if event.xdata and event.ydata:                 #avoid out of frame errors
     1709        if event.xdata and event.ydata and len(G2frame.ImageZ):                 #avoid out of frame errors
    17101710            Page.canvas.SetCursor(wx.CROSS_CURSOR)
    17111711            item = G2frame.itemPicked
    … …  
    18261826        scaley = 1000./pixelSize[1]
    18271827        pixLimit = Data['pixLimit']
    1828         if G2frame.itemPicked is None and PickName == 'Image Controls':
     1828        if G2frame.itemPicked is None and PickName == 'Image Controls' and len(G2frame.ImageZ):
    18291829#            sizexy = Data['size']
    18301830            Xpos = event.xdata
    … …  
    19861986        G2frame.PatternTree.SetItemPyData(G2frame.Image,[size,imagefile])
    19871987        G2frame.ImageZ = G2IO.GetImageData(G2frame,imagefile,imageOnly=True)
    1988 #        print G2frame.ImageZ.shape,G2frame.ImageZ.size,Data['size'] #might be useful debugging line
    19891988        G2frame.oldImagefile = imagefile
    19901989
    … …  
    22812280    Mydir = generalData['Mydir']
    22822281    atomData = data['Atoms']
     2282    if 'Map Peaks' in data:
     2283        mapPeaks = data['Map Peaks']                       
    22832284    drawingData = data['Drawing']
    22842285    drawAtoms = drawingData['Atoms']
    … …  
    23272328        if Add:
    23282329            Indx = GetSelectedAtoms()
    2329         for i,atom in enumerate(drawAtoms):
    2330             x,y,z = atom[cx:cx+3]
    2331             X,Y,Z = gluProject(x,y,z,Model,Proj,View)
    2332             XY = [int(X),int(View[3]-Y)]
    2333             if np.allclose(xy,XY,atol=10) and Z < Zmax:
    2334                 Zmax = Z
    2335                 try:
    2336                     Indx.remove(i)
    2337                     ClearSelectedAtoms()
    2338                     for id in Indx:
    2339                         SetSelectedAtoms(id,Add)
    2340                 except:
    2341                     SetSelectedAtoms(i,Add)
    2342                    
     2330        if G2frame.dataDisplay.GetPageText(getSelection()) == 'Map peaks':
     2331            for i,peak in enumerate(mapPeaks):
     2332                x,y,z = peak[1:]
     2333                X,Y,Z = gluProject(x,y,z,Model,Proj,View)
     2334                XY = [int(X),int(View[3]-Y)]
     2335                if np.allclose(xy,XY,atol=10) and Z < Zmax:
     2336                    Zmax = Z
     2337                    try:
     2338                        Indx.remove(i)
     2339                        ClearSelectedAtoms()
     2340                        for id in Indx:
     2341                            SetSelectedAtoms(id,Add)
     2342                    except:
     2343                        SetSelectedAtoms(i,Add)
     2344        else:
     2345            for i,atom in enumerate(drawAtoms):
     2346                x,y,z = atom[cx:cx+3]
     2347                X,Y,Z = gluProject(x,y,z,Model,Proj,View)
     2348                XY = [int(X),int(View[3]-Y)]
     2349                if np.allclose(xy,XY,atol=10) and Z < Zmax:
     2350                    Zmax = Z
     2351                    try:
     2352                        Indx.remove(i)
     2353                        ClearSelectedAtoms()
     2354                        for id in Indx:
     2355                            SetSelectedAtoms(id,Add)
     2356                    except:
     2357                        SetSelectedAtoms(i,Add)
     2358                                       
    23432359    def OnMouseDown(event):
    23442360        xy = event.GetPosition()
    … …  
    24182434            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
    24192435                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Draw Atoms
     2436            elif G2frame.dataDisplay.GetPageText(page) == 'Map peaks':
     2437                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Atoms
    24202438            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
    24212439                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Atoms
     2440               
    24222441                   
    24232442    def SetSelectedAtoms(ind,Add=False):
    … …  
    24262445            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
    24272446                G2frame.dataDisplay.GetPage(page).SelectRow(ind,Add)      #this is the Atoms grid in Draw Atoms
     2447            elif G2frame.dataDisplay.GetPageText(page) == 'Map peaks':
     2448                G2frame.dataDisplay.GetPage(page).SelectRow(ind,Add)                 
    24282449            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
    24292450                Id = drawAtoms[ind][-2]
    … …  
    24382459            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
    24392460                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Draw Atoms
     2461            elif G2frame.dataDisplay.GetPageText(page) == 'Map peaks':
     2462                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()
    24402463            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
    24412464                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Atoms
    … …  
    26592682            Dx = np.inner(Amat,bond)
    26602683            Z = np.sqrt(np.sum(Dx**2))
    2661             azm = atan2d(-Dx[1],-Dx[0])
    2662             phi = acosd(Dx[2]/Z)
    2663             glRotate(-azm,0,0,1)
    2664             glRotate(phi,1,0,0)
    2665             q = gluNewQuadric()
    2666             gluCylinder(q,radius,radius,Z,slice,2)
     2684            if Z:
     2685                azm = atan2d(-Dx[1],-Dx[0])
     2686                phi = acosd(Dx[2]/Z)
     2687                glRotate(-azm,0,0,1)
     2688                glRotate(phi,1,0,0)
     2689                q = gluNewQuadric()
     2690                gluCylinder(q,radius,radius,Z,slice,2)
    26672691            glPopMatrix()           
    26682692        glPopMatrix()
    … …  
    26982722            glEnd()
    26992723        glPopMatrix()
     2724
     2725    def RenderMapPeak(x,y,z,color):
     2726        vecs = np.array([[[-.01,0,0],[.01,0,0]],[[0,-.01,0],[0,.01,0]],[[0,0,-.01],[0,0,.01]]])
     2727        xyz = np.array([x,y,z])
     2728        glEnable(GL_COLOR_MATERIAL)
     2729        glLineWidth(3)
     2730        glColor3fv(color)
     2731        glPushMatrix()
     2732        glBegin(GL_LINES)
     2733        for vec in vecs:
     2734            glVertex3fv(vec[0]+xyz)
     2735            glVertex3fv(vec[1]+xyz)
     2736        glEnd()
     2737        glColor4ubv([0,0,0,0])
     2738        glPopMatrix()
     2739        glDisable(GL_COLOR_MATERIAL)
    27002740       
    27012741    def RenderBackbone(Backbone,BackboneColor,radius):
    … …  
    27372777    def Draw():
    27382778        mapData = generalData['Map']
     2779        pageName = ''
     2780        page = getSelection()
     2781        if page:
     2782            pageName = G2frame.dataDisplay.GetPageText(page)
    27392783        rhoXYZ = []
    27402784        if len(mapData['rho']):
    … …  
    28152859            CL = atom[cs+2]
    28162860            color = np.array(CL)/255.
    2817             if iat in Ind:
     2861            if iat in Ind and G2frame.dataDisplay.GetPageText(getSelection()) != 'Map peaks':
    28182862                color = np.array(Gr)/255.
    28192863#            color += [.25,]
    … …  
    28932937        if len(rhoXYZ):
    28942938            RenderMap(rho,rhoXYZ,indx,Rok)
     2939        if len(mapPeaks):
     2940            for ind,[mag,x,y,z] in enumerate(mapPeaks):
     2941                if ind in Ind and pageName == 'Map peaks':
     2942                    RenderMapPeak(x,y,z,Gr)
     2943                else:
     2944                    RenderMapPeak(x,y,z,Wt)
    28952945        if Backbone:
    28962946            RenderBackbone(Backbone,BackboneColor,bondR)
Note: See TracChangeset for help on using the changeset viewer.