Changeset 2073


Ignore:
Timestamp:
Nov 30, 2015 8:33:57 AM (6 years ago)
Author:
vondreele
Message:

numeric derivatives for ZigZag/Block? waves

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r2072 r2073  
    962962            XmodZ[iatm][0] += posZigZag(glTau,Tmm,XYZmax).T
    963963        elif 'Block' in waveTypes[iatm]:
    964 
    965             nx = 1
    966964            Tmm = Ax[iatm][0][:2]                       
    967965            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
    968966            XmodZ[iatm][0] += posBlock(glTau,Tmm,XYZmax).T
    969967        tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
    970         XmodA[iatm][nx:] = Ax[iatm,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X 32
    971         XmodB[iatm][nx:] = Bx[iatm,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
     968        if nx:   
     969            XmodA[iatm][:-nx] = Ax[iatm,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X 32
     970            XmodB[iatm][:-nx] = Bx[iatm,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
     971        else:
     972            XmodA[iatm] = Ax[iatm,:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X 32
     973            XmodB[iatm] = Bx[iatm,:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
    972974    Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X 3 X 32; sum waves
    973975    Xmod = np.swapaxes(Xmod,1,2)
     
    10191021    '''
    10201022    glTau,glWt = pwd.pygauleg(0.,1.,32)         #get Gauss-Legendre intervals & weights
     1023    dT = 1/32.
     1024    dX = 0.0001
    10211025    waveShapes = [FSSdata.T.shape,XSSdata.T.shape,USSdata.T.shape]
    10221026    Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
     
    10361040        if 'ZigZag' in waveTypes[iatm]:
    10371041            nx = 1
    1038 #            Tmm = Ax[iatm][0][:2]                       
    1039 #            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
    1040 #            DT = Tmm[1]-Tmm[0]
    1041 #            slopeUp = 2.*XYZmax/DT
    1042 #            slopeDn = 2.*XYZmax/(1.-DT)
    1043 #            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
     1042            Tmm = Ax[iatm][0][:2]                       
     1043            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])           
     1044            for i in range(2):
     1045                Tmm[i] += dT
     1046                Zp = posZigZag(glTau,Tmm,XYZmax).T
     1047                Tmm[i] -= 2.*dT
     1048                Zm = posZigZag(glTau,Tmm,XYZmax).T
     1049                Tmm[i] += dT
     1050                ZtauXt[iatm][i] = (Zp-Zm)[i]/(2.*dT)
     1051            for i in range(3):
     1052                XYZmax[i] += dX
     1053                Zp = posZigZag(glTau,Tmm,XYZmax).T
     1054                XYZmax[i] -= 2.*dX
     1055                Zm = posZigZag(glTau,Tmm,XYZmax).T
     1056                XYZmax[i] += dX
     1057                ZtauXx[iatm][i] = (Zp-Zm)[i]/(2.*dX)
    10441058        elif 'Block' in waveTypes[iatm]:
    10451059            nx = 1
    10461060            Tmm = Ax[iatm][0][:2]                       
    1047             XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
    1048             ZtauXt[iatm] = np.ones(2)[:,nxs]
    1049             ZtauXx[iatm] += np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-np.ones(3),np.ones(3)) for t in glTau]).T                   
     1061            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])           
     1062            for i in range(2):
     1063                Tmm[i] += dT
     1064                Zp = posBlock(glTau,Tmm,XYZmax).T
     1065                Tmm[i] -= 2.*dT
     1066                Zm = posBlock(glTau,Tmm,XYZmax).T
     1067                Tmm[i] += dT
     1068                ZtauXt[iatm][i] = (Zp-Zm)[i]/(2.*dT)
     1069            for i in range(3):
     1070                XYZmax[i] += dX
     1071                Zp = posBlock(glTau,Tmm,XYZmax).T
     1072                XYZmax[i] -= 2.*dX
     1073                Zm = posBlock(glTau,Tmm,XYZmax).T
     1074                XYZmax[i] += dX
     1075                ZtauXx[iatm][i] = (Zp-Zm)[i]/(2.*dX)
    10501076        tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
    1051         StauX[iatm][nx:] = np.ones_like(Ax)[iatm,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X 32
    1052         CtauX[iatm][nx:] = np.ones_like(Bx)[iatm,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
     1077        if nx:   
     1078            StauX[iatm][:-nx] = np.ones_like(Ax)[iatm,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X 32
     1079            CtauX[iatm][:-nx] = np.ones_like(Bx)[iatm,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
     1080        else:
     1081            StauX[iatm] = np.ones_like(Ax)[iatm,:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X 32
     1082            CtauX[iatm] = np.ones_like(Bx)[iatm,:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
    10531083#    GSASIIpath.IPyBreak()
    10541084    if nWaves[0]:
     
    11091139        dGdMuSb = np.sum(part1[:,:,:,:,nxs]*dUdBu*np.sin(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    11101140        dGdMuS = np.concatenate((dGdMuSa,dGdMuSb),axis=-1)   #ops x atoms x waves x 12uij
    1111 #        #GSASIIpath.IPyBreak()
    11121141    else:
    11131142        HbH = np.ones_like(HdotX)
     
    11301159    dGdMzSx = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
    11311160    dGdMzS = np.concatenate((dGdMzSt,dGdMzSx),axis=-1)
     1161#    GSASIIpath.IPyBreak()
    11321162# ops x atoms x waves x 2xyz - imaginary part - good
    1133 #    GSASIIpath.IPyBreak()
    11341163    return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS]
    11351164   
Note: See TracChangeset for help on using the changeset viewer.