Changes in / [eb63414:a57b31d] in sasmodels


Ignore:
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • example/fit.py

    rf4b36fa r1182da5  
    2424            % section) 
    2525data = radial_data if section != "tangential" else tan_data 
    26 theta = 89.9 if section != "tangential" else 0 
    27 phi = 90 
     26phi = 0 if section != "tangential" else 90 
    2827kernel = load_model(name, dtype="single") 
    2928cutoff = 1e-3 
    … …  
    3130if name == "ellipsoid": 
    3231    model = Model(kernel, 
    33         scale=0.08, background=35, 
    34         radius_polar=15, radius_equatorial=800, 
     32        scale=0.08, 
     33        r_polar=15, r_equatorial=800, 
    3534        sld=.291, sld_solvent=7.105, 
    36         theta=theta, phi=phi, 
    37         theta_pd=0, theta_pd_n=0, theta_pd_nsigma=3, 
     35        background=0, 
     36        theta=90, phi=phi, 
     37        theta_pd=15, theta_pd_n=40, theta_pd_nsigma=3, 
     38        r_polar_pd=0.222296, r_polar_pd_n=1, r_polar_pd_nsigma=0, 
     39        r_equatorial_pd=.000128, r_equatorial_pd_n=1, r_equatorial_pd_nsigma=0, 
    3840        phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3, 
    39         radius_polar_pd=0.222296, radius_polar_pd_n=1, radius_polar_pd_nsigma=0, 
    40         radius_equatorial_pd=.000128, radius_equatorial_pd_n=1, radius_equatorial_pd_nsigma=0, 
    4141        ) 
    4242 
     43 
    4344    # SET THE FITTING PARAMETERS 
    44     model.radius_polar.range(15, 1000) 
    45     model.radius_equatorial.range(15, 1000) 
    46     #model.theta.range(0, 90) 
    47     #model.theta_pd.range(0,10) 
    48     model.phi_pd.range(0,20) 
    49     model.phi.range(0, 180) 
     45    model.r_polar.range(15, 1000) 
     46    model.r_equatorial.range(15, 1000) 
     47    model.theta_pd.range(0, 360) 
    5048    model.background.range(0,1000) 
    5149    model.scale.range(0, 10) 
    5250 
    5351 
     52 
    5453elif name == "lamellar": 
    5554    model = Model(kernel, 
    56         scale=0.08, background=0.003, 
     55        scale=0.08, 
    5756        thickness=19.2946, 
    5857        sld=5.38,sld_sol=7.105, 
     58        background=0.003, 
    5959        thickness_pd= 0.37765, thickness_pd_n=10, thickness_pd_nsigma=3, 
    6060        ) 
     61 
    6162 
    6263    # SET THE FITTING PARAMETERS 
    … …  
    7677        radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3, 
    7778        length_pd=0.493, length_pd_n=10, length_pd_nsigma=3, 
    78         phi_pd=0, phi_pd_n=5 phi_pd_nsigma=3,) 
     79        phi_pd=0, phi_pd_n=5, phi_pd_nsigma=3,) 
    7980        """ 
    8081    pars = dict( 
    8182        scale=.01, background=35, 
    8283        sld=.291, sld_solvent=5.77, 
    83         radius=250, length=178, 
     84        radius=250, length=178,  
     85        theta=90, phi=phi, 
    8486        radius_pd=0.1, radius_pd_n=5, radius_pd_nsigma=3, 
    8587        length_pd=0.1,length_pd_n=5, length_pd_nsigma=3, 
    86         theta=theta, phi=phi, 
    87         theta_pd=0, theta_pd_n=0, theta_pd_nsigma=3, 
    88         phi_pd=10, phi_pd_n=20, phi_pd_nsigma=3) 
     88        theta_pd=10, theta_pd_n=50, theta_pd_nsigma=3, 
     89        phi_pd=0, phi_pd_n=10, phi_pd_nsigma=3) 
    8990    model = Model(kernel, **pars) 
    9091 
    … …  
    9293    model.radius.range(1, 500) 
    9394    model.length.range(1, 5000) 
    94     #model.theta.range(0, 90) 
    95     model.phi.range(0, 180) 
    96     model.phi_pd.range(0, 30) 
     95    model.theta.range(-90,100) 
     96    model.theta_pd.range(0, 30) 
     97    model.theta_pd_n = model.theta_pd + 5 
    9798    model.radius_pd.range(0, 1) 
    98     model.length_pd.range(0, 1) 
     99    model.length_pd.range(0, 2) 
    99100    model.scale.range(0, 10) 
    100101    model.background.range(0, 100) 
    … …  
    103104elif name == "core_shell_cylinder": 
    104105    model = Model(kernel, 
    105         scale= .031, background=0, 
    106         radius=19.5, thickness=30, length=22, 
    107         sld_core=7.105, sld_shell=.291, sld_solvent=7.105, 
     106        scale= .031, radius=19.5, thickness=30, length=22, 
     107        sld_core=7.105, sld_shell=.291, sdl_solvent=7.105, 
     108        background=0, theta=0, phi=phi, 
     109 
    108110        radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3, 
    109111        length_pd=0.26, length_pd_n=10, length_pd_nsigma=3, 
    110112        thickness_pd=1, thickness_pd_n=1, thickness_pd_nsigma=1, 
    111         theta=theta, phi=phi, 
    112         theta_pd=1, theta_pd_n=1, theta_pd_nsigma=3, 
    113         phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3, 
     113        theta_pd=1, theta_pd_n=10, theta_pd_nsigma=3, 
     114        phi_pd=0.1, phi_pd_n=1, phi_pd_nsigma=1, 
    114115        ) 
    115116 
    116117    # SET THE FITTING PARAMETERS 
    117     model.radius.range(115, 1000) 
    118     model.length.range(0, 2500) 
     118    #model.radius.range(115, 1000) 
     119    #model.length.range(0, 2500) 
    119120    #model.thickness.range(18, 38) 
    120121    #model.thickness_pd.range(0, 1) 
    121122    #model.phi.range(0, 90) 
    122     model.phi_pd.range(0,20) 
    123123    #model.radius_pd.range(0, 1) 
    124124    #model.length_pd.range(0, 1) 
    … …  
    131131elif name == "capped_cylinder": 
    132132    model = Model(kernel, 
    133         scale=.08, background=35, 
    134         radius=20, cap_radius=40, length=400, 
     133        scale=.08, radius=20, cap_radius=40, length=400, 
    135134        sld=1, sld_solvent=6.3, 
     135        background=0, theta=0, phi=phi, 
    136136        radius_pd=.1, radius_pd_n=5, radius_pd_nsigma=3, 
    137137        cap_radius_pd=.1, cap_radius_pd_n=5, cap_radius_pd_nsigma=3, 
    138138        length_pd=.1, length_pd_n=1, length_pd_nsigma=0, 
    139         theta=theta, phi=phi, 
    140         theta_pd=0, theta_pd_n=1, theta_pd_nsigma=0, 
    141         phi_pd=10, phi_pd_n=20, phi_pd_nsigma=0, 
     139        theta_pd=.1, theta_pd_n=1, theta_pd_nsigma=0, 
     140        phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
    142141        ) 
    143142 
    144     model.radius.range(115, 1000) 
    145     model.length.range(0, 2500) 
    146     #model.thickness.range(18, 38) 
    147     #model.thickness_pd.range(0, 1) 
    148     #model.phi.range(0, 90) 
    149     model.phi_pd.range(0,20) 
    150     #model.radius_pd.range(0, 1) 
    151     #model.length_pd.range(0, 1) 
    152     #model.theta_pd.range(0, 360) 
    153     #model.background.range(0,5) 
    154143    model.scale.range(0, 1) 
    155144 
    … …  
    157146elif name == "triaxial_ellipsoid": 
    158147    model = Model(kernel, 
    159         scale=0.08, background=35, 
    160         radius_equat_minor=15, radius_equat_major=20, radius_polar=500, 
     148        scale=0.08, req_minor=15, req_major=20, rpolar=500, 
    161149        sld=7.105, solvent_sld=.291, 
    162         radius_equat_minor_pd=.1, radius_equat_minor_pd_n=1, radius_equat_minor_pd_nsigma=0, 
    163         radius_equat_major_pd=.1, radius_equat_major_pd_n=1, radius_equat_major_pd_nsigma=0, 
    164         radius_polar_pd=.1, radius_polar_pd_n=1, radius_polar_pd_nsigma=0, 
    165         theta=theta, phi=phi, psi=0, 
     150        background=5, theta=0, phi=phi, psi=0, 
    166151        theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3, 
    167152        phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
    168153        psi_pd=30, psi_pd_n=1, psi_pd_nsigma=0, 
     154        req_minor_pd=.1, req_minor_pd_n=1, req_minor_pd_nsigma=0, 
     155        req_major_pd=.1, req_major_pd_n=1, req_major_pd_nsigma=0, 
     156        rpolar_pd=.1, rpolar_pd_n=1, rpolar_pd_nsigma=0, 
    169157        ) 
    170158 
    171159    # SET THE FITTING PARAMETERS 
    172     model.radius_equat_minor.range(15, 1000) 
    173     model.radius_equat_major.range(15, 1000) 
    174     #model.radius_polar.range(15, 1000) 
     160    model.req_minor.range(15, 1000) 
     161    model.req_major.range(15, 1000) 
     162    #model.rpolar.range(15, 1000) 
    175163    #model.background.range(0,1000) 
    176164    #model.theta_pd.range(0, 360) 
    … …  
    185173M = Experiment(data=data, model=model) 
    186174if section == "both": 
    187    tan_model = Model(model.sasmodel, **model.parameters()) 
     175   tan_model = Model(model.kernel, **model.parameters()) 
    188176   tan_model.phi = model.phi - 90 
    189177   tan_model.cutoff = cutoff 

  • sasmodels/models/rpa.c

    r20c856a rfa4a994  
    3939  double S31,S32,S33,S34,S41,S42,S43,S44; 
    4040  double Lad,Lbd,Lcd,Nav,Intg; 
    41    
    42   // Set values for non existent parameters (eg. no A or B in case 0 and 1 etc) 
     41 
    4342  //icase was shifted to N-1 from the original code 
    4443  if (icase <= 1){ 
    … …  
    5857    Kab = Kac = Kad = -0.0004; 
    5958  } 
    60   
    61   // Set volume fraction of component D based on constraint that sum of vol frac =1 
    62   Phi[3]=1.0-Phi[0]-Phi[1]-Phi[2]; 
    63  
    64   //set up values for cross terms in case of block copolymers (1,3,4,6,7,8,9) 
     59 
    6560  Nab=sqrt(N[0]*N[1]); 
    6661  Nac=sqrt(N[0]*N[2]); 
    … …  
    8479  Phicd=sqrt(Phi[2]*Phi[3]); 
    8580 
    86   // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk  
    8781  Xa=q*q*b[0]*b[0]*N[0]/6.0; 
    8882  Xb=q*q*b[1]*b[1]*N[1]/6.0; 
    … …  
    9084  Xd=q*q*b[3]*b[3]*N[3]/6.0; 
    9185 
    92   //calculate all partial structure factors Pij and normalize n^2 
    93   Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); // free A chain form factor 
    94   S0aa=N[0]*Phi[0]*v[0]*Paa; // Phi * Vp * P(Q)= I(Q0)/delRho^2 
    95   Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); //AB diblock (anchored Paa * anchored Pbb) partial form factor 
     86  Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); 
     87  S0aa=N[0]*Phi[0]*v[0]*Paa; 
     88  Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); 
    9689  S0ab=(Phiab*vab*Nab)*Pab; 
    97   Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); //ABC triblock AC partial form factor 
     90  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); 
    9891  S0ac=(Phiac*vac*Nac)*Pac; 
    99   Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); //ABCD four block 
     92  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); 
    10093  S0ad=(Phiad*vad*Nad)*Pad; 
    10194 
    10295  S0ba=S0ab; 
    103   Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); // free B chain 
     96  Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); 
    10497  S0bb=N[1]*Phi[1]*v[1]*Pbb; 
    105   Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); // BC diblock 
     98  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); 
    10699  S0bc=(Phibc*vbc*Nbc)*Pbc; 
    107   Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); // BCD triblock 
     100  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); 
    108101  S0bd=(Phibd*vbd*Nbd)*Pbd; 
    109102 
    110103  S0ca=S0ac; 
    111104  S0cb=S0bc; 
    112   Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); // Free C chain 
     105  Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); 
    113106  S0cc=N[2]*Phi[2]*v[2]*Pcc; 
    114   Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); // CD diblock 
     107  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); 
    115108  S0cd=(Phicd*vcd*Ncd)*Pcd; 
    116109 
    … …  
    118111  S0db=S0bd; 
    119112  S0dc=S0cd; 
    120   Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); // free D chain 
     113  Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); 
    121114  S0dd=N[3]*Phi[3]*v[3]*Pdd; 
    122  
    123   // Reset all unused partial structure factors to 0 (depends on case) 
    124115  //icase was shifted to N-1 from the original code 
    125116  switch(icase){ 
    … …  
    202193  S0dc=S0cd; 
    203194 
    204   // self chi parameter is 0 ... of course 
    205195  Kaa=0.0; 
    206196  Kbb=0.0; 
    … …  
    253243  ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd); 
    254244 
    255   // D is considered the matrix or background component so enters here 
    256245  m=1.0/(S0dd-ZZ); 
    257246 
    … …  
    308297  S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23; 
    309298 
    310   //calculate contrast where L[i] is the scattering length of i and D is the matrix 
    311   //Note that should multiply by Nav to get units of SLD which will become 
    312   // Nav*2 in the next line (SLD^2) but then normalization in that line would 
    313   //need to divide by Nav leaving only Nav or sqrt(Nav) before squaring.  
    314299  Nav=6.022045e+23; 
    315300  Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav); 
    … …  
    318303 
    319304  Intg=Lad*Lad*S11+Lbd*Lbd*S22+Lcd*Lcd*S33+2.0*Lad*Lbd*S12+2.0*Lbd*Lcd*S23+2.0*Lad*Lcd*S13; 
    320  
    321   //rescale for units of Lij^2 (fm^2 to cm^2) 
    322   Intg *= 1.0e-26;     
    323305 
    324306  return Intg; 

  • sasmodels/models/rpa.py

    r20c856a rfa4a994  
    11r""" 
    2 Definition 
    3 ---------- 
    4  
    52Calculates the macroscopic scattering intensity for a multi-component 
    63homogeneous mixture of polymers using the Random Phase Approximation. 
    … …  
    2724Case 9: A-B-C-D tetra-block copolymer 
    2825 
    29 .. note:: 
    30     These case numbers are different from those in the NIST SANS package! 
     26**NB: these case numbers are different from those in the NIST SANS package!** 
    3127 
    32 USAGE NOTES: 
     28Only one case can be used at any one time. 
    3329 
    34 * Only one case can be used at any one time. 
    35 * The RPA (mean field) formalism only applies only when the multicomponent 
    36   polymer mixture is in the homogeneous mixed-phase region. 
    37 * **Component D is assumed to be the "background" component (ie, all contrasts 
    38   are calculated with respect to component D).** So the scattering contrast 
    39   for a C/D blend = [SLD(component C) - SLD(component D)]\ :sup:`2`. 
    40 * Depending on which case is being used, the number of fitting parameters can  
    41   vary.  Note that in general the degrees of polymerization, the volume 
    42   fractions, the molar volumes, and the neutron scattering lengths for each 
    43   component are obtained from other methods and held fixed while the segment 
    44   lengths (b\ :sub:`a`, b\ :sub:`b`, etc) and $\chi$ parameters (K\ :sub:`ab`, 
    45   K\ :sub:`ac`, etc). The *scale* parameter should be held equal to unity. 
     30The RPA (mean field) formalism only applies only when the multicomponent 
     31polymer mixture is in the homogeneous mixed-phase region. 
     32 
     33**Component D is assumed to be the "background" component (ie, all contrasts 
     34are calculated with respect to component D).** So the scattering contrast 
     35for a C/D blend = [SLD(component C) - SLD(component D)]\ :sup:`2`. 
     36 
     37Depending on which case is being used, the number of fitting parameters - the 
     38segment lengths (ba, bb, etc) and $\chi$ parameters (Kab, Kac, etc) - vary. 
     39The *scale* parameter should be held equal to unity. 
     40 
     41The input parameters are the degrees of polymerization, the volume fractions, 
     42the specific volumes, and the neutron scattering length densities for each 
     43component. 
    4644 
    4745 
    … …  
    4947---------- 
    5048 
    51 .. [#] A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136 
     49A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136 
    5250""" 
    5351 
    … …  
    5553 
    5654name = "rpa" 
    57 title = "Random Phase Approximation" 
     55title = "Random Phase Approximation - unfinished work in progress" 
    5856description = """ 
    5957This formalism applies to multicomponent polymer mixtures in the 
    … …  
    9290    ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    9391    ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 
    94     ["v[4]", "mL/mol", 100.0, [0, inf], "", "molar volume"], 
     92    ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
    9593    ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    9694    ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], 
    … …  
    109107 
    110108control = "case_num" 
    111 HIDE_ALL = set("Phi4".split()) 
    112 HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()).union(HIDE_ALL) 
     109HIDE_NONE = set() 
     110HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()) 
    113111HIDE_AB = set("N2 Phi2 v2 L2 b2 K23 K24".split()).union(HIDE_A) 
    114112def hidden(case_num): 
    … …  
    122120        return HIDE_A 
    123121    else: 
    124         return HIDE_ALL 
     122        return HIDE_NONE 
    125123 
Note: See TracChangeset for help on using the changeset viewer.