Changeset 1310

Timestamp:

Apr 30, 2014 2:22:33 PM (9 years ago)

Author:

vondreele

Message:

Add sticky hard spheres to structure factor list

File:

• trunk/GSASIIsasd.py (modified) (1 diff)

trunk/GSASIIsasd.py

@@ -519 +519 @@
     returns numpy array S(Q)
     '''
-    R,VolFr,epis,tau = args
-#//      Input (fitting) variables are:
-#       //radius = w[0]
-#       //volume fraction = w[1]
-#       //epsilon (perurbation param) = w[2]
-#       //tau (stickiness) = w[3]
-#       Variable rad,phi,eps,tau,eta
-#       Variable sig,aa,etam1,qa,qb,qc,radic
-#       Variable lam,lam2,test,mu,alpha,BetaVar
-#       Variable qv,kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq
-#       rad = w[0]
-#       phi = w[1]
-#       eps = w[2]
-#       tau = w[3]
-#       
-#       eta = phi/(1.0-eps)/(1.0-eps)/(1.0-eps)
-#       
-#       sig = 2.0 * rad
-#       aa = sig/(1.0 - eps)
-#       etam1 = 1.0 - eta
-#//C
-#//C  SOLVE QUADRATIC FOR LAMBDA
-#//C
-#       qa = eta/12.0
-#       qb = -1.0*(tau + eta/(etam1))
-#       qc = (1.0 + eta/2.0)/(etam1*etam1)
-#       radic = qb*qb - 4.0*qa*qc
-#       if(radic<0)
-#               if(x>0.01 && x<0.015)
-#                       Print "Lambda unphysical - both roots imaginary"
-#               endif
-#               return(-1)
-#       endif
-#//C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
-#       lam = (-1.0*qb-sqrt(radic))/(2.0*qa)
-#       lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa)
-#       if(lam2<lam)
-#               lam = lam2
-#       endif
-#       test = 1.0 + 2.0*eta
-#       mu = lam*eta*etam1
-#       if(mu>test)
-#               if(x>0.01 && x<0.015)
-#                Print "Lambda unphysical mu>test"
-#               endif
-#               return(-1)
-#       endif
-#       alpha = (1.0 + 2.0*eta - mu)/(etam1*etam1)
-#       BetaVar = (mu - 3.0*eta)/(2.0*etam1*etam1)
-#       
-#//C
-#//C   CALCULATE THE STRUCTURE FACTOR
-#//C
-#
-#       qv = x
-#       kk = qv*aa
-#       k2 = kk*kk
-#       k3 = kk*k2
-#       ds = sin(kk)
-#       dc = cos(kk)
-#       aq1 = ((ds - kk*dc)*alpha)/k3
-#       aq2 = (BetaVar*(1.0-dc))/k2
-#       aq3 = (lam*ds)/(12.0*kk)
-#       aq = 1.0 + 12.0*eta*(aq1+aq2-aq3)
-#//
-#       bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3)
-#       bq2 = BetaVar*(1.0/kk - ds/k2)
-#       bq3 = (lam/12.0)*((1.0 - dc)/kk)
-#       bq = 12.0*eta*(bq1+bq2-bq3)
-#//
-#       sq = 1.0/(aq*aq +bq*bq)
-#
-#       Return (sq)
-    pass
-    
+    R,VolFr,epis,sticky = args
+    eta = VolFr/(1.0-epis)/(1.0-epis)/(1.0-epis)        
+    sig = 2.0 * R
+    aa = sig/(1.0 - epis)
+    etam1 = 1.0 - eta
+#  SOLVE QUADRATIC FOR LAMBDA
+    qa = eta/12.0
+    qb = -1.0*(sticky + eta/(etam1))
+    qc = (1.0 + eta/2.0)/(etam1*etam1)
+    radic = qb*qb - 4.0*qa*qc
+#   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
+    lam1 = (-1.0*qb-np.sqrt(radic))/(2.0*qa)
+    lam2 = (-1.0*qb+np.sqrt(radic))/(2.0*qa)
+    lam = min(lam1,lam2)
+    test = 1.0 + 2.0*eta
+    mu = lam*eta*etam1
+    alp = (1.0 + 2.0*eta - mu)/(etam1*etam1)
+    bet = (mu - 3.0*eta)/(2.0*etam1*etam1)
+#   CALCULATE THE STRUCTURE FACTOR<P></P>
+    kk = Q*aa
+    k2 = kk*kk
+    k3 = kk*k2
+    ds = np.sin(kk)
+    dc = np.cos(kk)
+    aq1 = ((ds - kk*dc)*alp)/k3
+    aq2 = (bet*(1.0-dc))/k2
+    aq3 = (lam*ds)/(12.0*kk)
+    aq = 1.0 + 12.0*eta*(aq1+aq2-aq3)
+
+    bq1 = alp*(0.5/kk - ds/k2 + (1.0 - dc)/k3)
+    bq2 = bet*(1.0/kk - ds/k2)
+    bq3 = (lam/12.0)*((1.0 - dc)/kk)
+    bq = 12.0*eta*(bq1+bq2-bq3)
+    sq = 1.0/(aq*aq +bq*bq)
+
+    return sq
+
 #def HayterPenfoldSF(Q,args): #big & ugly function - do later (if ever)
 #    pass