Changeset 1310 for trunk/GSASIIsasd.py


Ignore:
Timestamp:
Apr 30, 2014 2:22:33 PM (7 years ago)
Author:
vondreele
Message:

Add sticky hard spheres to structure factor list

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIsasd.py

    r1309 r1310  
    519519    returns numpy array S(Q)
    520520    '''
    521     R,VolFr,epis,tau = args
    522 #//      Input (fitting) variables are:
    523 #       //radius = w[0]
    524 #       //volume fraction = w[1]
    525 #       //epsilon (perurbation param) = w[2]
    526 #       //tau (stickiness) = w[3]
    527 #       Variable rad,phi,eps,tau,eta
    528 #       Variable sig,aa,etam1,qa,qb,qc,radic
    529 #       Variable lam,lam2,test,mu,alpha,BetaVar
    530 #       Variable qv,kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq
    531 #       rad = w[0]
    532 #       phi = w[1]
    533 #       eps = w[2]
    534 #       tau = w[3]
    535 #       
    536 #       eta = phi/(1.0-eps)/(1.0-eps)/(1.0-eps)
    537 #       
    538 #       sig = 2.0 * rad
    539 #       aa = sig/(1.0 - eps)
    540 #       etam1 = 1.0 - eta
    541 #//C
    542 #//C  SOLVE QUADRATIC FOR LAMBDA
    543 #//C
    544 #       qa = eta/12.0
    545 #       qb = -1.0*(tau + eta/(etam1))
    546 #       qc = (1.0 + eta/2.0)/(etam1*etam1)
    547 #       radic = qb*qb - 4.0*qa*qc
    548 #       if(radic<0)
    549 #               if(x>0.01 && x<0.015)
    550 #                       Print "Lambda unphysical - both roots imaginary"
    551 #               endif
    552 #               return(-1)
    553 #       endif
    554 #//C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
    555 #       lam = (-1.0*qb-sqrt(radic))/(2.0*qa)
    556 #       lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa)
    557 #       if(lam2<lam)
    558 #               lam = lam2
    559 #       endif
    560 #       test = 1.0 + 2.0*eta
    561 #       mu = lam*eta*etam1
    562 #       if(mu>test)
    563 #               if(x>0.01 && x<0.015)
    564 #                Print "Lambda unphysical mu>test"
    565 #               endif
    566 #               return(-1)
    567 #       endif
    568 #       alpha = (1.0 + 2.0*eta - mu)/(etam1*etam1)
    569 #       BetaVar = (mu - 3.0*eta)/(2.0*etam1*etam1)
    570 #       
    571 #//C
    572 #//C   CALCULATE THE STRUCTURE FACTOR
    573 #//C
    574 #
    575 #       qv = x
    576 #       kk = qv*aa
    577 #       k2 = kk*kk
    578 #       k3 = kk*k2
    579 #       ds = sin(kk)
    580 #       dc = cos(kk)
    581 #       aq1 = ((ds - kk*dc)*alpha)/k3
    582 #       aq2 = (BetaVar*(1.0-dc))/k2
    583 #       aq3 = (lam*ds)/(12.0*kk)
    584 #       aq = 1.0 + 12.0*eta*(aq1+aq2-aq3)
    585 #//
    586 #       bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3)
    587 #       bq2 = BetaVar*(1.0/kk - ds/k2)
    588 #       bq3 = (lam/12.0)*((1.0 - dc)/kk)
    589 #       bq = 12.0*eta*(bq1+bq2-bq3)
    590 #//
    591 #       sq = 1.0/(aq*aq +bq*bq)
    592 #
    593 #       Return (sq)
    594     pass
    595    
     521    R,VolFr,epis,sticky = args
     522    eta = VolFr/(1.0-epis)/(1.0-epis)/(1.0-epis)       
     523    sig = 2.0 * R
     524    aa = sig/(1.0 - epis)
     525    etam1 = 1.0 - eta
     526#  SOLVE QUADRATIC FOR LAMBDA
     527    qa = eta/12.0
     528    qb = -1.0*(sticky + eta/(etam1))
     529    qc = (1.0 + eta/2.0)/(etam1*etam1)
     530    radic = qb*qb - 4.0*qa*qc
     531#   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
     532    lam1 = (-1.0*qb-np.sqrt(radic))/(2.0*qa)
     533    lam2 = (-1.0*qb+np.sqrt(radic))/(2.0*qa)
     534    lam = min(lam1,lam2)
     535    test = 1.0 + 2.0*eta
     536    mu = lam*eta*etam1
     537    alp = (1.0 + 2.0*eta - mu)/(etam1*etam1)
     538    bet = (mu - 3.0*eta)/(2.0*etam1*etam1)
     539#   CALCULATE THE STRUCTURE FACTOR<P></P>
     540    kk = Q*aa
     541    k2 = kk*kk
     542    k3 = kk*k2
     543    ds = np.sin(kk)
     544    dc = np.cos(kk)
     545    aq1 = ((ds - kk*dc)*alp)/k3
     546    aq2 = (bet*(1.0-dc))/k2
     547    aq3 = (lam*ds)/(12.0*kk)
     548    aq = 1.0 + 12.0*eta*(aq1+aq2-aq3)
     549
     550    bq1 = alp*(0.5/kk - ds/k2 + (1.0 - dc)/k3)
     551    bq2 = bet*(1.0/kk - ds/k2)
     552    bq3 = (lam/12.0)*((1.0 - dc)/kk)
     553    bq = 12.0*eta*(bq1+bq2-bq3)
     554    sq = 1.0/(aq*aq +bq*bq)
     555
     556    return sq
     557
    596558#def HayterPenfoldSF(Q,args): #big & ugly function - do later (if ever)
    597559#    pass
Note: See TracChangeset for help on using the changeset viewer.