Changeset 1928


Ignore:
Timestamp:
Jul 15, 2015 2:34:19 PM (7 years ago)
Author:
vondreele
Message:

fix error in Fourier calc.
new getRho routine
new add OH hydrogens looks in delt-F map
put xyz & rho of view point on status line

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1927 r1928  
    490490            Hpos = np.array([[a,b*cosd(x),b*sind(x)] for x in azm])
    491491            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
    492             return Hpos
     492            Rhos = np.array([getRho(pos,mapData) for pos in Hpos])
     493            imax = np.argmax(Rhos)
     494            return [Hpos[imax],]
    493495    return []
    494496       
     
    20492051                hkl = np.asarray(hkl,dtype='i')
    20502052                dp = 360.*Phi[i]                #and phi
    2051                 a = cosd(ph)
    2052                 b = sind(ph)
    2053                 phasep = complex(a,b)+dp
    2054                 phasem = complex(a,-b)-dp
     2053                a = cosd(ph+dp)
     2054                b = sind(ph+dp)
     2055                phasep = complex(a,b)
     2056                phasem = complex(a,-b)
    20552057                if 'Fobs' in mapData['MapType']:
    20562058                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
     
    24562458    return mapData,map4DData
    24572459   
     2460def getRho(xyz,mapData):
     2461    ''' get scattering density at a point by 8-point interpolation
     2462    param xyz: coordinate to be probed
     2463    param: mapData: dict of map data
     2464   
     2465    :returns: density at xyz
     2466    '''
     2467    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
     2468    if not len(mapData):
     2469        return 0.0
     2470    rho = copy.copy(mapData['rho'])     #don't mess up original
     2471    if not len(rho):
     2472        return 0.0
     2473    mapShape = np.array(rho.shape)
     2474    mapStep = 1./mapShape
     2475    X = np.array(xyz)%1.    #get into unit cell
     2476    I = np.array(X*mapShape,dtype='int')
     2477    return rho[I[0],I[1],I[2]]
     2478    D = X-I*mapStep         #position inside map cell
     2479    D12 = D[0]*D[1]
     2480    D13 = D[0]*D[2]
     2481    D23 = D[1]*D[2]
     2482    D123 = np.prod(D)
     2483    Rho = rollMap(rho,I)    #shifts map so point is in corner
     2484    R = Rho[0,0,0]*(1.-np.sum(D))+Rho[1,0,0]*D[0]+Rho[0,1,0]*D[1]+Rho[0,0,1]*D[2]+  \
     2485        Rho[1,1,1]*D123+Rho[0,1,1]*(D23-D123)+Rho[1,0,1]*(D13-D123)+Rho[1,1,0]*(D12-D123)+  \
     2486        Rho[0,0,0]*(D12+D13+D23-D123)-Rho[0,0,1]*(D13+D23-D123)-    \
     2487        Rho[0,1,0]*(D23+D12-D123)-Rho[1,0,0]*(D13+D12-D123)   
     2488    return R
     2489       
    24582490def SearchMap(generalData,drawingData,Neg=False):
    24592491    '''Does a search of a density map for peaks meeting the criterion of peak
  • trunk/GSASIIphsGUI.py

    r1927 r1928  
    15911591                    neigh[1][1].append(nextneigh[1][1][0])
    15921592                neigh[2] = max(0,nH)  #set expected no. H's needed
    1593                 AddHydIds.append(neigh[1][1])
    1594                 Neigh.append(neigh)
     1593                if len(neigh[1][0]):
     1594                    AddHydIds.append(neigh[1][1])
     1595                    Neigh.append(neigh)
    15951596            if Neigh:
     1597                mapError = False
    15961598                dlg = G2gd.AddHatomDialog(G2frame,Neigh,data)
    15971599                if dlg.ShowModal() == wx.ID_OK:
     
    15991601                    Neigh = dlg.GetData()
    16001602                    mapData = generalData['Map']
    1601                     mapError = False
    1602                     for ineigh,neigh in enumerate(Neigh):                       
     1603                    for ineigh,neigh in enumerate(Neigh):
    16031604                        AddHydIds[ineigh].append(neigh[2])
    1604                         if 'O' in neigh[0] and not len(mapData['rho']) or not 'delt-F' in mapData['MapType']:
     1605                        if 'O' in neigh[0] and (not len(mapData['rho']) or not 'delt-F' in mapData['MapType']):
    16051606                            mapError = True
    16061607                            continue                           
  • trunk/GSASIIplot.py

    r1909 r1928  
    45214521                    SetTranslation(newxy)
    45224522                    Tx,Ty,Tz = drawingData['viewPoint'][0]
    4523                     G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f'%(Tx,Ty,Tz),1)
     4523                    rho = G2mth.getRho([Tx,Ty,Tz],mapData)
     4524                    G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f; density: %.4f'%(Tx,Ty,Tz,rho),1)
    45244525                elif event.MiddleIsDown():
    45254526                    SetRotationZ(newxy)
     
    52295230        choice += ['+: increase tau','-: decrease tau','0: set tau = 0']
    52305231
     5232    Tx,Ty,Tz = drawingData['viewPoint'][0]
     5233    rho = G2mth.getRho([Tx,Ty,Tz],mapData)
     5234    G2frame.G2plotNB.status.SetStatusText('View point: %.4f, %.4f, %.4f; density: %.4f'%(Tx,Ty,Tz,rho),1)
     5235
    52315236    cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,choices=choice)
    52325237    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
Note: See TracChangeset for help on using the changeset viewer.