Changeset 520


Ignore:
Timestamp:
Mar 14, 2012 3:17:17 PM (12 years ago)
Author:
vondreele
Message:

separate getHKLmax

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIlattice.py

    r432 r520  
    590590                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
    591591    return sortHKLd(HKL,True,False)
     592
     593def getHKLmax(dmin,SGData,A):
     594    #finds maximum allowed hkl for given A within dmin
     595    SGLaue = SGData['SGLaue']
     596    if SGLaue in ['3R','3mR']:        #Rhombohedral axes
     597        Hmax = [0,0,0]
     598        cell = A2cell(A)
     599        aHx = cell[0]*math.sqrt(2.0*(1.0-cosd(cell[3])))
     600        cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3])))
     601        Hmax[0] = Hmax[1] = int(round(aHx/dmin))
     602        Hmax[2] = int(round(cHx/dmin))
     603        #print Hmax,aHx,cHx
     604    else:                           # all others
     605        Hmax = MaxIndex(dmin,A)
     606    return Hmax
    592607   
    593608def GenHLaue(dmin,SGData,A):
    … …  
    614629    SGUniq = SGData['SGUniq']
    615630    #finds maximum allowed hkl for given A within dmin
    616     if SGLaue in ['3R','3mR']:        #Rhombohedral axes
    617         Hmax = [0,0,0]
    618         cell = A2cell(A)
    619         aHx = cell[0]*math.sqrt(2.0*(1.0-cosd(cell[3])))
    620         cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3])))
    621         Hmax[0] = Hmax[1] = int(round(aHx/dmin))
    622         Hmax[2] = int(round(cHx/dmin))
    623         #print Hmax,aHx,cHx
    624     else:                           # all others
    625         Hmax = MaxIndex(dmin,A)
     631    Hmax = getHKLmax(dmin,SGData,A)
    626632       
    627633    dminsq = 1./(dmin**2)

  • trunk/GSASIIphsGUI.py

    r519 r520  
    252252            generalData['POhkl'] = [0,0,1]
    253253        if 'Map' not in generalData:
    254             generalData['Map'] = {'MapType':'','RefList':'','Resolution':4.0}
     254            generalData['Map'] = {'MapType':'','RefList':'','Resolution':1.0,'rhoMax':100.,'rho':[],'rhoMax':0.}
    255255#        if 'SH Texture' not in generalData:
    256256#            generalData['SH Texture'] = data['SH Texture']
    … …  
    274274                Info = G2elem.GetAtomInfo(atom[ct])
    275275                generalData['AtomTypes'].append(atom[ct])
     276                generalData['Z'] = Info['Z']
    276277                generalData['Isotopes'][atom[ct]] = Info['Isotopes']
    277278                generalData['BondRadii'].append(Info['Drad'])
    … …  
    285286                generalData['NoAtoms'][atom[ct]] = atom[cs-1]*float(atom[cs+1])
    286287                generalData['Color'].append(Info['Color'])
     288        F000X = 0.
     289        F000N = 0.
     290        for i,elem in enumerate(generalData['AtomTypes']):
     291            F000X += generalData['NoAtoms'][elem]*generalData['Z']
     292            isotope = generalData['Isotope'][elem]
     293            F000N += generalData['NoAtoms'][elem]*generalData['Isotopes'][elem][isotope][1]
     294        generalData['F000X'] = F000X
     295        generalData['F000N'] = F000N
     296       
    287297
    288298################################################################################
    … …  
    308318        SetupGeneral()
    309319        generalData = data['General']
    310         Map = generalData['Map']  # {'MapType':'','RefList':'','Resolution':4.0}
     320        Map = generalData['Map']  # {'MapType':'','RefList':'','Resolution':1.0,'rho':[],'rhoMax':0.}
    311321       
    312322        def NameSizer():
    … …  
    614624                except ValueError:
    615625                    pass
    616                 mapRes.SetValue("%.1f"%(Map['Resolution']))          #reset in case of error
     626                mapRes.SetValue("%.2f"%(Map['Resolution']))          #reset in case of error
    617627           
    618628            mapTypes = ['Fobs','Fcalc','delt-F','2*Fo-Fc','Patterson']
    … …  
    633643            mapSizer.Add(refList,0,wx.ALIGN_CENTER_VERTICAL)
    634644            mapSizer.Add(wx.StaticText(dataDisplay,label=' Resolution: '),0,wx.ALIGN_CENTER_VERTICAL)
    635             mapRes =  wx.TextCtrl(dataDisplay,value='%.1f'%(Map['Resolution']),style=wx.TE_PROCESS_ENTER)
     645            mapRes =  wx.TextCtrl(dataDisplay,value='%.2f'%(Map['Resolution']),style=wx.TE_PROCESS_ENTER)
    636646            mapRes.Bind(wx.EVT_TEXT_ENTER,OnResVal)       
    637647            mapRes.Bind(wx.EVT_KILL_FOCUS,OnResVal)
    … …  
    11901200            'L','K','M','F','P','S','T','W','Y','V','M',' ',' ',' ']
    11911201        defaultDrawing = {'viewPoint':[[0.5,0.5,0.5],[]],'showHydrogen':True,'backColor':[0,0,0],'depthFog':False,
    1192             'Zclip':50.0,'cameraPos':50.,'radiusFactor':0.85,
     1202            'Zclip':50.0,'cameraPos':50.,'radiusFactor':0.85,'contourLevel':1.,
    11931203            'bondRadius':0.1,'ballScale':0.33,'vdwScale':0.67,'ellipseProb':50,'sizeH':0.50,
    11941204            'unitCellBox':False,'showABC':True,'selectedAtoms':[],
    … …  
    12021212            drawingData = copy.copy(defaultDrawing)
    12031213            drawingData['Atoms'] = []
     1214        if 'contourLevel' not in drawingData:
     1215            drawingData['contourLevel'] = 1.
    12041216        cx,ct,cs = [0,0,0]
    12051217        if generalData['Type'] == 'nuclear':
    … …  
    20072019                bondRadiusTxt.SetLabel('Bond radius, A: '+'%.2f'%(drawingData['bondRadius']))
    20082020                G2plt.PlotStructure(G2frame,data)
     2021               
     2022            def OnContourLevel(event):
     2023                drawingData['contourLevel'] = contourLevel.GetValue()/100.
     2024                contourLevelTxt.SetLabel('Contour level: '+'%.2f'%(drawingData['contourLevel']*generalData['Map']['rhoMax']))
     2025                G2plt.PlotStructure(G2frame,data)
    20092026           
    20102027            slopSizer = wx.BoxSizer(wx.HORIZONTAL)
    2011             slideSizer = wx.FlexGridSizer(6,2)
     2028            slideSizer = wx.FlexGridSizer(7,2)
    20122029            slideSizer.AddGrowableCol(1,1)
    20132030   
    … …  
    20522069            bondRadius.Bind(wx.EVT_SLIDER, OnBondRadius)
    20532070            slideSizer.Add(bondRadius,1,wx.EXPAND|wx.RIGHT)
     2071           
     2072            if generalData['Map']['rhoMax']:
     2073                contourLevelTxt = wx.StaticText(dataDisplay,-1,' Contour level: '+'%.2f'%(drawingData['contourLevel']*generalData['Map']['rhoMax']))
     2074                slideSizer.Add(contourLevelTxt,0,wx.ALIGN_CENTER_VERTICAL)
     2075                contourLevel = wx.Slider(dataDisplay,style=wx.SL_HORIZONTAL,value=int(100*drawingData['contourLevel']))
     2076                contourLevel.SetRange(1,100)
     2077                contourLevel.Bind(wx.EVT_SLIDER, OnContourLevel)
     2078                slideSizer.Add(contourLevel,1,wx.EXPAND|wx.RIGHT)
    20542079           
    20552080            slopSizer.Add(slideSizer,1,wx.EXPAND|wx.RIGHT)
    … …  
    33143339   
    33153340    def OnFourierMaps(event):
     3341        import scipy.fftpack as fft
    33163342        generalData = data['General']
     3343#        for item in generalData:
     3344#            print item,generalData[item]
    33173345        if not generalData['Map']['MapType']:
    33183346            print '**** ERROR - Fourier map not defined'
    33193347            return
    3320         print 'Calculate Fourier maps'
    3321         print generalData['Map']
     3348        mapData = generalData['Map']
     3349        dmin = mapData['Resolution']
     3350        phaseName = generalData['Name']
     3351        SGData = generalData['SGData']
     3352        cell = generalData['Cell'][1:8]       
     3353        A = G2lat.cell2A(cell[:6])
     3354        Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
     3355        Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
     3356        reflName = mapData['RefList']
     3357        if 'PWDR' in reflName:
     3358            PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, reflName)
     3359            reflSets = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
     3360            reflData = reflSets[phaseName]
     3361        elif 'HKLF' in reflName:
     3362            print 'single crystal reflections'
     3363            return
     3364#        Fhkl[0,0,0] = generalData['F000X']
     3365        time0 = time.time()
     3366        for ref in reflData:
     3367            if ref[4] >= dmin:
     3368                for i,hkl in enumerate(ref[11]):
     3369                    hkl = np.asarray(hkl,dtype='i')
     3370                    Fosq,Fcsq,ph = ref[8:11]
     3371                    dp = 360.*ref[12][i]
     3372                    a = cosd(ph+dp)
     3373                    b = sind(ph+dp)
     3374                    phasep = complex(a,b)
     3375                    phasem = complex(a,-b)
     3376                    if 'Fobs' in mapData['MapType']:
     3377                        F = np.sqrt(Fosq)
     3378                        h,k,l = hkl+Hmax
     3379                        Fhkl[h,k,l] = F*phasep
     3380                        h,k,l = -hkl+Hmax
     3381                        Fhkl[h,k,l] = F*phasem
     3382                    elif 'Fcalc' in mapData['MapType']:
     3383                        F = np.sqrt(Fcsq)
     3384                        h,k,l = hkl+Hmax
     3385                        Fhkl[h,k,l] = F*phasep
     3386                        h,k,l = -hkl+Hmax
     3387                        Fhkl[h,k,l] = F*phasem
     3388                    elif 'delt-F' in mapData['MapType']:
     3389                        dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
     3390                        h,k,l = hkl+Hmax
     3391                        Fhkl[h,k,l] = dF*phasep
     3392                        h,k,l = -hkl+Hmax
     3393                        Fhkl[h,k,l] = dF*phasem
     3394                    elif '2*Fo-Fc' in mapData['MapType']:
     3395                        F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
     3396                        h,k,l = hkl+Hmax
     3397                        Fhkl[h,k,l] = F*phasep
     3398                        h,k,l = -hkl+Hmax
     3399                        Fhkl[h,k,l] = F*phasem
     3400                    elif 'Patterson' in mapData['MapType']:
     3401                        h,k,l = hkl+Hmax
     3402                        Fhkl[h,k,l] = complex(Fosq,0.)
     3403                        h,k,l = -hkl+Hmax
     3404                        Fhkl[h,k,l] = complex(Fosq,0.)
     3405        Fhkl = fft.fftshift(Fhkl)
     3406        rho = fft.fftn(Fhkl)/cell[6]
     3407        print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
     3408        mapData['rho'] = np.real(rho)
     3409        mapData['rhoMax'] = np.max(np.real(rho))
     3410        data['Drawing']['contourLevel'] = 1.
     3411        print mapData['MapType']+' computed: rhomax = %.3f rhomin = %.3f'%(mapData['rhoMax'],np.min(np.real(rho)))
     3412## map printing for testing purposes
     3413#        ix,jy,kz = mapData['rho'].shape
     3414#        for k in range(kz):
     3415#            print 'k = ',k
     3416#            for j in range(jy):
     3417#                line = ''
     3418#                if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
     3419#                    line += (jy-j)*'  '
     3420#                for i in range(ix):
     3421#                    r = int(100*mapData['rho'][i,j,k]/mapData['rhoMax'])
     3422#                    line += '%4d'%(r)
     3423#                print line+'\n'
     3424### keep this               
     3425               
    33223426       
    33233427    def OnTextureRefine(event):

  • trunk/GSASIIplot.py

    r515 r520  
    1212import os.path
    1313import numpy as np
     14import numpy.ma as ma
    1415import numpy.linalg as nl
    1516import wx
    … …  
    435436                if G2frame.PickId:
    436437                    found = []
    437                     if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Index Peak List','Unit Cells List','Reflection Lists']:
     438                    if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Index Peak List','Unit Cells List','Reflection Lists'] or \
     439                        'PWDR' in G2frame.PatternTree.GetItemText(PickId):
    438440                        if len(HKL):
    439441                            view = Page.toolbar._views.forward()[0][:2]
    … …  
    499501            else:                                                   #picked a limit line
    500502                G2frame.itemPicked = pick
    501         elif G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists':
     503        elif G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists' or \
     504            'PWDR' in G2frame.PatternTree.GetItemText(PickId):
    502505            G2frame.itemPicked = pick
    503506            pick = str(pick)
    … …  
    540543                G2frame.PatternTree.SetItemPyData(PeakId,data)
    541544                G2pdG.UpdatePeakGrid(G2frame,data)
    542         elif G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists' and xpos:
     545        elif (G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists' or \
     546            'PWDR' in G2frame.PatternTree.GetItemText(PickId)) and xpos:
    543547            Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
    544548            pick = str(G2frame.itemPicked).split('(')[1].strip(')')
    … …  
    738742                else:
    739743                    Plot.axvline(hkl[5],color='r',dashes=(5,5))
    740         elif G2frame.PatternTree.GetItemText(PickId) in ['Reflection Lists']:
     744        elif G2frame.PatternTree.GetItemText(PickId) in ['Reflection Lists'] or \
     745            'PWDR' in G2frame.PatternTree.GetItemText(PickId):
    741746            Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
    742747            for pId,phase in enumerate(Phases):
    … …  
    22712276    drawingData = data['Drawing']
    22722277    drawAtoms = drawingData['Atoms']
     2278    mapData = {}
     2279    rhoXYZ = []
     2280    if 'Map' in generalData:
     2281        mapData = generalData['Map']
     2282        contLevel = drawingData['contourLevel']*mapData['rhoMax']
     2283        if 'delt-F' in mapData['MapType']:
     2284            rho = ma.array(mapData['rho'],mask=(np.abs(mapData['rho'])<contLevel))
     2285        else:
     2286            rho = ma.array(mapData['rho'],mask=(mapData['rho']<contLevel))
     2287        indx = np.array(ma.nonzero(rho)).T
     2288        steps = 1./np.array(rho.shape)
     2289        rhoXYZ = indx*steps
    22732290    cx,ct,cs = drawingData['atomPtrs']
    2274     Wt = [255,255,255]
    2275     Rd = [255,0,0]
    2276     Gr = [0,255,0]
    2277     Bl = [0,0,255]
     2291    Wt = np.array([255,255,255])
     2292    Rd = np.array([255,0,0])
     2293    Gr = np.array([0,255,0])
     2294    Bl = np.array([0,0,255])
    22782295    uBox = np.array([[0,0,0],[1,0,0],[1,1,0],[0,1,0],[0,0,1],[1,0,1],[1,1,1],[0,1,1]])
    22792296    uEdges = np.array([
    … …  
    26072624        glPopMatrix()
    26082625       
     2626    def RenderSmallSphere(x,y,z,radius,color):
     2627        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
     2628        glPushMatrix()
     2629        glTranslate(x,y,z)
     2630        glMultMatrixf(B4mat.T)
     2631        q = gluNewQuadric()
     2632        gluSphere(q,radius,4,2)
     2633        glPopMatrix()
     2634       
    26092635    def RenderEllipsoid(x,y,z,ellipseProb,E,R4,color):
    26102636        s1,s2,s3 = E
    … …  
    26912717        glEnable(GL_LIGHTING)
    26922718        glPopMatrix()
     2719       
     2720    def RenderMap(rhoXYZ,indx,rho,cLevel):
     2721        for i,xyz in enumerate(rhoXYZ):
     2722            x,y,z = xyz
     2723            I,J,K = indx[i]
     2724            alpha = 1.0
     2725            if cLevel < 1.:
     2726                alpha = (abs(rho[I,J,K])/mapData['rhoMax']-cLevel)/(1.-cLevel)
     2727            if rho[I,J,K] < 0.:
     2728                RenderSmallSphere(x,y,z,0.1*alpha,Rd)
     2729            else:
     2730                RenderSmallSphere(x,y,z,0.1*alpha,Bl)
    26932731                           
    26942732    def Draw():
    … …  
    27392777        BackboneColor = []
    27402778        time0 = time.time()
     2779#        glEnable(GL_BLEND)
     2780#        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
    27412781        for iat,atom in enumerate(drawingData['Atoms']):
    27422782            x,y,z = atom[cx:cx+3]
    … …  
    27512791            if iat in Ind:
    27522792                color = np.array(Gr)/255.
     2793#            color += [.25,]
    27532794            radius = 0.5
    27542795            if atom[cs] != '':
    2755                 glLoadName(atom[-3])                   
     2796                try:
     2797                    glLoadName(atom[-3])
     2798                except: #problem with old files - missing code
     2799                    pass                   
    27562800            if 'balls' in atom[cs]:
    27572801                vdwScale = drawingData['vdwScale']
    … …  
    28202864            elif atom[cs+1] == 'chain' and atom[ct-1] == 'CA':
    28212865                RenderLabel(x,y,z,atom[ct-2],radius)
     2866#        glDisable(GL_BLEND)
     2867        if len(rhoXYZ):
     2868            RenderMap(rhoXYZ,indx,rho,drawingData['contourLevel'])
    28222869        if Backbone:
    28232870            RenderBackbone(Backbone,BackboneColor,bondR)
Note: See TracChangeset for help on using the changeset viewer.