Changeset 1996


Ignore:
Timestamp:
Oct 9, 2015 2:57:18 PM (7 years ago)
Author:
vondreele
Message:

fix Fade problem
ModulationDeriv? work
twin sf error fixed

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImath.py

    r1995 r1996  
    965965    else:
    966966        Fmod = 1.0           
    967     if Ax.shape[1]:
    968         tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
    969         XmodA = Ax[:,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X pos X 32
    970         XmodB = Bx[:,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
    971         Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
     967    tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
     968    XmodA = Ax[:,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X pos X 32
     969    XmodB = Bx[:,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
     970    Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
    972971    if Au.shape[1]:
    973972        tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
    … …  
    982981    sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1)       #refBlk X ops x atoms; sum for G-L integration
    983982    cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1)       #ditto
    984 #        GSASIIpath.IPyBreak()                 
    985 
     983#    GSASIIpath.IPyBreak()
    986984    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
    987985   
    … …  
    994992    Mast: array orthogonalization matrix for Uij
    995993    '''
    996     import scipy.misc as sm
    997    
     994   
     995    nxs = np.newaxis
     996    numeric = False   
    998997    cosHA,sinHA = Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast)
    999     deltx = 0.00001
    1000     deltf = 0.0001
    1001     deltu = 0.000001
    1002     Mf = FSSdata.T.shape
     998    Mf = [H.shape[0],]+list(FSSdata.T.shape)    #ops x atoms x waves x 2 (sin+cos frac mods)
    1003999    dGdMfC = np.zeros(Mf)
    10041000    dGdMfS = np.zeros(Mf)
    1005     Mx = XSSdata.T.shape
    1006     Mx = [H.shape[0],]+list(Mx)   #atoms x waves x sin+cos pos mods
     1001    Mx = [H.shape[0],]+list(XSSdata.T.shape)   #ops x atoms x waves x 6 (sin+cos pos mods)
    10071002    dGdMxC = np.zeros(Mx)
    10081003    dGdMxS = np.zeros(Mx)
    1009     for i in range(Mx[1]):  #atoms
    1010         for j in range(Mx[2]):  #waves
    1011             for k in range(Mx[3]):  #coeff
    1012                 XSSdata[k,j,i] += deltx
    1013                 Cap,Sap = np.squeeze(Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast))
    1014                 XSSdata[k,j,i] -= 2*deltx
    1015                 Cam,Sam = np.squeeze(Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast))
    1016                 XSSdata[k,j,i] += deltx
    1017                 dGdMxC[:,:,j,k] += (Cap-Cam)/(2*deltx)
    1018                 dGdMxS[:,:,j,k] += (Sap-Sam)/(2*deltx)
    1019     Mu = USSdata.T.shape
     1004    Mu = [H.shape[0],]+list(USSdata.T.shape)    #ops x atoms x waves x 12 (sin+cos Uij mods)
    10201005    dGdMuC = np.zeros(Mu)
    10211006    dGdMuS = np.zeros(Mu)
    1022 #    GSASIIpath.IPyBreak()                 
    1023                
    1024    
    1025 #    nxs = np.newaxis
    1026 #    glTau,glWt = pwd.pygauleg(0.,1.,32)         #get Gauss-Legendre intervals & weights
    1027 #    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
    1028 #    dGdBx = np.zeros_like(Bx)
    1029 #    Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
    1030 #    dGdAf = np.zeros_like(Af)
    1031 #    Bf = np.array(FSSdata[1]).T    #cos frac mods...
    1032 #    dGdBf = np.zeros_like(Bf)
    1033 #    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
    1034 #    dGdAu = np.zeros_like(Au)
    1035 #    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
    1036 #    dGdBu = np.zeros_like(Bu)
    1037 #   
    1038 #    GSASIIpath.IPyBreak()                 
    1039 #    if 'Fourier' in waveTypes:
    1040 #        nf = 0
    1041 #        nx = 0
    1042 #        XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32))
    1043 #        FmodZ = np.zeros((Af.shape[0],Af.shape[1],32))
    1044 #        if 'Crenel' in waveTypes:
    1045 #            nC = np.where('Crenel' in waveTypes)
    1046 #            nf = 1
    1047 #            #FmodZ = 0   replace
    1048 #    else:
    1049 #        nx = 1
    1050 #        if 'Sawtooth' in waveTypes:
    1051 #            nS = np.where('Sawtooth' in waveTypes)
    1052 #            #XmodZ = 0   replace
    1053 #    if Ax.shape[1]:
    1054 #        tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
    1055 #        StauX = np.ones_like(Ax)[:,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #also dXmod/dAx
    1056 #        CtauX = np.ones_like(Bx)[:,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #also dXmod/dBx
    1057 #        XmodA = Ax[:,nx:,:,nxs]*StauX #atoms X waves X pos X 32
    1058 #        XmodB = Bx[:,nx:,:,nxs]*CtauX #ditto
    1059 #        Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
    1060 #        D = H[:,3][:,nxs]*glTau[nxs,:]              #m*tau; ops X 32
    1061 #        HdotX = np.inner(H[:,:3],np.swapaxes(Xmod,1,2))+D[:,nxs,:]     #ops x atoms X 32
    1062 #    if Af.shape[1]:
    1063 #        tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau  #Fwaves x 32
    1064 #        StauF = np.ones_like(Af)[:,nf:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf
    1065 #        CtauF = np.ones_like(Bf)[:,nf:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf
    1066 #        FmodA = Af[:,nf:,nxs]*StauF   #atoms X Fwaves X 32
    1067 #        FmodB = Bf[:,nf:,nxs]*CtauF
    1068 #        Fmod = np.sum(FmodA+FmodB+FmodC,axis=1)             #atoms X 32; sum waves
    1069 #    else:
    1070 #        Fmod = np.ones_like(HdotX)           
    1071 #    if Au.shape[1]:
    1072 #        tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
    1073 #        StauU = np.ones_like(Au)[:,:,:,:,nxs]*np.sin(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmod/dAu
    1074 #        CtauU = np.ones_like(Bu)[:,:,:,:,nxs]*np.cos(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmod/dBu
    1075 #        UmodA = Au[:,:,:,:,nxs]*StauU #atoms x waves x 3x3 x 32
    1076 #        UmodB = Bu[:,:,:,:,nxs]*CtauU #ditto
    1077 #        Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x 32; sum waves
    1078 #        HbH = np.exp(-np.sum(H[:,nxs,nxs,:3]*np.inner(H[:,:3],Umod),axis=-1)) # ops x atoms x 32 add Overhauser corr.?
    1079 #    else:
    1080 #        HbH = np.ones_like(HdotX)
    1081 #    HdotXA = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodA,-1,-2)[nxs,:,:,:,:]  #ops x atoms x waves x 32 x xyz
    1082 #    HdotXB = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodB,-1,-2)[nxs,:,:,:,:]
     1007    if numeric:        #numeric derivatives - slow!! works for pos waves
     1008        deltx = 0.00001
     1009        deltf = 0.0001
     1010        deltu = 0.000001
     1011        for i in range(Mx[1]):  #atoms
     1012            for j in range(Mx[2]):  #waves
     1013                for k in range(Mx[3]):  #coeff
     1014                    XSSdata[k,j,i] += deltx
     1015                    Cap,Sap = np.squeeze(Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast))
     1016                    XSSdata[k,j,i] -= 2*deltx
     1017                    Cam,Sam = np.squeeze(Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast))
     1018                    XSSdata[k,j,i] += deltx
     1019                    dGdMxC[:,:,j,k] += (Cap-Cam)/(2*deltx)
     1020                    dGdMxS[:,:,j,k] += (Sap-Sam)/(2*deltx)
     1021    else:   #analytic derivatives                 
     1022        glTau,glWt = pwd.pygauleg(0.,1.,32)         #get Gauss-Legendre intervals & weights
     1023        Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
     1024        Bf = np.array(FSSdata[1]).T    #cos frac mods...
     1025        Ax = np.array(XSSdata[:3]).T   #...cos pos mods
     1026        Bx = np.array(XSSdata[3:]).T   #...cos pos mods
     1027        Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
     1028        Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
     1029       
     1030        if 'Fourier' in waveTypes:
     1031            nf = 0
     1032            nx = 0
     1033            XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32))
     1034            FmodZ = np.zeros((Af.shape[0],Af.shape[1],32))
     1035            if 'Crenel' in waveTypes:
     1036                nC = np.where('Crenel' in waveTypes)
     1037                nf = 1
     1038                #FmodZ = 0   replace
     1039        else:
     1040            nx = 1
     1041            if 'Sawtooth' in waveTypes:
     1042                nS = np.where('Sawtooth' in waveTypes)
     1043                #XmodZ = 0   replace
     1044        tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
     1045        StauX = np.ones_like(Ax)[:,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #also dXmod/dAx
     1046        CtauX = np.ones_like(Bx)[:,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #also dXmod/dBx
     1047        XmodA = Ax[:,nx:,:,nxs]*StauX #atoms X waves X pos X 32
     1048        XmodB = Bx[:,nx:,:,nxs]*CtauX #ditto
     1049        Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
     1050        D = H[:,3][:,nxs]*glTau[nxs,:]              #m*tau; ops X 32
     1051        HdotX = np.inner(H[:,:3],np.swapaxes(Xmod,1,2))+D[:,nxs,:]     #ops x atoms X 32
     1052        if Af.shape[1]:
     1053            tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau  #Fwaves x 32
     1054            StauF = np.ones_like(Af)[:,nf:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf
     1055            CtauF = np.ones_like(Bf)[:,nf:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf
     1056            FmodA = Af[:,nf:,nxs]*StauF   #atoms X Fwaves X 32
     1057            FmodB = Bf[:,nf:,nxs]*CtauF
     1058            Fmod = np.sum(FmodA+FmodB+FmodC,axis=1)             #atoms X 32; sum waves
     1059        else:
     1060            Fmod = np.ones_like(HdotX)           
     1061        if Au.shape[1]:
     1062            tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
     1063            StauU = np.ones_like(Au)[:,:,:,:,nxs]*np.sin(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmod/dAu
     1064            CtauU = np.ones_like(Bu)[:,:,:,:,nxs]*np.cos(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmod/dBu
     1065            UmodA = Au[:,:,:,:,nxs]*StauU #atoms x waves x 3x3 x 32
     1066            UmodB = Bu[:,:,:,:,nxs]*CtauU #ditto
     1067            Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x 32; sum waves
     1068            HbH = np.exp(-np.sum(H[:,nxs,nxs,:3]*np.inner(H[:,:3],Umod),axis=-1)) # ops x atoms x 32 add Overhauser corr.?
     1069        else:
     1070            HbH = np.ones_like(HdotX)
     1071        HdotXA = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodA,-1,-2)[nxs,:,:,:,:]  #ops x atoms x waves x 32 x xyz
     1072        HdotXB = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodB,-1,-2)[nxs,:,:,:,:]
     1073        sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1)       #ops x atoms x waves x xyz; sum for G-L integration
     1074        cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1)       #ditto
     1075#    GSASIIpath.IPyBreak()
    10831076#    dHdXA = H[:,nxs,nxs,nxs,:3]*np.swapaxes(StauX,-1,-2)[nxs,:,:,:,:]    #ops x atoms x waves x 32 x xyz
    10841077#    dHdXB = H[:,nxs,nxs,nxs,:3]*np.swapaxes(CtauX,-1,-2)[nxs,:,:,:,:]              #ditto
    1085 #    sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1)       #ops x atoms x waves x xyz; sum for G-L integration
    1086 #    cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1)       #ditto
    10871078#    dGdAx = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(twopi*dHdXA*np.sin(twopi*HdotXA)+np.cos(twopi*HdotXA))*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    10881079## ops x atoms x waves x xyz - imaginary part

  • trunk/GSASIIplot.py

    r1993 r1996  
    50505050        glShadeModel(GL_SMOOTH)
    50515051                           
    5052     def Draw(caller='',Fade=None):
     5052    def Draw(caller='',Fade=[]):
    50535053#useful debug?       
    50545054#        if caller:
    … …  
    50625062            pageName = G2frame.dataDisplay.GetPageText(page)
    50635063        rhoXYZ = []
     5064        rho = []
    50645065        if len(D4mapData.get('rho',[])):        #preferentially select 4D map if there
    50655066            rho = D4mapData['rho'][:,:,:,int(G2frame.tau*10)]   #pick current tau 3D slice
    … …  
    51295130#        glEnable(GL_BLEND)
    51305131#        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
    5131         if Fade == None:
     5132        if not len(Fade):
    51325133            atmFade = np.ones(len(drawingData['Atoms']))
    51335134        else:

  • trunk/GSASIIstrIO.py

    r1978 r1996  
    22892289                    else:
    22902290                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
    2291                         controlDict[pfx+'TwinNMN'] = twin[1][2]
     2291                        controlDict[pfx+'TwinNMN'] = twin[1][1]
    22922292                    if Twins[0][1][1]:
    22932293                        hapVary.append(pfx+'TwinFr:'+str(it))

  • trunk/GSASIIstrMath.py

    r1993 r1996  
    754754        else:
    755755            if len(TwinLaw) > 1:
    756                 refl.T[9] = np.sum(fas[:,:,0]**2,axis=0)+np.sum(fbs[:,:,0]**2,axis=0)   #FcT from primary twin element
     756                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
    757757                refl.T[7] = np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fas,axis=0)**2,axis=-1)+   \
    758758                    np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
    … …  
    10671067        fag = fa*GfpuA[0]-fb*GfpuA[1]
    10681068        fbg = fb*GfpuA[0]+fa*GfpuA[1]
    1069 #        GSASIIpath.IPyBreak()
    10701069        fas = np.sum(np.sum(fag,axis=-1),axis=-1)
    10711070        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
     1071#        GSASIIpath.IPyBreak()
    10721072        if 'P' in calcControls[hfx+'histType']:
    10731073            refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
    … …  
    10751075        else:
    10761076            if len(TwinLaw) > 1:
    1077                 refl.T[10] = np.sum(fas[:,:,0]**2,axis=0)+np.sum(fbs[:,:,0]**2,axis=0)   #FcT from primary twin element
     1077                refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
    10781078                refl.T[8] = np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fas,axis=0)**2,axis=-1)+   \
    10791079                    np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
Note: See TracChangeset for help on using the changeset viewer.