Changeset 1966


Ignore:
Timestamp:
Aug 28, 2015 3:02:12 PM (7 years ago)
Author:
vondreele
Message:

new 3D HKL plot commands for selecting view axes & set a zone plot - trouble with rotations tho.
SS structure factor by integration, not Bessel fxns.

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIgrid.py

    r1958 r1966  
    28312831        Hmax = np.array([int(np.max(refList.T[0])),int(np.max(refList.T[1])),int(np.max(refList.T[2]))])
    28322832        Vpoint = [int(np.mean(refList.T[0])),int(np.mean(refList.T[1])),int(np.mean(refList.T[2]))]
    2833         controls = {'Type' : 'Fosq','Iscale' : False,'HKLmax' : Hmax,'HKLmin' : Hmin,
     2833        controls = {'Type' : 'Fosq','Iscale' : False,'HKLmax' : Hmax,'HKLmin' : Hmin,'Zone':False,
    28342834            'FoMax' : FoMax,'Scale' : 1.0,'Drawing':{'viewPoint':[Vpoint,[]],'default':Vpoint[:],
    2835             'backColor':[0,0,0],'depthFog':False,'Zclip':10.0,'cameraPos':10.,'Zstep':0.05,
    2836             'Scale':1.0,'oldxy':[],'viewDir':[1,0,0]},'Super':Super,'SuperVec':SuperVec}
     2835            'backColor':[0,0,0],'depthFog':False,'Zclip':10.0,'cameraPos':10.,'Zstep':0.05,'viewUp':[0,1,0],
     2836            'Scale':1.0,'oldxy':[],'viewDir':[0,0,1]},'Super':Super,'SuperVec':SuperVec}
    28372837        G2plt.Plot3DSngl(G2frame,newPlot=True,Data=controls,hklRef=refList,Title=phaseName)
    28382838       
    … …  
    28612861        Hmax = np.array([int(np.max(refList.T[0])),int(np.max(refList.T[1])),int(np.max(refList.T[2]))])
    28622862        Vpoint = [int(np.mean(refList.T[0])),int(np.mean(refList.T[1])),int(np.mean(refList.T[2]))]
    2863         controls = {'Type' : 'Fosq','Iscale' : False,'HKLmax' : Hmax,'HKLmin' : Hmin,
     2863        controls = {'Type' : 'Fosq','Iscale' : False,'HKLmax' : Hmax,'HKLmin' : Hmin,'Zone':False,
    28642864            'FoMax' : FoMax,'Scale' : 1.0,'Drawing':{'viewPoint':[Vpoint,[]],'default':Vpoint[:],
    2865             'backColor':[0,0,0],'depthFog':False,'Zclip':10.0,'cameraPos':10.,'Zstep':0.05,
     2865            'backColor':[0,0,0],'depthFog':False,'Zclip':10.0,'cameraPos':10.,'Zstep':0.05,'viewUp':[0,1,0],
    28662866            'Scale':1.0,'oldxy':[],'viewDir':[1,0,0]},'Super':Super,'SuperVec':SuperVec}
    28672867        G2plt.Plot3DSngl(G2frame,newPlot=True,Data=controls,hklRef=refList,Title=phaseName)

  • trunk/GSASIImath.py

    r1958 r1966  
    924924################################################################################
    925925   
    926 def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM):
     926def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast):
    927927    import scipy.special as sp
    928928    import scipy.integrate as si
     929   
     930    def cosMod(tau,H,A,B,S):
     931        uTau = posFourier(tau,A,B,S)
     932        SdotU = twopi*(np.inner(uTau.T,H[:3].T))+H[3]*tau
     933        return list(np.cos(SdotU))[0]
     934       
     935    def sinMod(tau,H,A,B,S):
     936        uTau = posFourier(tau,A,B,S)
     937        SdotU = twopi*(np.inner(uTau.T,H[:3]))+H[3]*tau
     938        return list(np.sin(SdotU))[0]
     939               
     940    def expModInt(H,A,B,S):
     941        intReal = np.array([[si.romberg(cosMod,0,1,args=(h,a,b,s),rtol=1.e-1) for a,b in zip(A,B)] for h,s in zip(H,S)])
     942        intImag = np.array([[si.romberg(sinMod,0,1,args=(h,a,b,s),rtol=1.e-1) for a,b in zip(A,B)] for h,s in zip(H,S)])
     943#        return np.sqrt(intReal**2+intImag**2)
     944        return intReal,intImag
     945
    929946    Smult,TauT = SStauM             # both atoms x SGops
    930947    m = SSUniq.T[3]
    931     A = np.array(XSSdata[:3])   #sin mods x waves x atoms
    932     B = np.array(XSSdata[3:])   #cos mods...
    933     if XSSdata.ndim > 2:
    934         nh = np.arange(XSSdata.shape[1])        #0..max harmonic order-1
    935         M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis])     # waves x SGops
    936         HdotA = (np.inner(A.T,SSUniq.T[:3].T))    #atoms X waves x SGops
    937         HdotB = (np.inner(B.T,SSUniq.T[:3].T))    # ditto
    938         Hdot = np.sqrt(HdotA**2+HdotB**2)
    939         Hphi = np.sum(np.arctan2(HdotB,HdotA)*M,axis=1)     
    940         GpA = np.sum(sp.jn(-M,twopi*Hdot),axis=1)                      # sum on waves
     948    Ax = np.array(XSSdata[:3]).T   #atoms x waves x sin pos mods
     949    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
     950    Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
     951    Bf = np.array(FSSdata[1]).T    #cos frac mods...
     952    Ab = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
     953    Bb = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
     954#    GSASIIpath.IPyBreak()       
     955    if Ax.ndim > 2:
     956        GpA = np.array(expModInt(SSUniq,Ax,Bx,Smult))
    941957    else:
    942         HdotA = np.inner(A.T,SSUniq.T[:3].T)    #atoms X SGops
    943         HdotB = np.inner(B.T,SSUniq.T[:3].T)    # ditto
    944         Hdot = np.sqrt(HdotA**2+HdotB**2)
    945         Hphi = np.arctan2(HdotB,HdotA)       
    946         GpA = sp.jn(-m,twopi*Hdot)                      # atoms x SGops
    947        
    948 #    GSASIIpath.IPyBreak()
    949     return GpA.T,Hphi.T             # SGops x atoms
    950    
    951 def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM):
     958        GpA = np.array(expModInt(SSUniq,Ax[:,np.newaxis,:],Bx[:,np.newaxis,:],Smult))       
     959    return GpA             # SGops x atoms
     960   
     961def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast):
    952962    import scipy.special as sp
    953963    Smult,TauT = SStauM             # both atoms x SGops
    … …  
    10631073    return drawAtoms,Fade
    10641074   
     1075def BessJn(nmax,x):
     1076    ''' compute Bessel function J(n,x) from scipy routine & recurrance relation
     1077    returns sequence of J(n,x) for n in range [-nmax...0...nmax]
     1078   
     1079    :param  integer nmax: maximul order for Jn(x)
     1080    :param  float x: argument for Jn(x)
     1081   
     1082    :returns numpy array: [J(-nmax,x)...J(0,x)...J(nmax,x)]
     1083   
     1084    '''
     1085    import scipy.special as sp
     1086    bessJn = np.zeros(2*nmax+1)
     1087    bessJn[nmax] = sp.j0(x)
     1088    bessJn[nmax+1] = sp.j1(x)
     1089    bessJn[nmax-1] = -bessJn[nmax+1]
     1090    for i in range(2,nmax+1):
     1091        bessJn[i+nmax] = 2*(i-1)*bessJn[nmax+i-1]/x-bessJn[nmax+i-2]
     1092        bessJn[nmax-i] = bessJn[i+nmax]*(-1)**i
     1093    return bessJn
     1094   
     1095def BessIn(nmax,x):
     1096    ''' compute modified Bessel function I(n,x) from scipy routines & recurrance relation
     1097    returns sequence of I(n,x) for n in range [-nmax...0...nmax]
     1098   
     1099    :param  integer nmax: maximul order for In(x)
     1100    :param  float x: argument for In(x)
     1101   
     1102    :returns numpy array: [I(-nmax,x)...I(0,x)...I(nmax,x)]
     1103   
     1104    '''
     1105    import scipy.special as sp
     1106    bessIn = np.zeros(2*nmax+1)
     1107    bessIn[nmax] = sp.i0(x)
     1108    bessIn[nmax+1] = sp.i1(x)
     1109    bessIn[nmax-1] = bessIn[nmax+1]
     1110    for i in range(2,nmax+1):
     1111        bessIn[i+nmax] = bessIn[nmax+i-2]-2*(i-1)*bessIn[nmax+i-1]/x
     1112        bessIn[nmax-i] = bessIn[i+nmax]
     1113    return bessIn
     1114       
    10651115   
    10661116################################################################################

  • trunk/GSASIIplot.py

    r1958 r1966  
    572572        global ifBox
    573573        Choice = {'F':'Fo','S':'Fosq','U':'Unit','D':'dFsq','W':'dFsq/sig'}
     574        viewChoice = {'H':[[1,0,0],[0,0,1]],'K':[[0,1,0],[1,0,0]],'L':[[0,0,1],[0,1,0]]}
    574575        try:
    575576            keyCode = event.GetKeyCode()
    … …  
    579580        except AttributeError:       #if from OnKeyBox above
    580581            key = str(event.key).upper()
    581         if key in ['C']:
     582        if key in ['C','H','K','L']:
     583            if key == 'C':
     584                Data['Zone'] = False
     585                key = 'L'
    582586            drawingData['viewPoint'][0] = drawingData['default']
    583             drawingData['viewDir'] = [0,0,1]
     587            drawingData['viewDir'] = viewChoice[key][0]
     588            drawingData['viewUp'] = viewChoice[key][1]
    584589            drawingData['oldxy'] = []
    585             V0 = np.array([0,0,1])
    586             V = np.inner(Amat,V0)
    587             V /= np.sqrt(np.sum(V**2))
    588             A = np.arccos(np.sum(V*V0))
    589             Q = G2mth.AV2Q(A,[0,1,0])
    590             drawingData['Quaternion'] = Q
    591             Q = drawingData['Quaternion']
     590#            V0 = np.array(viewChoice[key][0])
     591#            V = np.inner(Amat,V0)
     592#            V /= np.sqrt(np.sum(V**2))
     593#            A = np.arccos(np.sum(V*V0))
     594#            Q = G2mth.AV2Q(A,viewChoice[key][1])
     595            drawingData['Quaternion'] = [-1,0,0,0]  #the old stuff always gave this...
     596        elif key in 'Z':
     597            Data['Zone'] = not Data['Zone']
    592598        elif key in 'B':
    593599            ifBox = not ifBox
    … …  
    883889        cPos = drawingData['cameraPos']
    884890        Zclip = drawingData['Zclip']*cPos/20.
     891        if Data['Zone']:
     892            Zclip = 0.1
    885893        Q = drawingData['Quaternion']
    886894        Tx,Ty,Tz = drawingData['viewPoint'][0]
    … …  
    901909        glViewport(0,0,VS[0],VS[1])
    902910        gluPerspective(20.,aspect,cPos-Zclip,cPos+Zclip)
    903         gluLookAt(0,0,cPos,0,0,0,0,1,0)
     911        vDir = drawingData['viewDir']
     912        vUp = drawingData['viewUp']
     913#        gluLookAt(0,0,cPos,0,0,0,0,1,0)
     914        gluLookAt(vDir[0]*cPos,vDir[1]*cPos,vDir[2]*cPos, 0,0,0, vUp[0],vUp[1],vUp[2])
    904915        SetLights()           
    905916           
    … …  
    936947    Page.SetFocus()
    937948    Page.Choice = None
    938     choice = [' save as/key:','jpeg','tiff','bmp','c: recenter to default','b: toggle box ','+: increase scale',
    939     '-: decrease scale','f: Fobs','s: Fobs**2','u: unit','d: Fo-Fc','w: DF/sig','i: toggle intensity scaling']
     949    choice = [' save as/key:','jpeg','tiff','bmp','h: view down h','k: view down k','l: view down l',
     950    'z: zero zone toggle','c: reset to default','b: toggle box ','+: increase scale','-: decrease scale',
     951    'f: Fobs','s: Fobs**2','u: unit','d: Fo-Fc','w: DF/sig','i: toggle intensity scaling']
    940952    cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,choices=choice)
    941953    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)

  • trunk/GSASIIpwdGUI.py

    r1959 r1966  
    31213121        Hmax = np.array([int(np.max(refList.T[0])),int(np.max(refList.T[1])),int(np.max(refList.T[2]))])
    31223122        Vpoint = [int(np.mean(refList.T[0])),int(np.mean(refList.T[1])),int(np.mean(refList.T[2]))]
    3123         controls = {'Type':'Fosq','Iscale':False,'HKLmax':Hmax,'HKLmin':Hmin,
     3123        controls = {'Type':'Fosq','Iscale':False,'HKLmax':Hmax,'HKLmin':Hmin,'Zone':False,
    31243124            'FoMax' : FoMax,'Scale' : 1.0,'Drawing':{'viewPoint':[Vpoint,[]],'default':Vpoint[:],
    3125             'backColor':[0,0,0],'depthFog':False,'Zclip':10.0,'cameraPos':10.,'Zstep':0.05,
    3126             'Scale':1.0,'oldxy':[],'viewDir':[1,0,0]},'Super':Super,'SuperVec':SuperVec}
     3125            'backColor':[0,0,0],'depthFog':False,'Zclip':10.0,'cameraPos':10.,'Zstep':0.05,'viewUp':[0,1,0],
     3126            'Scale':1.0,'oldxy':[],'viewDir':[0,0,1]},'Super':Super,'SuperVec':SuperVec}
    31273127        G2plt.Plot3DSngl(G2frame,newPlot=True,Data=controls,hklRef=refList,Title=phaseName)
    31283128       

  • trunk/GSASIIstrMath.py

    r1958 r1966  
    814814    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
    815815        Flack = 1.-2.*parmDict[phfx+'Flack']
     816    time0 = time.time()
     817    nref = len(refDict['RefList'])/100   
    816818    for iref,refl in enumerate(refDict['RefList']):
    817819        if 'T' in calcControls[hfx+'histType']:
    … …  
    915917                dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
    916918                    2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     919    print ' derivative time %.4f\r'%(time.time()-time0)
    917920        #loop over atoms - each dict entry is list of derivatives for all the reflections
    918921    if nTwin > 1:
    … …  
    9981001            dat = G2el.getFFvalues(FFtables,0.)       
    9991002        refDict['FF']['El'] = dat.keys()
    1000         refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))   
     1003        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
     1004    time0 = time.time()
     1005    nref = len(refDict['RefList'])/100   
    10011006    for iref,refl in enumerate(refDict['RefList']):
    10021007        if 'NT' in calcControls[hfx+'histType']:
    … …  
    10281033            Phi = np.hstack((Phi,-Phi))
    10291034            SSPhi = np.hstack((SSPhi,-SSPhi))
    1030         GfpuA,Hphi = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM)       
     1035        GfpuA = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast)
    10311036        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
    10321037        sinp = np.sin(phase)
    … …  
    10401045        fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
    10411046        fa *= GfpuA
    1042         fb *= GfpuA        
     1047        fb *= GfpuA       
    10431048        fas = np.real(np.sum(np.sum(fa,axis=1),axis=1))
    10441049        fbs = np.real(np.sum(np.sum(fb,axis=1),axis=1))
    … …  
    10481053        refl[7+im] = np.sum(fasq)+np.sum(fbsq)
    10491054        refl[10+im] = atan2d(fbs[0],fas[0])
     1055        if not iref%nref: print 'ref no. %d time %.4f\r'%(iref,time.time()-time0),
     1056    print '\nref no. %d time %.4f\r'%(iref,time.time()-time0)
    10501057#        GSASIIpath.IPyBreak()
    10511058   
    … …  
    10941101        Phi = np.inner(H[:3],SGT)
    10951102        SSPhi = np.inner(H,SSGT)
    1096         GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM)
    1097         dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM)
     1103        GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast)
     1104        dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast)
    10981105        #need ModulationDerv here dGAdXsin, etc 
    10991106        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
Note: See TracChangeset for help on using the changeset viewer.