Changeset 2070


Ignore:
Timestamp:
Nov 25, 2015 2:19:24 PM (6 years ago)
Author:
vondreele
Message:

remove extraneous argument from G2mth.posFourier
work on ZigZag/Block? SS functions; function OK, derivatives not.

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r2067 r2070  
    957957        nx = 0
    958958        if 'ZigZag' in waveTypes[iatm]:
    959             nt = 1
     959            nx = 1
    960960            Tmm = Ax[iatm][0][:2]                       
    961961            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
     
    963963            slopeUp = 2.*XYZmax/DT
    964964            slopeDn = 2.*XYZmax/(1.-DT)
    965             XmodZ[iatm][0] = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax+slopeUp*((t-Tmm[0])%1.),XYZmax-slopeDn*((t-Tmm[1])%1.)) for t in glTau]).T
     965            XmodZ[iatm][0] += np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax+slopeUp*((t-Tmm[0])%1.),XYZmax-slopeDn*((t-Tmm[1])%1.)) for t in glTau]).T
    966966        elif 'Block' in waveTypes[iatm]:
    967             nt = 1
     967            nx = 1
    968968            Tmm = Ax[iatm][0][:2]                       
    969969            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
    970             XmodZ[iatm][0] = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax,XYZmax) for t in glTau]).T                   
     970            XmodZ[iatm][0] += np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax,XYZmax) for t in glTau]).T                   
    971971        tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
    972         XmodA[iatm] = Ax[iatm,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X 32
    973         XmodB[iatm] = Bx[iatm,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
     972        XmodA[iatm][nx:] = Ax[iatm,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X 32
     973        XmodB[iatm][nx:] = Bx[iatm,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
    974974    Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X 3 X 32; sum waves
    975975    Xmod = np.swapaxes(Xmod,1,2)
     
    981981    else:
    982982        Umod = 1.0
     983#    GSASIIpath.IPyBreak()
    983984    return nWaves,Fmod,Xmod,Umod,glTau,glWt
    984985       
     
    10231024    Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
    10241025    Bf = np.array(FSSdata[1]).T    #cos frac mods...
    1025     Ax = np.array(XSSdata[:3]).T   #...cos pos mods
     1026    Ax = np.array(XSSdata[:3]).T   #...cos pos mods x waves x atoms
    10261027    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
    10271028    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
    10281029    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
    10291030    nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]]
    1030     nf = 0
    10311031    nx = 0
    1032 #    if 'Fourier' in waveTypes:
    1033 #        nf = 0
    1034 #        nx = 0
    1035 #        XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32))
    1036 #        FmodZ = np.zeros((Af.shape[0],Af.shape[1],32))
    1037 #        if 'Crenel' in waveTypes:
    1038 #            nC = np.where('Crenel' in waveTypes)
    1039 #            nf = 1
    1040 #            #FmodZ = 0   replace
    1041 #    else:
    1042 #        nx = 1
    1043 #        if 'Sawtooth' in waveTypes:
    1044 #            nS = np.where('Sawtooth' in waveTypes)
    1045 #            #XmodZ = 0   replace
    1046     tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
    1047     StauX = np.ones_like(Ax)[:,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X 32
    1048     CtauX = np.ones_like(Bx)[:,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
     1032    StauX = np.zeros((Ax.shape[0],Ax.shape[1],3,32))    #atoms x waves x 3 x 32
     1033    CtauX = np.zeros((Ax.shape[0],Ax.shape[1],3,32))
     1034    ZtauXx = np.zeros((Ax.shape[0],3,32))               #atoms x XYZmax x 32
     1035    ZtauXt = np.zeros((Ax.shape[0],2,32))               #atoms x Tminmax x 32
     1036    for iatm in range(Ax.shape[0]):
     1037        nx = 0
     1038        if 'ZigZag' in waveTypes[iatm]:
     1039            nx = 1
     1040#            Tmm = Ax[iatm][0][:2]                       
     1041#            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
     1042#            DT = Tmm[1]-Tmm[0]
     1043#            slopeUp = 2.*XYZmax/DT
     1044#            slopeDn = 2.*XYZmax/(1.-DT)
     1045#            XmodZ[iatm][0] += np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax+slopeUp*((t-Tmm[0])%1.),XYZmax-slopeDn*((t-Tmm[1])%1.)) for t in glTau]).T
     1046        elif 'Block' in waveTypes[iatm]:
     1047            nx = 1
     1048            Tmm = Ax[iatm][0][:2]                       
     1049            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
     1050            ZtauXt[iatm] = np.ones(2)[:,nxs]
     1051            ZtauXx[iatm] += np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-np.ones(3),np.ones(3)) for t in glTau]).T                   
     1052        tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
     1053        StauX[iatm][nx:] = np.ones_like(Ax)[iatm,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X 32
     1054        CtauX[iatm][nx:] = np.ones_like(Bx)[iatm,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
     1055#    GSASIIpath.IPyBreak()
    10491056    if nWaves[0]:
    10501057        tauF = np.arange(1.,nWaves[0]+1-nf)[:,nxs]*glTau  #Fwaves x 32
    1051         StauF = np.ones_like(Af)[:,nf:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf
    1052         CtauF = np.ones_like(Bf)[:,nf:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf
     1058        StauF = np.ones_like(Af)[:,:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf
     1059        CtauF = np.ones_like(Bf)[:,:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf
    10531060    else:
    10541061        StauF = 1.0
     
    10681075        UmodA = 0.
    10691076        UmodB = 0.
    1070     return waveShapes,[StauF,CtauF],[StauX,CtauX],[StauU,CtauU],UmodA+UmodB
     1077    return waveShapes,[StauF,CtauF],[StauX,CtauX,ZtauXt,ZtauXx],[StauU,CtauU],UmodA+UmodB
    10711078   
    10721079def ModulationDerv(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt):
     
    11161123    dGdMxSb = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   
    11171124    dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1)
     1125# ZigZag/Block waves
     1126    dHdXZt = twopi*np.ones(HP.shape[0])[:,nxs,nxs,nxs]*np.swapaxes(SCtauX[2],-1,-2)[nxs,:,:,:]          #??ops x atoms x 32 x 3(ZigZag/Block Tminmax)
     1127    dHdXZx = twopi*HP[:,nxs,nxs,:]*np.swapaxes(SCtauX[3],-1,-2)[nxs,:,:,:]          #ops x atoms x 32 x 3(ZigZag/Block XYZmax)
     1128    dGdMzCt = -np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
     1129    dGdMzCx = -np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.sin(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
     1130    dGdMzC = np.concatenate((dGdMzCt,dGdMzCx),axis=-1)
     1131    dGdMzSt = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
     1132    dGdMzSx = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
     1133    dGdMzS = np.concatenate((dGdMzSt,dGdMzSx),axis=-1)
    11181134# ops x atoms x waves x 2xyz - imaginary part - good
    1119     return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS]
    1120    
    1121 def posFourier(tau,psin,pcos,smul):
    1122     A = np.array([ps[:,np.newaxis]*np.sin(2*np.pi*(i+1)*tau) for i,ps in enumerate(psin)])  #*smul
     1135#    GSASIIpath.IPyBreak()
     1136    return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS]
     1137   
     1138def posFourier(tau,psin,pcos):
     1139    A = np.array([ps[:,np.newaxis]*np.sin(2*np.pi*(i+1)*tau) for i,ps in enumerate(psin)])
    11231140    B = np.array([pc[:,np.newaxis]*np.cos(2*np.pi*(i+1)*tau) for i,pc in enumerate(pcos)])
    11241141    return np.sum(A,axis=0)+np.sum(B,axis=0)
     
    11751192            sop,ssop,icent = G2spc.OpsfromStringOps(opr,SGData,SSGData)
    11761193            sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(tau,sop,ssop,atxyz)
    1177             smul = sdet         # n-fold vs m operator on wave
    11781194            tauT *= icent       #invert wave on -1
    11791195            wave = np.zeros(3)
     
    12061222                        ccof.append(spos[0][3:])
    12071223                if len(scof):
    1208                     wave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof),smul),axis=1)
     1224                    wave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1)
    12091225            if len(Sadp):
    12101226                scof = []
     
    12131229                    scof.append(sadp[0][:6])
    12141230                    ccof.append(sadp[0][6:])
    1215                 uwave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof),smul),axis=1)
     1231                uwave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1)
    12161232            if atom[cia] == 'A':                   
    12171233                X,U = G2spc.ApplyStringOps(opr,SGData,atxyz+wave,atuij+uwave)
  • trunk/GSASIIphsGUI.py

    r2064 r2070  
    24172417            waveSizer.Add(waveHead)
    24182418            if len(waveBlk):
    2419                 nFour = 0
     2419                nx = 0
    24202420                for iwave,wave in enumerate(waveBlk):
    2421                     if waveType == 'Fourier':
    2422                         nFour += 1
    24232421                    if not iwave:
    2424                         CSI = G2spc.GetSSfxuinel(waveType,nFour,xyz,SGData,SSGData)
     2422                        if waveType in ['ZigZag','Block']:
     2423                            nx = 1
     2424                        CSI = G2spc.GetSSfxuinel(waveType,1,xyz,SGData,SSGData)
    24252425                    else:
    2426                         CSI = G2spc.GetSSfxuinel('Fourier',nFour,xyz,SGData,SSGData)
     2426                        CSI = G2spc.GetSSfxuinel('Fourier',iwave+1-nx,xyz,SGData,SSGData)
    24272427                    waveName = 'Fourier'
    24282428                    if Stype == 'Sfrac':
  • trunk/GSASIIplot.py

    r2065 r2070  
    31213121                scof.append(spos[0][:3])
    31223122                ccof.append(spos[0][3:])
    3123         wave += G2mth.posFourier(tau,np.array(scof),np.array(ccof),1)
     3123        wave += G2mth.posFourier(tau,np.array(scof),np.array(ccof))
    31243124    if mapData['Flip']:
    31253125        Title = 'Charge flip'
     
    31503150    Plot.set_xlabel('t')
    31513151    Plot.set_ylabel(r'$\mathsf{\Delta}$%s'%(Ax))
    3152     Slab = np.hstack((slab,slab,slab))
     3152    Slab = np.hstack((slab,slab,slab))   
    31533153    acolor = mpl.cm.get_cmap('RdYlGn')
    31543154    if 'delt' in MapType:
  • trunk/GSASIIstrIO.py

    r2064 r2070  
    11681168                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
    11691169                        Waves = AtomSS[Stype]
     1170                        nx = 0
    11701171                        for iw,wave in enumerate(Waves):
    11711172                            if not iw:
    1172                                 CSI = G2spc.GetSSfxuinel(waveType,iw+1,at[cx:cx+3],SGData,SSGData)
     1173                                if waveType in ['ZigZag','Block']:
     1174                                    nx = 1
     1175                                CSI = G2spc.GetSSfxuinel(waveType,1,at[cx:cx+3],SGData,SSGData)
    11731176                            else:
    1174                                 CSI = G2spc.GetSSfxuinel('Fourier',iw+1,at[cx:cx+3],SGData,SSGData)
     1177                                CSI = G2spc.GetSSfxuinel('Fourier',iw+1-nx,at[cx:cx+3],SGData,SSGData)
    11751178                            uId,uCoef = CSI[Stype]
    11761179                            stiw = str(i)+':'+str(iw)
  • trunk/GSASIIstrMath.py

    r2067 r2070  
    11641164        dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
    11651165        dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
     1166        dFdGz = np.zeros((nRef,mSize,5))
    11661167        dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
    11671168    Flack = 1.0
     
    12101211        fotp = FPP*Tcorr            #ops x atoms
    12111212        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
    1212         dGdf,dGdx,dGdu = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
     1213        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
    12131214        # GfpuA is 2 x ops x atoms
    1214         # derivs are: ops x atoms x waves x 2,6,or 12 parms as [real,imag] parts
     1215        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
    12151216        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
    12161217        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
     
    12411242            dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
    12421243            dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
     1244            dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:],axis=1)
     1245            dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:],axis=1)
    12431246            dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
    12441247            dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
     
    12531256            dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
    12541257            dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
     1258            dfadGz = np.sum(fa[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]-fb[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
     1259            dfbdGz = np.sum(fb[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]+fa[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
    12551260            dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
    12561261            dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
     
    12981303                dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
    12991304                dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
     1305                dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
    13001306                dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
    13011307            else:   #these are good for no twin single crystals
     
    13081314                dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
    13091315                dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
     1316                dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1])      #array(nRef,natom,5)
    13101317                dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
    13111318        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
     
    13301337            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
    13311338            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
    1332         for j in range(XSSdata.shape[1]):       #loop over waves 'Tzero','Xslope','Yslope','Zslope'?
     1339        for j in range(XSSdata.shape[1]):       #loop over waves
    13331340            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j)] = dFdGx.T[0][j][i]
    13341341            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j)] = dFdGx.T[1][j][i]
     
    13371344            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j)] = dFdGx.T[4][j][i]
    13381345            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j)] = dFdGx.T[5][j][i]
     1346        dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]
     1347        dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[0][i]
     1348        dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[0][i]
     1349        dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[0][i]
     1350        dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[0][i]
    13391351        for j in range(USSdata.shape[1]):       #loop over waves
    13401352            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
Note: See TracChangeset for help on using the changeset viewer.