Changeset 2cefd79 in sasmodels for explore/beta/sasfit_compare.py


Ignore:
Timestamp:
Jul 5, 2018 11:01:54 AM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
757a3ff
Parents:
b1a0f3e
Message:

restructure beta approx tests to allow command line options

File:
1 edited

Legend:

Unmodified
Added
Removed

  • explore/beta/sasfit_compare.py

    r707cbdb r2cefd79  
    11from __future__ import division, print_function 
    22# Make sasmodels available on the path 
    3 import sys,os 
     3import sys, os 
    44BETA_DIR = os.path.dirname(os.path.realpath(__file__)) 
    55SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR)) 
    66sys.path.insert(0, SASMODELS_DIR) 
    7 import os 
     7 
     8from collections import namedtuple 
     9 
    810from matplotlib import pyplot as plt 
    911import numpy as np 
    … …  
    1416from scipy.special import gammaln  # type: ignore 
    1517 
    16 # THE FOLLOWING 7 BOOLEANS TAILOR WHICH GRAPHS ARE PRINTED WHEN THE PROGRAM RUNS 
    17 #RICHARD, YOU WILL WANT SASVIEW, SASFIT, SCHULZ, ELLIPSOID, AND SPHERE ALL TRUE. 
    18 ELLIPSOID = False 
    19 SPHERE = True 
    20  
    21 GAUSSIAN = True 
    22 SCHULZ = False 
    23  
    24 SASVIEW=True 
    25 SASFIT=True 
    26 YUN = True 
    27  
    28 def data_file(name): 
    29     return os.path.join(BETA_DIR + '\\data_files', name) 
     18Theory = namedtuple('Theory', 'Q F1 F2 P S I Seff Ibeta') 
     19Theory.__new__.__defaults__ = (None,) * len(Theory._fields) 
    3020 
    3121#Used to calculate F(q) for the cylinder, sphere, ellipsoid models 
    … …  
    4737    with np.errstate(all='ignore'): 
    4838        # GSL bessel_j1 taylor expansion 
    49         index = (x < 0.25)   
     39        index = (x < 0.25) 
    5040        y = x[index]**2 
    5141        c1 = -1.0/10.0 
    … …  
    7565 
    7666#gives the hardsphere structure factor that sasview uses 
    77 def hardsphere_simple(q, radius_effective, volfraction):  
    78     CUTOFFHS=0.05  
     67def _hardsphere_simple(q, radius_effective, volfraction): 
     68    CUTOFFHS=0.05 
    7969    if fabs(radius_effective) < 1.E-12: 
    8070        HARDSPH=1.0 
    … …  
    9686        FF = 8.0*A +6.0*B + 4.0*G + ( -0.8*A -B/1.5 -0.5*G +(A/35. +0.0125*B +0.02*G)*X2)*X2 
    9787        HARDSPH= 1./(1. + volfraction*FF ) 
    98         return HARDSPH    
     88        return HARDSPH 
    9989    X4=X2*X2 
    10090    S, C = sin(X), cos(X) 
    101     FF=  (( G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )/X2 + B*(2.*X*S -(X2-2.)*C -2.) )/X + A*(S-X*C))/X ; 
    102     HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ); 
     91    FF=  (( G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )/X2 + B*(2.*X*S -(X2-2.)*C -2.) )/X + A*(S-X*C))/X 
     92    HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ) 
    10393    return HARDSPH 
     94 
     95def hardsphere_simple(q, radius_effective, volfraction): 
     96    SQ = [_hardsphere_simple(qk, radius_effective, volfraction) for qk in q] 
     97    return np.array(SQ) 
    10498 
    10599#Used in gaussian quadrature for polydispersity 
    106100#returns values and the probability of those values based on gaussian distribution 
    107 def gaussian_distribution(center, sigma,lb,ub): 
    108     #3 standard deviations covers approx. 99.7%  
     101N_GAUSS = 35 
     102NSIGMA_GAUSS = 3 
     103def gaussian_distribution(center, sigma, lb, ub): 
     104    #3 standard deviations covers approx. 99.7% 
    109105    if sigma != 0: 
    110         nsigmas=3 
    111         x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=35) 
     106        nsigmas = NSIGMA_GAUSS 
     107        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_GAUSS) 
    112108        x= x[(x >= lb) & (x <= ub)] 
    113109        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma)) 
    … …  
    116112        return np.array([center]), np.array([1]) 
    117113 
     114N_SCHULZ = 80 
     115NSIGMA_SCHULZ = 8 
    118116def schulz_distribution(center, sigma, lb, ub): 
    119117    if sigma != 0: 
    120         nsigmas=8 
    121         x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=80) 
     118        nsigmas = NSIGMA_SCHULZ 
     119        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_SCHULZ) 
    122120        x= x[(x >= lb) & (x <= ub)] 
    123121        R = x/center 
    … …  
    160158#SQ is monodisperse approach for structure factor 
    161159#SQ_EFF is the effective structure factor from beta approx 
    162 def ellipsoid_Theta(q, radius_polar, radius_equatorial, sld, sld_solvent,volfraction=0,effective_radius=0): 
     160def ellipsoid_theta(q, radius_polar, radius_equatorial, sld, sld_solvent, 
     161                    volfraction=0, radius_effective=None): 
    163162    #creates values z and corresponding probabilities w from legendre-gauss quadrature 
     163    volume = ellipsoid_volume(radius_polar, radius_equatorial) 
    164164    z, w = leggauss(76) 
    165165    F1 = np.zeros_like(q) 
    166166    F2 = np.zeros_like(q) 
    167167    #use a u subsition(u=cos) and then u=(z+1)/2 to change integration from 
    168     #0->2pi with respect to alpha to -1->1 with respect to z, allowing us to use  
     168    #0->2pi with respect to alpha to -1->1 with respect to z, allowing us to use 
    169169    #legendre-gauss quadrature 
    170170    for k, qk in enumerate(q): 
    171171        r = sqrt(radius_equatorial**2*(1-((z+1)/2)**2)+radius_polar**2*((z+1)/2)**2) 
    172         F2i = ((sld-sld_solvent)*ellipsoid_volume(radius_polar,radius_equatorial)*sas_3j1x_x(qk*r))**2 
     172        F2i = ((sld-sld_solvent)*volume*sas_3j1x_x(qk*r))**2 
    173173        F2[k] = np.sum(w*F2i) 
    174         F1i = (sld-sld_solvent)*ellipsoid_volume(radius_polar,radius_equatorial)*sas_3j1x_x(qk*r) 
     174        F1i = (sld-sld_solvent)*volume*sas_3j1x_x(qk*r) 
    175175        F1[k] = np.sum(w*F1i) 
    176176    #the 1/2 comes from the change of variables mentioned above 
    177177    F2 = F2/2.0 
    178178    F1 = F1/2.0 
    179     if effective_radius == 0: 
    180         effective_radius = ER_ellipsoid(radius_polar,radius_equatorial) 
    181     else: 
    182         effective_radius = effective_radius 
    183     SQ = np.array([hardsphere_simple(qk, effective_radius, volfraction) for qk in q])    
    184     SQ_EFF = 1 + F1**2/F2*(SQ - 1)  
    185     IQM = 1e-4/ellipsoid_volume(radius_polar,radius_equatorial)*F2 
     179    if radius_effective is None: 
     180        radius_effective = ER_ellipsoid(radius_polar,radius_equatorial) 
     181    SQ = hardsphere_simple(q, radius_effective, volfraction) 
     182    SQ_EFF = 1 + F1**2/F2*(SQ - 1) 
     183    IQM = 1e-4*F2/volume 
    186184    IQSM = volfraction*IQM*SQ 
    187185    IQBM = volfraction*IQM*SQ_EFF 
    188     return F1, F2, IQM, IQSM, IQBM, SQ, SQ_EFF 
    189  
    190 #IQD is I(q) polydispursed, IQSD is I(q)S(q) polydispursed, etc.  
     186    return Theory(Q=q, F1=F1, F2=F2, P=IQM, S=SQ, I=IQSM, Seff=SQ_EFF, Ibeta=IQBM) 
     187 
     188#IQD is I(q) polydispursed, IQSD is I(q)S(q) polydispursed, etc. 
    191189#IQBD HAS NOT BEEN CROSS CHECKED AT ALL 
    192 def ellipsoid_pe(q, Rp, Re, sld, sld_solvent, volfraction=0,effective_radius=0,radius_polar_pd=0.1,radius_equatorial_pd=0.1,distribution='gaussian'): 
    193     if distribution == 'gaussian': 
    194         Rp_val, Rp_prob = gaussian_distribution(Rp, radius_polar_pd*Rp, 0, inf) 
    195         Re_val, Re_prob = gaussian_distribution(Re, radius_equatorial_pd*Re, 0, inf) 
    196     elif distribution == 'schulz': 
    197         Rp_val, Rp_prob = schulz_distribution(Rp, radius_polar_pd*Rp, 0, inf) 
    198         Re_val, Re_prob = schulz_distribution(Re, radius_equatorial_pd*Re, 0, inf) 
     190def ellipsoid_pe(q, radius_polar, radius_equatorial, sld, sld_solvent, 
     191                 radius_polar_pd=0.1, radius_equatorial_pd=0.1, 
     192                 radius_polar_pd_type='gaussian', 
     193                 radius_equatorial_pd_type='gaussian', 
     194                 volfraction=0, radius_effective=None, 
     195                 background=0, scale=1, 
     196                 norm='sasview'): 
     197    if norm not in ['sasview', 'sasfit', 'yun']: 
     198        raise TypeError("unknown norm "+norm) 
     199    if radius_polar_pd_type == 'gaussian': 
     200        Rp_val, Rp_prob = gaussian_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf) 
     201    elif radius_polar_pd_type == 'schulz': 
     202        Rp_val, Rp_prob = schulz_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf) 
     203    if radius_equatorial_pd_type == 'gaussian': 
     204        Re_val, Re_prob = gaussian_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf) 
     205    elif radius_equatorial_pd_type == 'schulz': 
     206        Re_val, Re_prob = schulz_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf) 
    199207    Normalization = 0 
    200     total_weight = 0 
    201     PQ = np.zeros_like(q) 
    202     F12,F21 = np.zeros_like(q), np.zeros_like(q) 
    203     radius_eff = 0 
     208    F1,F2 = np.zeros_like(q), np.zeros_like(q) 
     209    radius_eff = total_weight = 0 
    204210    for k, Rpk in enumerate(Rp_val): 
    205211        for i, Rei in enumerate(Re_val): 
    206             F1i, F2i, PQi, IQSM, IQBM, SQ, SQ_EFF = ellipsoid_Theta(q,Rpk,Rei,sld,sld_solvent) 
     212            theory = ellipsoid_theta(q,Rpk,Rei,sld,sld_solvent) 
     213            volume = ellipsoid_volume(Rpk, Rei) 
     214            if norm == 'sasfit': 
     215                Normalization += Rp_prob[k]*Re_prob[i] 
     216            elif norm == 'sasview' or norm == 'yun': 
     217                Normalization += Rp_prob[k]*Re_prob[i]*volume 
     218            F1 += theory.F1*Rp_prob[k]*Re_prob[i] 
     219            F2 += theory.F2*Rp_prob[k]*Re_prob[i] 
    207220            radius_eff += Rp_prob[k]*Re_prob[i]*ER_ellipsoid(Rpk,Rei) 
    208             total_weight +=  Rp_prob[k]*Re_prob[i] 
    209             Normalization += Rp_prob[k]*Re_prob[i]*ellipsoid_volume(Rpk, Rei) 
    210             PQ += PQi*Rp_prob[k]*Re_prob[i]*ellipsoid_volume(Rpk, Rei) 
    211             F12 += F1i*Rp_prob[k]*Re_prob[i]*ellipsoid_volume(Rpk, Rei) 
    212             F21 += F2i*Rp_prob[k]*Re_prob[i]*ellipsoid_volume(Rpk, Rei) 
    213     F12 = (F12/Normalization)**2 
    214     F21 = F21/Normalization 
    215     if effective_radius == 0: 
    216         effective_radius = radius_eff/total_weight 
    217     else: 
    218         effective_radius = effective_radius 
    219     SQ = np.array([hardsphere_simple(qk, effective_radius, volfraction) for qk in q])    
    220     SQ_EFF = 1 + F12/F21*(SQ - 1)      
    221     IQD = PQ/Normalization 
    222     IQSD = volfraction*IQD*SQ 
    223     IQBD = volfraction*IQD*SQ_EFF 
    224     return IQD, IQSD, IQBD, SQ, SQ_EFF 
     221            total_weight += Rp_prob[k]*Re_prob[i] 
     222    F1 = F1/Normalization 
     223    F2 = F2/Normalization 
     224    if radius_effective is None: 
     225        radius_effective = radius_eff/total_weight 
     226    SQ = hardsphere_simple(q, radius_effective, volfraction) 
     227    SQ_EFF = 1 + F1**2/F2*(SQ - 1) 
     228    volume = ellipsoid_volume(radius_polar, radius_equatorial) 
     229    if norm == 'sasfit': 
     230        IQD = F2 
     231        IQSD = IQD*SQ 
     232        IQBD = IQD*SQ_EFF 
     233    elif norm == 'sasview': 
     234        IQD = F2*1e-4*volfraction 
     235        IQSD = IQD*SQ 
     236        IQBD = IQD*SQ_EFF 
     237    elif norm == 'yun': 
     238        SQ_EFF = 1 + Normalization*F1**2/F2*(SQ - 1) 
     239        F2 = F2/volume 
     240        IQD = F2 
     241        IQSD = IQD*SQ 
     242        IQBD = IQD*SQ_EFF 
     243    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD) 
    225244 
    226245#polydispersity for sphere 
    227 def sphere_r(q,radius,sld,sld_solvent,volfraction=0,radius_pd=0.1,distribution='gaussian',norm_type='volfraction'): 
    228     if distribution == 'gaussian':     
     246def sphere_r(q,radius,sld,sld_solvent, 
     247             radius_pd=0.1, radius_pd_type='gaussian', 
     248             volfraction=0, radius_effective=None, 
     249             background=0, scale=1, 
     250             norm='sasview'): 
     251    if norm not in ['sasview', 'sasfit', 'yun']: 
     252        raise TypeError("unknown norm "+norm) 
     253    if radius_pd_type == 'gaussian': 
    229254        radius_val, radius_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf) 
    230     elif distribution == 'schulz': 
     255    elif radius_pd_type == 'schulz': 
    231256        radius_val, radius_prob = schulz_distribution(radius, radius_pd*radius, 0, inf) 
    232257    Normalization=0 
     258    F1 = np.zeros_like(q) 
    233259    F2 = np.zeros_like(q) 
    234     F1 = np.zeros_like(q) 
    235     if norm_type == 'numdensity': 
    236         for k, rk in enumerate(radius_val): 
    237             volume = 4./3.*pi*rk**3 
     260    for k, rk in enumerate(radius_val): 
     261        volume = 4./3.*pi*rk**3 
     262        if norm == 'sasfit': 
     263            Normalization += radius_prob[k] 
     264        elif norm == 'sasview' or norm == 'yun': 
    238265            Normalization += radius_prob[k]*volume 
    239             F2i = ((sld-sld_solvent)*volume*sas_3j1x_x(q*rk))**2/volume 
    240             F1i = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk)/volume 
    241              
    242             F2 += radius_prob[k]*volume*F2i 
    243             F1 += radius_prob[k]*volume*F1i 
    244         F21 = F2/Normalization 
    245         F12 = (F1/Normalization)**2 
    246         SQ = np.array([hardsphere_simple(qk, radius, volfraction) for qk in q]) 
    247         SQ_EFF = 1 + F12/F21*(SQ - 1)  
    248         IQD = 1e-4*F21 
    249         IQSD = volfraction*IQD*SQ 
    250         IQBD = volfraction*IQD*SQ_EFF 
    251     elif norm_type == 'volfraction': 
    252         for k, rk in enumerate(radius_val): 
    253             volume = 4./3.*pi*rk**3 
    254             Normalization += radius_prob[k] 
    255             F2i = ((sld-sld_solvent)*volume*sas_3j1x_x(q*rk))**2 
    256             F1i = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk) 
    257             F2 += radius_prob[k]*F2i 
    258             F1 += radius_prob[k]*F1i 
    259      
    260         F21 = F2/Normalization 
    261         F12 = (F1/Normalization)**2 
    262         SQ = np.array([hardsphere_simple(qk, radius, volfraction) for qk in q]) 
    263         SQ_EFF = 1 + F12/F21*(SQ - 1)  
    264         IQD = 1e-4/(4./3.*pi*radius**3)*F21 
    265         IQSD = volfraction*IQD*SQ 
    266         IQBD = volfraction*IQD*SQ_EFF 
    267     return IQD, IQSD, IQBD, SQ, SQ_EFF 
     266        F2k = ((sld-sld_solvent)*volume*sas_3j1x_x(q*rk))**2 
     267        F1k = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk) 
     268        F2 += radius_prob[k]*F2k 
     269        F1 += radius_prob[k]*F1k 
     270 
     271    F2 = F2/Normalization 
     272    F1 = F1/Normalization 
     273    if radius_effective is None: 
     274        radius_effective = radius 
     275    SQ = hardsphere_simple(q, radius_effective, volfraction) 
     276    SQ_EFF = 1 + F1**2/F2*(SQ - 1) 
     277    volume = 4./3.*pi*radius**3 
     278    if norm == 'sasfit': 
     279        IQD = F2 
     280        IQSD = IQD*SQ 
     281        IQBD = IQD*SQ_EFF 
     282    elif norm == 'sasview': 
     283        IQD = F2*1e-4*volfraction 
     284        IQSD = IQD*SQ 
     285        IQBD = IQD*SQ_EFF 
     286    elif norm == 'yun': 
     287        SQ_EFF = 1 + Normalization*F1**2/F2*(SQ - 1) 
     288        F2 = F2/volume 
     289        IQD = F2 
     290        IQSD = IQD*SQ 
     291        IQBD = IQD*SQ_EFF 
     292    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD) 
    268293 
    269294############################################################################### 
    … …  
    276301############################################################################### 
    277302############################################################################### 
    278 # SASVIEW     
    279 if (ELLIPSOID == True) & (GAUSSIAN == True) & (SASVIEW == True): 
     303 
     304def popn(d, keys): 
     305    """ 
     306    Splits a dict into two, with any key of *d* which is in *keys* removed 
     307    from *d* and added to *b*. Returns *b*. 
     308    """ 
     309    b = {} 
     310    for k in keys: 
     311        try: 
     312            b[k] = d.pop(k) 
     313        except KeyError: 
     314            pass 
     315    return b 
     316 
     317def sasmodels_theory(q, Pname, **pars): 
     318    """ 
     319    Call sasmodels to compute the model with and without a hard sphere 
     320    structure factor. 
     321    """ 
     322    #mono_pars = {k: (0 if k.endswith('_pd') else v) for k, v in pars.items()} 
     323    Ppars = pars.copy() 
     324    Spars = popn(Ppars, ['radius_effective', 'volfraction']) 
     325    Ipars = pars.copy() 
     326 
     327    # Autofill npts and nsigmas for the given pd type 
     328    for k, v in pars.items(): 
     329        if k.endswith("_pd_type"): 
     330            if v == "gaussian": 
     331                n, nsigmas = N_GAUSS, NSIGMA_GAUSS 
     332            elif v == "schulz": 
     333                n, nsigmas = N_SCHULZ, NSIGMA_SCHULZ 
     334            Ppars.setdefault(k.replace("_pd_type", "_pd_n"), n) 
     335            Ppars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas) 
     336            Ipars.setdefault(k.replace("_pd_type", "_pd_n"), n) 
     337            Ipars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas) 
     338 
     339    #Ppars['scale'] = Spars.get('volfraction', 1) 
     340    P = build_model(Pname, q) 
     341    S = build_model("hardsphere", q) 
     342    I = build_model(Pname+"@hardsphere", q) 
     343    Pq = P(**Ppars)*pars.get('volfraction', 1) 
     344    #Sq = S(**Spars) 
     345    Iq = I(**Ipars) 
     346    #Iq = Pq*Sq*pars.get('volfraction', 1) 
     347    Sq = Iq/Pq 
     348    return Theory(Q=q, F1=None, F2=None, P=Pq, S=Sq, I=Iq, Seff=None, Ibeta=None) 
     349 
     350def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'): 
     351    """ 
     352    Plot fields in common between target and actual, along with relative error. 
     353    """ 
     354    available = [s for s in fields.split() 
     355                 if getattr(target, s) is not None and getattr(actual, s) is not None] 
     356    rows = len(available) 
     357    for row, field in enumerate(available): 
     358        Q = target.Q 
     359        I1, I2 = getattr(target, field), getattr(actual, field) 
     360        plt.subplot(rows, 2, 2*row+1) 
     361        plt.loglog(Q, abs(I1), label="target "+field) 
     362        plt.loglog(Q, abs(I2), label="value "+field) 
     363        #plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
     364        plt.legend(loc='lower left') 
     365        plt.subplot(rows, 2, 2*row+2) 
     366        #plt.semilogx(Q, I2/I1 - 1, label="relative error") 
     367        plt.semilogx(Q, I1/I2 - 1, label="relative error") 
     368    plt.tight_layout() 
     369    plt.suptitle(title) 
     370    plt.show() 
     371 
     372def data_file(name): 
     373    return os.path.join(BETA_DIR, 'data_files', name) 
     374 
     375def load_sasfit(path): 
     376    data = np.loadtxt(path, dtype=str, delimiter=';').T 
     377    data = np.vstack((map(float, v) for v in data[0:2])) 
     378    return data 
     379 
     380COMPARISON = {}  # Type: Dict[(str,str,str)] -> Callable[(), None] 
     381 
     382def compare_sasview_sphere(pd_type='schulz'): 
     383    q = np.logspace(-5, 0, 250) 
     384    model = 'sphere' 
     385    pars = dict( 
     386        radius=20,sld=4,sld_solvent=1, 
     387        background=0, 
     388        radius_pd=.1, radius_pd_type=pd_type, 
     389        volfraction=0.15, 
     390        #radius_effective=12.59921049894873,  # equivalent average sphere radius 
     391        ) 
     392    target = sasmodels_theory(q, model, **pars) 
     393    actual = sphere_r(q, norm='sasview', **pars) 
     394    title = " ".join(("sasmodels", model, pd_type)) 
     395    compare(title, target, actual) 
     396COMPARISON[('sasview','sphere','gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian') 
     397COMPARISON[('sasview','sphere','schulz')] = lambda: compare_sasview_sphere(pd_type='schulz') 
     398 
     399def compare_sasview_ellipsoid(pd_type='gaussian'): 
    280400    q = np.logspace(-5, 0, 50) 
    281     volfraction = 0.2 
    282     model = build_model("ellipsoid",q) 
    283     Sq = model(radius_polar=20,radius_equatorial=400,sld=4,sld_solvent=1,\ 
    284                background=0,radius_polar_pd=.1,radius_polar_pd_n=35,radius_equatorial_pd=.1,radius_equatorial_pd_n=35) 
    285     IQD, IQSD, IQBD, SQ, SQ_EFF = ellipsoid_pe(q,20,400,4,1,volfraction) 
    286     Sq2 = model(radius_polar=20,radius_equatorial=10,sld=2,sld_solvent=1,background=0) 
    287     IQM = ellipsoid_Theta(q,20,10,2,1)[2] 
    288     model2 = build_model("ellipsoid@hardsphere", q) 
    289     Sq3 = model2(radius_polar=20,radius_equatorial=400,sld=4,sld_solvent=1,background=0,radius_polar_pd=.1,\ 
    290                  radius_polar_pd_n=35,radius_equatorial_pd=.1,radius_equatorial_pd_n=35) 
    291     Sqp = model(radius_polar=20,radius_equatorial=10,sld=4,sld_solvent=1,\ 
    292                    background=0,radius_polar_pd=0.1,radius_polar_pd_n=35,radius_equatorial_pd=0,radius_equatorial_pd_n=35) 
    293     IQDp = ellipsoid_pe(q,20,10,4,1,radius_polar_pd=0.1,radius_equatorial_pd=0)[0]   
    294     Sqp2 = model(radius_polar=20,radius_equatorial=10,sld=4,sld_solvent=1,\ 
    295                    background=0,radius_polar_pd=0,radius_polar_pd_n=35,radius_equatorial_pd=0.1,radius_equatorial_pd_n=35) 
    296     IQDe = ellipsoid_pe(q,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.1)[0] 
    297     print('\n\n\n\n ELLIPSOID COMPARISON WITH SASVIEW WITHOUT V') 
    298     plt.subplot(323) 
    299     plt.loglog(q, Sq,'r') 
    300     plt.loglog(q,IQD,'b') 
    301     plt.title('IQD') 
    302     plt.subplot(324) 
    303     plt.semilogx(q,Sq/IQD-1) 
    304     plt.subplot(321) 
    305     plt.title('IQM') 
    306     plt.loglog(q, Sq2,'r') 
    307     plt.loglog(q,IQM,'b') 
    308     plt.subplot(322) 
    309     plt.semilogx(q,IQM/Sq2-1) 
    310     plt.subplot(325) 
    311     plt.title('IQSD') 
    312     plt.loglog(q, Sq3, '-b') 
    313     plt.loglog(q, IQSD, '-r') 
    314     plt.subplot(326) 
    315     plt.semilogx(q, Sq3/IQSD-1, 'b') 
    316     plt.tight_layout() 
    317     plt.show() 
    318     plt.subplot(221) 
    319     plt.loglog(q,IQDp,'r') 
    320     plt.loglog(q,Sqp,'b') 
    321     plt.title('IQD Polar') 
    322     plt.subplot(222) 
    323     plt.plot(q,IQDp/Sqp-1) 
    324     plt.subplot(223) 
    325     plt.loglog(q,IQDe,'r') 
    326     plt.loglog(q,Sqp2,'b') 
    327     plt.title('IQD Equatorial') 
    328     plt.subplot(224) 
    329     plt.plot(q,IQDe/Sqp2-1) 
    330     plt.tight_layout() 
    331     plt.show() 
    332 #comparing monodisperse Beta approximation to version without 
    333     q=np.linspace(1e-5,1) 
    334     F1, F2, IQM, IQSM, IQBM, SQ, SQ_EFF = ellipsoid_Theta(q,20,10,2,1,0.15,12.59921049894873) 
    335     plt.subplot(321) 
    336     plt.loglog(q,IQBM,'g',label='IQBM') 
    337     plt.loglog(q,IQSM,'r',label= 'IQSM') 
    338     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    339     plt.subplot(323) 
    340     plt.plot(q,IQBM,'g',label='IQBM') 
    341     plt.plot(q,IQSM,'r',label= 'IQSM') 
    342     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    343     plt.subplot(325) 
    344     plt.plot(q,IQBM/IQSM,'r',label= 'IQBM/IQSM') 
    345     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    346     plt.tight_layout() 
    347     plt.show() 
    348 #comparing polydispersed Beta approximation to version without 
    349     IQD, IQSD, IQBD, SQ, SQ_EFF = ellipsoid_pe(q, 20,10, 2,1, 0.15,12.59921049894873) 
    350     plt.subplot(321) 
    351     plt.loglog(q, SQ) 
    352     plt.loglog(q, SQ_EFF,label='SQ, SQ_EFF') 
    353     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    354     plt.subplot(323) 
    355     plt.loglog(q,IQBD) 
    356     plt.loglog(q,IQSD,label='IQBD,IQSD') 
    357     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    358     plt.subplot(325) 
    359     plt.plot(q,IQBD) 
    360     plt.plot(q,IQSD,label='IQBD,IQSD') 
    361     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    362     plt.tight_layout() 
    363     plt.show() 
    364  
    365 # YUN 
    366 if (ELLIPSOID == True) & (GAUSSIAN == True) & (YUN == True): 
    367     #Yun's data for Beta Approximation 
    368     data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T    
    369     print('\n\n ELLIPSOID COMPARISON WITH YUN WITHOUT V') 
    370     Q=data[0] 
    371     F1,F2,IQM,IQSM,IQBM,SQ,SQ_EFF = ellipsoid_Theta(Q,20,10,2,1,0.15,12.59921049894873) 
    372     plt.subplot(211) 
    373     plt.loglog(Q,0.15*data[3]*data[5]*ellipsoid_volume(20,10)*1e-4) 
    374     plt.loglog(Q,IQSM,'-g',label='IQSM') 
    375     plt.loglog(data[0], data[3]*ellipsoid_volume(20,10)*1e-4,'-k') 
    376     plt.loglog(Q,IQM,'r',label = 'IQM') 
    377     plt.ylim(1e-5,4) 
    378     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    379     plt.show() 
    380     plt.subplot(221) 
    381     plt.loglog(data[0],SQ) 
    382     plt.loglog(data[0],data[5]) 
    383     plt.title('SQ') 
    384     plt.subplot(222) 
    385     plt.semilogx(data[0],SQ/data[5]-1) 
    386     plt.subplot(223) 
    387     plt.title('SQ_EFF') 
    388     plt.loglog(data[0],SQ_EFF) 
    389     plt.loglog(data[0],data[6]) 
    390     plt.subplot(224) 
    391     plt.semilogx(data[0],(SQ_EFF/data[6])-1,label='Sq effective ratio') 
    392     plt.tight_layout() 
    393     plt.show() 
    394  
    395 #SASFIT 
    396 if ELLIPSOID and GAUSSIAN and SASFIT: 
    397     print('\n\n ELLIPSOID COMPARISON WITH SASFIT WITHOUT V') 
    398     #Rp=20,Re=10,eta_core=4,eta_solv=1 
    399     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQM.txt'),dtype=str,delimiter=';').T 
    400     sasqdata=map(float,sasdata[0]) 
    401     sasvaluedata=map(float,sasdata[1]) 
    402     plt.subplot(321) 
    403     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0)[0],'b') 
    404     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata),'g') 
    405     plt.title('IQM') 
    406     plt.subplot(322) 
    407     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0)[0]-1) 
    408     plt.tight_layout() 
    409     plt.show() 
    410     print('========================================================================================') 
    411     print('========================================================================================') 
    412 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15 
    413     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq.txt'),dtype=str,delimiter=';').T 
    414     sasqdata=map(float,sasdata[0]) 
    415     sasvaluedata=map(float,sasdata[1]) 
    416     plt.subplot(321) 
    417     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]) 
    418     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    419     plt.title('SQ poly 10% 0.15') 
    420     plt.subplot(322) 
    421     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]-1) 
    422 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15 
    423     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq2.txt'),dtype=str,delimiter=';').T 
    424     sasqdata=map(float,sasdata[0]) 
    425     sasvaluedata=map(float,sasdata[1]) 
    426     plt.subplot(323) 
    427     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]) 
    428     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    429     plt.title('SQ poly 30% .15') 
    430     plt.subplot(324) 
    431     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]-1) 
    432 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15 
    433     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq3.txt'),dtype=str,delimiter=';').T 
    434     sasqdata=map(float,sasdata[0]) 
    435     sasvaluedata=map(float,sasdata[1]) 
    436     plt.subplot(325) 
    437     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]) 
    438     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    439     plt.title('SQ poly 60% .15') 
    440     plt.subplot(326) 
    441     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]-1) 
    442     plt.tight_layout() 
    443     plt.show()  
    444     print('========================================================================================') 
    445     print('========================================================================================') 
    446 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3 
    447     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq4.txt'),dtype=str,delimiter=';').T 
    448     sasqdata=map(float,sasdata[0]) 
    449     sasvaluedata=map(float,sasdata[1]) 
    450     plt.subplot(321) 
    451     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]) 
    452     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    453     plt.title('SQ poly 10% .3') 
    454     plt.subplot(322) 
    455     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]-1) 
    456 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3 
    457     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq5.txt'),dtype=str,delimiter=';').T 
    458     sasqdata=map(float,sasdata[0]) 
    459     sasvaluedata=map(float,sasdata[1]) 
    460     plt.subplot(323) 
    461     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]) 
    462     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    463     plt.title('SQ poly 30% .3') 
    464     plt.subplot(324) 
    465     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]-1) 
    466 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3 
    467     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq6.txt'),dtype=str,delimiter=';').T 
    468     sasqdata=map(float,sasdata[0]) 
    469     sasvaluedata=map(float,sasdata[1]) 
    470     plt.subplot(325) 
    471     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]) 
    472     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    473     plt.title('SQ poly 60% .3') 
    474     plt.subplot(326) 
    475     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]-1) 
    476     plt.tight_layout() 
    477     plt.show() 
    478     print('========================================================================================') 
    479     print('========================================================================================') 
    480 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6 
    481     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq7.txt'),dtype=str,delimiter=';').T 
    482     sasqdata=map(float,sasdata[0]) 
    483     sasvaluedata=map(float,sasdata[1]) 
    484     plt.subplot(321) 
    485     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]) 
    486     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    487     plt.title('SQ poly 10% .6') 
    488     plt.subplot(322) 
    489     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]-1) 
    490 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6 
    491     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq8.txt'),dtype=str,delimiter=';').T 
    492     sasqdata=map(float,sasdata[0]) 
    493     sasvaluedata=map(float,sasdata[1]) 
    494     plt.subplot(323) 
    495     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]) 
    496     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    497     plt.title('SQ poly 30% .6') 
    498     plt.subplot(324) 
    499     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]-1) 
    500 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6 
    501     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq9.txt'),dtype=str,delimiter=';').T 
    502     sasqdata=map(float,sasdata[0]) 
    503     sasvaluedata=map(float,sasdata[1]) 
    504     plt.subplot(325) 
    505     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]) 
    506     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    507     plt.title('SQ poly 60% .6') 
    508     plt.subplot(326) 
    509     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]-1) 
    510     plt.tight_layout() 
    511     plt.show() 
    512     print('========================================================================================') 
    513     print('========================================================================================') 
    514 #checks beta structre factor 
    515 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15 
    516     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt'),dtype=str,delimiter=';').T 
    517     sasqdata=map(float,sasdata[0]) 
    518     sasvaluedata=map(float,sasdata[1]) 
    519     plt.subplot(321) 
    520     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]) 
    521     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    522     plt.title('SQ_EFF poly 10% 0.15') 
    523     plt.subplot(322) 
    524     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]-1) 
    525 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15 
    526     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt'),dtype=str,delimiter=';').T 
    527     sasqdata=map(float,sasdata[0]) 
    528     sasvaluedata=map(float,sasdata[1]) 
    529     plt.subplot(323) 
    530     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]) 
    531     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    532     plt.title('SQ_EFF poly 30% .15') 
    533     plt.subplot(324) 
    534     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]-1) 
    535 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15 
    536     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt'),dtype=str,delimiter=';').T 
    537     sasqdata=map(float,sasdata[0]) 
    538     sasvaluedata=map(float,sasdata[1]) 
    539     plt.subplot(325) 
    540     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]) 
    541     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    542     plt.title('SQ_EFF poly 60% .15') 
    543     plt.subplot(326) 
    544     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]-1) 
    545     plt.tight_layout() 
    546     plt.show()  
    547     print('========================================================================================') 
    548     print('========================================================================================') 
    549 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3 
    550     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt'),dtype=str,delimiter=';').T 
    551     sasqdata=map(float,sasdata[0]) 
    552     sasvaluedata=map(float,sasdata[1]) 
    553     plt.subplot(321) 
    554     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]) 
    555     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    556     plt.title('SQ_EFF poly 10% .3') 
    557     plt.subplot(322) 
    558     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]-1) 
    559 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3 
    560     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt'),dtype=str,delimiter=';').T 
    561     sasqdata=map(float,sasdata[0]) 
    562     sasvaluedata=map(float,sasdata[1]) 
    563     plt.subplot(323) 
    564     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]) 
    565     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    566     plt.title('SQ_EFF poly 30% .3') 
    567     plt.subplot(324) 
    568     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]-1) 
    569 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3 
    570     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt'),dtype=str,delimiter=';').T 
    571     sasqdata=map(float,sasdata[0]) 
    572     sasvaluedata=map(float,sasdata[1]) 
    573     plt.subplot(325) 
    574     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]) 
    575     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    576     plt.title('SQ_EFF poly 60% .3') 
    577     plt.subplot(326) 
    578     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]-1) 
    579     plt.tight_layout() 
    580     plt.show() 
    581     print('========================================================================================') 
    582     print('========================================================================================') 
    583 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6 
    584     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt'),dtype=str,delimiter=';').T 
    585     sasqdata=map(float,sasdata[0]) 
    586     sasvaluedata=map(float,sasdata[1]) 
    587     plt.subplot(321) 
    588     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]) 
    589     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    590     plt.title('SQ_EFF poly 10% .6') 
    591     plt.subplot(322) 
    592     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]-1) 
    593 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6 
    594     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt'),dtype=str,delimiter=';').T 
    595     sasqdata=map(float,sasdata[0]) 
    596     sasvaluedata=map(float,sasdata[1]) 
    597     plt.subplot(323) 
    598     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]) 
    599     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    600     plt.title('SQ_EFF poly 30% .6') 
    601     plt.subplot(324) 
    602     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]-1) 
    603 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6 
    604     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt'),dtype=str,delimiter=';').T 
    605     sasqdata=map(float,sasdata[0]) 
    606     sasvaluedata=map(float,sasdata[1]) 
    607     plt.subplot(325) 
    608     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]) 
    609     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    610     plt.title('SQ_EFF poly 60% .6') 
    611     plt.subplot(326) 
    612     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]-1) 
    613     plt.tight_layout() 
    614     plt.show() 
    615     print('========================================================================================') 
    616     print('========================================================================================') 
    617 #checks polydispersion for single parameter and varying pd with sasfit 
    618 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
    619     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD.txt'),dtype=str,delimiter=';').T 
    620     sasqdata=map(float,sasdata[0]) 
    621     sasvaluedata=map(float,sasdata[1]) 
    622     plt.subplot(221) 
    623     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.1,radius_equatorial_pd=0)[0]) 
    624     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)) 
    625     plt.title('IQD polar 10%') 
    626     plt.subplot(222) 
    627     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.1,radius_equatorial_pd=0)[0]-1) 
    628 #N=1,s=1,X0=10,distr equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
    629     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD2.txt'),dtype=str,delimiter=';').T 
    630     sasqdata=map(float,sasdata[0]) 
    631     sasvaluedata=map(float,sasdata[1]) 
    632     plt.subplot(223) 
    633     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.1)[0]) 
    634     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)) 
    635     plt.title('IQD equatorial 10%') 
    636     plt.subplot(224) 
    637     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.1)[0]-1) 
    638     plt.tight_layout() 
    639     plt.show() 
    640     print('========================================================================================') 
    641     print('========================================================================================') 
    642     #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
    643     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD3.txt'),dtype=str,delimiter=';').T 
    644     sasqdata=map(float,sasdata[0]) 
    645     sasvaluedata=map(float,sasdata[1]) 
    646     plt.subplot(221) 
    647     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.3,radius_equatorial_pd=0)[0]) 
    648     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)) 
    649     plt.title('IQD polar 30%') 
    650     plt.subplot(222) 
    651     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.3,radius_equatorial_pd=0)[0]-1) 
    652 #N=1,s=3,X0=10,distr equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
    653     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD4.txt'),dtype=str,delimiter=';').T 
    654     sasqdata=map(float,sasdata[0]) 
    655     sasvaluedata=map(float,sasdata[1]) 
    656     plt.subplot(223) 
    657     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.3)[0]) 
    658     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)) 
    659     plt.title('IQD equatorial 30%') 
    660     plt.subplot(224) 
    661     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.3)[0]-1) 
    662     plt.tight_layout() 
    663     plt.show() 
    664     print('========================================================================================') 
    665     print('========================================================================================') 
    666     #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
    667     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD5.txt'),dtype=str,delimiter=';').T 
    668     sasqdata=map(float,sasdata[0]) 
    669     sasvaluedata=map(float,sasdata[1]) 
    670     plt.subplot(221) 
    671     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.6,radius_equatorial_pd=0)[0]) 
    672     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)) 
    673     plt.title('IQD polar 60%') 
    674     plt.subplot(222) 
    675     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.6,radius_equatorial_pd=0)[0]-1) 
    676 #N=1,s=6,X0=10,distr equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
    677     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD6.txt'),dtype=str,delimiter=';').T 
    678     sasqdata=map(float,sasdata[0]) 
    679     sasvaluedata=map(float,sasdata[1]) 
    680     plt.subplot(223) 
    681     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.6)[0]) 
    682     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)) 
    683     plt.title('IQD equatorial 60%') 
    684     plt.subplot(224) 
    685     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.6)[0]-1) 
    686     plt.tight_layout() 
    687     plt.show() 
    688     print('========================================================================================') 
    689     print('========================================================================================') 
    690  
    691 # TEST FOR RICHARD 
     401    model = 'ellipsoid' 
     402    pars = dict( 
     403        radius_polar=20,radius_equatorial=400,sld=4,sld_solvent=1, 
     404        background=0, 
     405        radius_polar_pd=.1, radius_polar_pd_type=pd_type, 
     406        radius_equatorial_pd=.1, radius_equatorial_pd_type=pd_type, 
     407        volfraction=0.15, 
     408        #radius_effective=12.59921049894873, 
     409        ) 
     410    target = sasmodels_theory(q, model, **pars) 
     411    actual = ellipsoid_pe(q, norm='sasview', **pars) 
     412    title = " ".join(("sasmodels", model, pd_type)) 
     413    compare(title, target, actual) 
     414COMPARISON[('sasview','ellipsoid','gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian') 
     415COMPARISON[('sasview','ellipsoid','schulz')] = lambda: compare_sasview_sphere(pd_type='schulz') 
     416 
     417def compare_yun_ellipsoid_mono(): 
     418    pars = { 
     419        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian', 
     420        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian', 
     421        'sld': 2, 'sld_solvent': 1, 
     422        'volfraction': 0.15, 
     423        # Yun uses radius for same volume sphere for effective radius 
     424        # whereas sasview uses the average curvature. 
     425        'radius_effective': 12.59921049894873, 
     426    } 
     427    volume = ellipsoid_volume(pars['radius_polar'], pars['radius_equatorial']) 
     428 
     429    data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T 
     430    Q = data[0] 
     431    F1 = data[1] 
     432    F2 = data[3] 
     433    S = data[5] 
     434    Seff = data[6] 
     435    P = F2 
     436    I = P*S 
     437    Ibeta = P*Seff 
     438    P = I = Ibeta = None 
     439    target = Theory(Q=Q, F1=F1, F2=F2, P=P, S=S, I=I, Seff=Seff, Ibeta=Ibeta) 
     440    actual = ellipsoid_pe(Q, norm='yun', **pars) 
     441    title = " ".join(("yun", "ellipsoid", "no pd")) 
     442    #compare(title, target, actual, fields="P S I Seff Ibeta") 
     443    compare(title, target, actual) 
     444COMPARISON[('yun','ellipsoid','gaussian')] = compare_yun_ellipsoid_mono 
     445COMPARISON[('yun','ellipsoid','schulz')] = compare_yun_ellipsoid_mono 
     446 
     447def compare_sasfit_sphere_gauss(): 
     448    #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3 
     449    pars = { 
     450        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian', 
     451        'sld': 4, 'sld_solvent': 1, 
     452        'volfraction': 0.3, 
     453    } 
     454    volume = 4./3.*pi*pars['radius']**3 
     455    Q, IQ = load_sasfit(data_file('sasfit_sphere_IQD.txt')) 
     456    Q, IQSD = load_sasfit(data_file('sasfit_sphere_IQSD.txt')) 
     457    Q, IQBD = load_sasfit(data_file('sasfit_sphere_IQBD.txt')) 
     458    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_sphere_sq.txt')) 
     459    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_sphere_sqeff.txt')) 
     460    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD) 
     461    actual = sphere_r(Q, norm="sasfit", **pars) 
     462    title = " ".join(("sasfit", "sphere", "pd=10% gaussian")) 
     463    compare(title, target, actual) 
     464    #compare(title, target, actual, fields="P") 
     465COMPARISON[('sasfit','sphere','gaussian')] = compare_sasfit_sphere_gauss 
     466 
     467def compare_sasfit_sphere_schulz(): 
    692468    #radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1 
    693469    #We have scaled the output from sasfit by 1e-4*volume*volfraction 
    694470    #0.10050378152592121 
    695 if (SPHERE == True) & (SCHULZ == True) & (SASVIEW == True) & (SASFIT == True): 
    696     print('                   *****SPHERE MODEL*****') 
    697     sasdata = np.loadtxt(data_file('richard_test.txt'),dtype=str,delimiter=';').T 
    698     sasqdata=map(float,sasdata[0]) 
    699     model=build_model('sphere',sasqdata) 
    700     Iq = model(radius=20,radius_pd=0.1,radius_pd_n=80,sld=4,sld_solvent=1,background=0,radius_pd_type='schulz',radius_pd_nsigma=8) 
    701     model2=build_model('sphere@hardsphere',sasqdata) 
    702     IQS = model2(radius=20,radius_pd=0.1,radius_pd_n=80,sld=4,sld_solvent=1,background=0,radius_pd_type='schulz',radius_pd_nsigma=8,volfraction=0.3) 
    703      
    704     sasvaluedata=map(float,sasdata[1]) 
    705     sasdata2 = np.loadtxt(data_file('richard_test2.txt'),dtype=str,delimiter=';').T 
    706     sasqdata2=map(float,sasdata2[0]) 
    707     sasvaluedata2=map(float,sasdata2[1]) 
    708     sasdata3 = np.loadtxt(data_file('richard_test3.txt'),dtype=str,delimiter=';').T 
    709     sasqdata3=map(float,sasdata3[0]) 
    710     sasvaluedata3=map(float,sasdata3[1]) 
    711     IQD, IQSD, IQBD, SQ, SQ_EFF = sphere_r(np.array(sasqdata),20,4,1,.3,0.1,distribution='schulz') 
    712     plt.grid(True) 
    713     plt.title('IQD(blue) vs. SASVIEW(red) vs. SASFIT(green)') 
    714     plt.loglog(sasqdata,Iq,'b') 
    715     plt.loglog(sasqdata,IQD,'r') 
    716     plt.loglog(sasqdata,1e-4/(4./3.*pi*20**3)*np.array(sasvaluedata),'g') 
    717     plt.show() 
    718     plt.grid(True) 
    719     plt.title('IQSD(blue) vs. SASVIEW(red) vs. SASFIT(green)') 
    720     plt.loglog(sasqdata,IQSD,'b') 
    721     plt.loglog(sasqdata,IQS,'r') 
    722     plt.loglog(sasqdata,1e-4/(4./3.*pi*20**3)*0.3*np.array(sasvaluedata2),'g') 
    723     plt.show() 
    724     plt.grid(True) 
    725     plt.title('IQBD(blue) vs. SASFIT(green)') 
    726     plt.loglog(sasqdata,IQBD,'b') 
    727     plt.loglog(sasqdata,1e-4/(4./3.*pi*20**3)*0.3*np.array(sasvaluedata3),'g') 
    728     plt.show() 
    729  
    730 # TEST FOR RICHARD 
    731 if (ELLIPSOID == True) & (SCHULZ == True) & (SASVIEW == True) & (SASFIT == True):    
     471    pars = { 
     472        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'schulz', 
     473        'sld': 4, 'sld_solvent': 1, 
     474        'volfraction': 0.3, 
     475    } 
     476    volume = 4./3.*pi*pars['radius']**3 
     477 
     478    Q, IQ = load_sasfit(data_file('richard_test.txt')) 
     479    Q, IQSD = load_sasfit(data_file('richard_test2.txt')) 
     480    Q, IQBD = load_sasfit(data_file('richard_test3.txt')) 
     481    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
     482    actual = sphere_r(Q, norm="sasfit", **pars) 
     483    title = " ".join(("sasfit", "sphere", "pd=10% schulz")) 
     484    compare(title, target, actual) 
     485COMPARISON[('sasfit','sphere','schulz')] = compare_sasfit_sphere_schulz 
     486 
     487def compare_sasfit_ellipsoid_schulz(): 
    732488    #polarradius=20, equatorialradius=10, sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1 
    733          #Effective radius =13.1353356684 
    734         #We have scaled the output from sasfit by 1e-4*volume*volfraction 
    735         #0.10050378152592121 
    736     print('                   *****ELLIPSOID MODEL*****') 
    737     sasdata = np.loadtxt(data_file('richard_test4.txt'),dtype=str,delimiter=';').T 
    738     sasqdata=map(float,sasdata[0]) 
    739     model=build_model('ellipsoid',sasqdata) 
    740     Iq = model(radius_polar=20,radius_polar_pd=0.1,radius_polar_pd_n=80, radius_equatorial=10, radius_equatorial_pd=0, sld=4,sld_solvent=1,background=0,radius_polar_pd_type='schulz',radius_polar_pd_nsigma=8) 
    741     model2=build_model('ellipsoid@hardsphere',sasqdata) 
    742     IQS =  model2(radius_polar=20,radius_equatorial=10,sld=4,sld_solvent=1,background=0,radius_polar_pd=.1,\ 
    743                      radius_polar_pd_n=80,radius_equatorial_pd=0,radius_polar_pd_type='schulz',radius_polar_pd_nsigma=8,volfraction=0.3) 
    744     sasvaluedata=map(float,sasdata[1]) 
    745     sasdata2 = np.loadtxt(data_file('richard_test5.txt'),dtype=str,delimiter=';').T 
    746     sasqdata2=map(float,sasdata2[0]) 
    747     sasvaluedata2=map(float,sasdata2[1]) 
    748     sasdata3 = np.loadtxt(data_file('richard_test6.txt'),dtype=str,delimiter=';').T 
    749     sasqdata3=map(float,sasdata3[0]) 
    750     sasvaluedata3=map(float,sasdata3[1]) 
    751     IQD, IQSD, IQBD, SQ, SQ_EFF = ellipsoid_pe(np.array(sasqdata),20,10,4,1,.3,radius_polar_pd=0.1,radius_equatorial_pd=0,distribution='schulz')     
    752     plt.grid(True) 
    753     plt.title('IQD(blue) vs. SASVIEW(red) vs. SASFIT(green)') 
    754     plt.loglog(sasqdata,Iq,'b') 
    755     plt.loglog(sasqdata,IQD,'r') 
    756     plt.loglog(sasqdata,1e-4/(4./3.*pi*20*10**2)*np.array(sasvaluedata),'g') 
    757     plt.show() 
    758     plt.grid(True) 
    759     plt.title('IQSD(blue) vs. SASVIEW(red) vs. SASFIT(green)') 
    760     plt.loglog(sasqdata,IQSD,'b') 
    761     plt.loglog(sasqdata,IQS,'r') 
    762     plt.loglog(sasqdata,1e-4/(4./3.*pi*20*10**2)*0.3*np.array(sasvaluedata2),'g') 
    763     plt.show() 
    764     plt.grid(True) 
    765     plt.title('IQBD(blue) vs. SASFIT(green)') 
    766     plt.loglog(sasqdata,IQBD,'b') 
    767     plt.loglog(sasqdata,1e-4/(4./3.*pi*20*10**2)*0.3*np.array(sasvaluedata3),'g') 
    768     plt.show()     
    769      
    770 # SASVIEW GAUSS 
    771 if (SPHERE == True) & (GAUSSIAN == True) & (SASVIEW == True): 
    772     q = np.logspace(-5, 0, 200) 
    773     IQD, IQSD, IQBD, SQ, SQ_EFF = sphere_r(q,20,4,1,.3) 
    774     model = build_model("sphere", q) 
    775     Sq = model(radius=20,radius_pd=0.1,radius_pd_n=35,sld=4,sld_solvent=1,background=0) 
    776     model2 = build_model("sphere@hardsphere", q) 
    777     S=build_model("hardsphere",q)(radius_effective=20,volfraction=.3,background=0) 
    778     Sq2 = model2(radius=20,radius_pd=0.1,radius_pd_n=35,sld=4,sld_solvent=1,background=0,volfraction=.3) 
    779     print('\n\n SPHERE COMPARISON WITH SASVIEW WITHOUT V') 
    780     plt.subplot(221) 
    781     plt.title('IQD') 
    782     plt.loglog(q, IQD, '-b') 
    783     plt.loglog(q, Sq, '-r') 
    784     plt.subplot(222) 
    785     plt.semilogx(q, Sq/IQD-1, '-g') 
    786     plt.tight_layout() 
    787     plt.show() 
    788     plt.subplot(221) 
    789     plt.title('SQ') 
    790     plt.plot(q, SQ, '-r') 
    791     plt.plot(q,S,'-k') 
    792     plt.subplot(222) 
    793     plt.plot(SQ/S-1) 
    794     plt.tight_layout() 
    795     plt.show() 
    796     plt.subplot(221) 
    797     plt.title('IQSD') 
    798     plt.loglog(q, IQSD, '-b') 
    799     plt.loglog(q, Sq2, '-r',label='IQSD') 
    800     plt.subplot(222) 
    801     plt.semilogx(q, Sq2/IQSD-1, '-g',label='IQSD ratio') 
    802     plt.tight_layout() 
    803     plt.show() 
    804     IQD, IQSD, IQBD, SQ, SQ_EFF = sphere_r(q,20,4,1,0.15) 
    805     plt.subplot(211) 
    806     plt.title('SQ vs SQ_EFF') 
    807     plt.plot(q, SQ) 
    808     plt.plot(q, SQ_EFF) 
    809     plt.tight_layout() 
    810     plt.show()  
    811     plt.subplot(221) 
    812     plt.title('IQSD vs IQBD') 
    813     plt.loglog(q,IQBD) 
    814     plt.loglog(q,IQSD,label='IQBD,IQSD') 
    815     plt.subplot(222) 
    816     plt.plot(q,IQBD) 
    817     plt.plot(q,IQSD,) 
    818     plt.tight_layout() 
    819     plt.show()  
    820  
    821 # SASFIT GAUSS 
    822 if (SPHERE == True) & (GAUSSIAN == True) & (SASFIT == True): 
    823     print('\n\n SPHERE COMPARISON WITH SASFIT WITHOUT V') 
    824 #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3 
    825     sasdata = np.loadtxt(data_file('sasfit_sphere_IQD.txt'),dtype=str,delimiter=';').T 
    826     sasqdata=map(float,sasdata[0]) 
    827     sasvaluedata=map(float,sasdata[1]) 
    828     IQD, IQSD, IQBD, SQ, SQ_EFF = sphere_r(np.array(sasqdata),20,4,1,.3) 
    829     plt.subplot(221) 
    830     plt.title('IQD') 
    831     plt.loglog(sasqdata,IQD ) 
    832     plt.loglog(sasqdata,1e-4/(4./3*pi*20**3)*np.array(sasvaluedata)) 
    833     plt.subplot(222) 
    834     plt.semilogx(sasqdata,1e-4/(4./3*pi*20**3)*np.array(sasvaluedata)/IQD-1) 
    835     plt.tight_layout() 
    836     plt.show() 
    837     sasdata = np.loadtxt(data_file('sasfit_sphere_IQSD.txt'),dtype=str,delimiter=';').T 
    838     sasqdata=map(float,sasdata[0]) 
    839     sasvaluedata=map(float,sasdata[1]) 
    840     plt.subplot(221) 
    841     plt.title('IQSD') 
    842     plt.loglog(sasqdata, IQSD) 
    843     plt.loglog(sasqdata,1e-4/(4./3*pi*20**3)*0.3*np.array(sasvaluedata)) 
    844     plt.subplot(222) 
    845     plt.semilogx(sasqdata,1e-4/(4./3*pi*20**3)*0.3*np.array(sasvaluedata)/IQSD-1) 
    846     plt.tight_layout() 
    847     plt.show() 
    848     sasdata = np.loadtxt(data_file('sasfit_sphere_IQBD.txt'),dtype=str,delimiter=';').T 
    849     sasqdata=map(float,sasdata[0]) 
    850     sasvaluedata=map(float,sasdata[1]) 
    851     plt.subplot(221) 
    852     plt.title('IQBD') 
    853     plt.loglog(sasqdata,IQBD) 
    854     plt.loglog(sasqdata,1e-4/(4./3*pi*20**3)*0.3*np.array(sasvaluedata)) 
    855     plt.subplot(222) 
    856     plt.semilogx(sasqdata,1e-4/(4./3*pi*20**3)*0.3*np.array(sasvaluedata)/IQBD-1) 
    857     plt.tight_layout() 
    858     plt.show() 
    859     sasdata = np.loadtxt(data_file('sasfit_polydisperse_sphere_sq.txt'),dtype=str,delimiter=';').T 
    860     sasqdata=map(float,sasdata[0]) 
    861     sasvaluedata=map(float,sasdata[1]) 
    862     plt.subplot(221) 
    863     plt.title('SQ') 
    864     plt.loglog(sasqdata, SQ) 
    865     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    866     plt.subplot(222) 
    867     plt.semilogx(sasqdata,np.array(sasvaluedata)/SQ-1) 
    868     plt.tight_layout() 
    869     plt.show() 
    870     sasdata = np.loadtxt(data_file('sasfit_polydisperse_sphere_sqeff.txt'),dtype=str,delimiter=';').T 
    871     sasqdata=map(float,sasdata[0]) 
    872     sasvaluedata=map(float,sasdata[1]) 
    873     plt.subplot(221) 
    874     plt.title('SQ_EFF') 
    875     plt.loglog(sasqdata,SQ_EFF) 
    876     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    877     plt.subplot(222) 
    878     plt.semilogx(sasqdata,np.array(sasvaluedata)/SQ_EFF-1) 
    879     plt.tight_layout() 
    880     plt.show() 
     489    #Effective radius =13.1353356684 
     490    #We have scaled the output from sasfit by 1e-4*volume*volfraction 
     491    #0.10050378152592121 
     492    pars = { 
     493        'radius_polar': 20, 'radius_polar_pd': 0.1, 'radius_polar_pd_type': 'schulz', 
     494        'radius_equatorial': 10, 'radius_equatorial_pd': 0., 'radius_equatorial_pd_type': 'schulz', 
     495        'sld': 4, 'sld_solvent': 1, 
     496        'volfraction': 0.3, 'radius_effective': 13.1353356684, 
     497    } 
     498    volume = ellipsoid_volume(pars['radius_polar'], pars['radius_equatorial']) 
     499    Q, IQ = load_sasfit(data_file('richard_test4.txt')) 
     500    Q, IQSD = load_sasfit(data_file('richard_test5.txt')) 
     501    Q, IQBD = load_sasfit(data_file('richard_test6.txt')) 
     502    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
     503    actual = ellipsoid_pe(Q, norm="sasfit", **pars) 
     504    title = " ".join(("sasfit", "ellipsoid", "pd=10% schulz")) 
     505    compare(title, target, actual) 
     506COMPARISON[('sasfit','ellipsoid','schulz')] = compare_sasfit_ellipsoid_schulz 
     507 
     508 
     509def compare_sasfit_ellipsoid_gaussian(): 
     510    pars = { 
     511        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian', 
     512        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian', 
     513        'sld': 4, 'sld_solvent': 1, 
     514        'volfraction': 0, 'radius_effective': None, 
     515    } 
     516 
     517    #Rp=20,Re=10,eta_core=4,eta_solv=1 
     518    Q, PQ0 = load_sasfit(data_file('sasfit_ellipsoid_IQM.txt')) 
     519    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0, radius_effective=None) 
     520    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     521    target = Theory(Q=Q, P=PQ0) 
     522    compare("sasfit ellipsoid no poly", target, actual); plt.show() 
     523 
     524    #N=1,s=2,X0=20,distr 10% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
     525    Q, PQ_Rp10 = load_sasfit(data_file('sasfit_ellipsoid_IQD.txt')) 
     526    pars.update(volfraction=0, radius_polar_pd=0.1, radius_equatorial_pd=0.0, radius_effective=None) 
     527    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     528    target = Theory(Q=Q, P=PQ_Rp10) 
     529    compare("sasfit ellipsoid P(Q) 10% Rp", target, actual); plt.show() 
     530    #N=1,s=1,X0=10,distr 10% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
     531    Q, PQ_Re10 = load_sasfit(data_file('sasfit_ellipsoid_IQD2.txt')) 
     532    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.1, radius_effective=None) 
     533    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     534    target = Theory(Q=Q, P=PQ_Re10) 
     535    compare("sasfit ellipsoid P(Q) 10% Re", target, actual); plt.show() 
     536    #N=1,s=6,X0=20,distr 30% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
     537    Q, PQ_Rp30 = load_sasfit(data_file('sasfit_ellipsoid_IQD3.txt')) 
     538    pars.update(volfraction=0, radius_polar_pd=0.3, radius_equatorial_pd=0.0, radius_effective=None) 
     539    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     540    target = Theory(Q=Q, P=PQ_Rp30) 
     541    compare("sasfit ellipsoid P(Q) 30% Rp", target, actual); plt.show() 
     542    #N=1,s=3,X0=10,distr 30% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
     543    Q, PQ_Re30 = load_sasfit(data_file('sasfit_ellipsoid_IQD4.txt')) 
     544    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.3, radius_effective=None) 
     545    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     546    target = Theory(Q=Q, P=PQ_Re30) 
     547    compare("sasfit ellipsoid P(Q) 30% Re", target, actual); plt.show() 
     548    #N=1,s=12,X0=20,distr 60% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
     549    Q, PQ_Rp60 = load_sasfit(data_file('sasfit_ellipsoid_IQD5.txt')) 
     550    pars.update(volfraction=0, radius_polar_pd=0.6, radius_equatorial_pd=0.0, radius_effective=None) 
     551    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     552    target = Theory(Q=Q, P=PQ_Rp60) 
     553    compare("sasfit ellipsoid P(Q) 60% Rp", target, actual); plt.show() 
     554    #N=1,s=6,X0=10,distr 60% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
     555    Q, PQ_Re60 = load_sasfit(data_file('sasfit_ellipsoid_IQD6.txt')) 
     556    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.6, radius_effective=None) 
     557    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     558    target = Theory(Q=Q, P=PQ_Re60) 
     559    compare("sasfit ellipsoid P(Q) 60% Re", target, actual); plt.show() 
     560 
     561    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15 
     562    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq.txt')) 
     563    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt')) 
     564    pars.update(volfraction=0.15, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254) 
     565    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     566    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     567    compare("sasfit ellipsoid P(Q) 10% Rp 15% Vf", target, actual); plt.show() 
     568    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15 
     569    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq2.txt')) 
     570    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt')) 
     571    pars.update(volfraction=0.15, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149) 
     572    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     573    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     574    compare("sasfit ellipsoid P(Q) 30% Rp 15% Vf", target, actual); plt.show() 
     575    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15 
     576    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq3.txt')) 
     577    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt')) 
     578    pars.update(volfraction=0.15, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917) 
     579    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     580    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     581    compare("sasfit ellipsoid P(Q) 60% Rp 15% Vf", target, actual); plt.show() 
     582 
     583    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3 
     584    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq4.txt')) 
     585    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt')) 
     586    pars.update(volfraction=0.3, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254) 
     587    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     588    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     589    compare("sasfit ellipsoid P(Q) 10% Rp 30% Vf", target, actual); plt.show() 
     590    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3 
     591    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq5.txt')) 
     592    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt')) 
     593    pars.update(volfraction=0.3, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149) 
     594    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     595    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     596    compare("sasfit ellipsoid P(Q) 30% Rp 30% Vf", target, actual); plt.show() 
     597    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3 
     598    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq6.txt')) 
     599    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt')) 
     600    pars.update(volfraction=0.3, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917) 
     601    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     602    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     603    compare("sasfit ellipsoid P(Q) 60% Rp 30% Vf", target, actual); plt.show() 
     604 
     605    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6 
     606    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq7.txt')) 
     607    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt')) 
     608    pars.update(volfraction=0.6, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254) 
     609    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     610    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     611    compare("sasfit ellipsoid P(Q) 10% Rp 60% Vf", target, actual); plt.show() 
     612    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6 
     613    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq8.txt')) 
     614    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt')) 
     615    pars.update(volfraction=0.6, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149) 
     616    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     617    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     618    compare("sasfit ellipsoid P(Q) 30% Rp 60% Vf", target, actual); plt.show() 
     619    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6 
     620    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq9.txt')) 
     621    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt')) 
     622    pars.update(volfraction=0.6, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917) 
     623    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     624    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     625    compare("sasfit ellipsoid P(Q) 60% Rp 60% Vf", target, actual); plt.show() 
     626COMPARISON[('sasfit','ellipsoid','gaussian')] = compare_sasfit_ellipsoid_schulz 
     627 
     628def main(): 
     629    key = tuple(sys.argv[1:]) 
     630    if key not in COMPARISON: 
     631        print("usage: sasfit_compare.py [sasview|sasfit|yun] [sphere|ellipsoid] [gaussian|schulz]") 
     632        return 
     633    comparison = COMPARISON[key] 
     634    comparison() 
     635 
     636if __name__ == "__main__": 
     637    main() 
Note: See TracChangeset for help on using the changeset viewer.