Changeset 1966 for trunk/GSASIImath.py


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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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################################################################################
Note: See TracChangeset for help on using the changeset viewer.