Changeset 1970 for trunk/GSASIImath.py


Ignore:
Timestamp:
Sep 2, 2015 4:01:41 PM (7 years ago)
Author:
vondreele
Message:

Add pygauleg to pypowder.for - recompile Win2.7 & Win2.7-64
work on SS structure factor

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1966 r1970  
    925925   
    926926def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast):
    927     import scipy.special as sp
    928     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
     927    import pypowder as pwd
     928   
     929    nxs = np.newaxis
     930    glTau,glWt = pwd.pygauleg(0.,1.,32)
     931   
     932    def expModInt(H,A,B):
     933        tau = np.arange(1.,A.shape[1]+1)[:,nxs]*glTau #waves x 32
     934        XmodA = np.swapaxes(A,1,2)[:,:,nxs]*np.sin(twopi*tau.T) #atoms X pos X 32 X waves
     935        XmodB = np.swapaxes(B,1,2)[:,:,nxs]*np.cos(twopi*tau.T)
     936        Xmod = np.sum(XmodA+XmodB,axis=-1)                    #atoms X pos X 32; sum waves
     937        #Xmod is ParSup in J2K - values match after Calcm2 with klic=-1
     938        D = H[:,3][:,nxs]*tau[0,nxs]
     939        HdotX = np.inner(np.swapaxes(Xmod,1,2),H[:,:3])+D.T         #atoms X 32 X ops
     940        sinHA = np.sum(np.sin(twopi*HdotX)*glWt[nxs,:,nxs],axis=1)
     941        cosHA = np.sum(np.cos(twopi*HdotX)*glWt[nxs,:,nxs],axis=1)
     942           
     943#        GSASIIpath.IPyBreak()     
     944        return cosHA.T,sinHA.T
    945945
    946946    Smult,TauT = SStauM             # both atoms x SGops
     
    952952    Ab = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
    953953    Bb = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
    954 #    GSASIIpath.IPyBreak()       
    955954    if Ax.ndim > 2:
    956         GpA = np.array(expModInt(SSUniq,Ax,Bx,Smult))
     955        GpA = np.array(expModInt(SSUniq,Ax,Bx))
    957956    else:
    958         GpA = np.array(expModInt(SSUniq,Ax[:,np.newaxis,:],Bx[:,np.newaxis,:],Smult))       
     957        GpA = np.array(expModInt(SSUniq,Ax[:,np.newaxis,:],Bx[:,np.newaxis,:]))       
    959958    return GpA             # SGops x atoms
    960959   
     
    10721071                drawatom[dcx:dcx+3] = X
    10731072    return drawAtoms,Fade
     1073   
     1074# gauleg.py Gauss Legendre numerical quadrature, x and w computation
     1075# integrate from a to b using n evaluations of the function f(x) 
     1076# usage: from gauleg import gaulegf         
     1077#        x,w = gaulegf( a, b, n)                               
     1078#        area = 0.0                                           
     1079#        for i in range(1,n+1):          #  yes, 1..n                   
     1080#          area += w[i]*f(x[i])                                   
     1081
     1082import math
     1083def gaulegf(a, b, n):
     1084  x = range(n+1) # x[0] unused
     1085  w = range(n+1) # w[0] unused
     1086  eps = 3.0E-14
     1087  m = (n+1)/2
     1088  xm = 0.5*(b+a)
     1089  xl = 0.5*(b-a)
     1090  for i in range(1,m+1):
     1091    z = math.cos(3.141592654*(i-0.25)/(n+0.5))
     1092    while True:
     1093      p1 = 1.0
     1094      p2 = 0.0
     1095      for j in range(1,n+1):
     1096        p3 = p2
     1097        p2 = p1
     1098        p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
     1099
     1100      pp = n*(z*p1-p2)/(z*z-1.0)
     1101      z1 = z
     1102      z = z1 - p1/pp
     1103      if abs(z-z1) <= eps:
     1104            break
     1105
     1106    x[i] = xm - xl*z
     1107    x[n+1-i] = xm + xl*z
     1108    w[i] = 2.0*xl/((1.0-z*z)*pp*pp)
     1109    w[n+1-i] = w[i]
     1110  return np.array(x), np.array(w)
     1111# end gaulegf
     1112   
    10741113   
    10751114def BessJn(nmax,x):
Note: See TracChangeset for help on using the changeset viewer.