Changes in / [a57b31d:eb63414] in sasmodels


Ignore:
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • example/fit.py

    r1182da5 rf4b36fa  
    2424            % section) 
    2525data = radial_data if section != "tangential" else tan_data 
    26 phi = 0 if section != "tangential" else 90 
     26theta = 89.9 if section != "tangential" else 0 
     27phi = 90 
    2728kernel = load_model(name, dtype="single") 
    2829cutoff = 1e-3 
    … …  
    3031if name == "ellipsoid": 
    3132    model = Model(kernel, 
    32         scale=0.08, 
    33         r_polar=15, r_equatorial=800, 
     33        scale=0.08, background=35, 
     34        radius_polar=15, radius_equatorial=800, 
    3435        sld=.291, sld_solvent=7.105, 
    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, 
     36        theta=theta, phi=phi, 
     37        theta_pd=0, theta_pd_n=0, theta_pd_nsigma=3, 
    4038        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  
    4443    # SET THE FITTING PARAMETERS 
    45     model.r_polar.range(15, 1000) 
    46     model.r_equatorial.range(15, 1000) 
    47     model.theta_pd.range(0, 360) 
     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) 
    4850    model.background.range(0,1000) 
    4951    model.scale.range(0, 10) 
    5052 
    5153 
    52  
    5354elif name == "lamellar": 
    5455    model = Model(kernel, 
    55         scale=0.08, 
     56        scale=0.08, background=0.003, 
    5657        thickness=19.2946, 
    5758        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  
    6261 
    6362    # SET THE FITTING PARAMETERS 
    … …  
    7776        radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3, 
    7877        length_pd=0.493, length_pd_n=10, length_pd_nsigma=3, 
    79         phi_pd=0, phi_pd_n=5, phi_pd_nsigma=3,) 
     78        phi_pd=0, phi_pd_n=5 phi_pd_nsigma=3,) 
    8079        """ 
    8180    pars = dict( 
    8281        scale=.01, background=35, 
    8382        sld=.291, sld_solvent=5.77, 
    84         radius=250, length=178,  
    85         theta=90, phi=phi, 
     83        radius=250, length=178, 
    8684        radius_pd=0.1, radius_pd_n=5, radius_pd_nsigma=3, 
    8785        length_pd=0.1,length_pd_n=5, length_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) 
     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) 
    9089    model = Model(kernel, **pars) 
    9190 
    … …  
    9392    model.radius.range(1, 500) 
    9493    model.length.range(1, 5000) 
    95     model.theta.range(-90,100) 
    96     model.theta_pd.range(0, 30) 
    97     model.theta_pd_n = model.theta_pd + 5 
     94    #model.theta.range(0, 90) 
     95    model.phi.range(0, 180) 
     96    model.phi_pd.range(0, 30) 
    9897    model.radius_pd.range(0, 1) 
    99     model.length_pd.range(0, 2) 
     98    model.length_pd.range(0, 1) 
    10099    model.scale.range(0, 10) 
    101100    model.background.range(0, 100) 
    … …  
    104103elif name == "core_shell_cylinder": 
    105104    model = Model(kernel, 
    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  
     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, 
    110108        radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3, 
    111109        length_pd=0.26, length_pd_n=10, length_pd_nsigma=3, 
    112110        thickness_pd=1, thickness_pd_n=1, thickness_pd_nsigma=1, 
    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, 
     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, 
    115114        ) 
    116115 
    117116    # SET THE FITTING PARAMETERS 
    118     #model.radius.range(115, 1000) 
    119     #model.length.range(0, 2500) 
     117    model.radius.range(115, 1000) 
     118    model.length.range(0, 2500) 
    120119    #model.thickness.range(18, 38) 
    121120    #model.thickness_pd.range(0, 1) 
    122121    #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, radius=20, cap_radius=40, length=400, 
     133        scale=.08, background=35, 
     134        radius=20, cap_radius=40, length=400, 
    134135        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_pd=.1, theta_pd_n=1, theta_pd_nsigma=0, 
    140         phi_pd=.1, phi_pd_n=1, phi_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, 
    141142        ) 
    142143 
     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) 
    143154    model.scale.range(0, 1) 
    144155 
    … …  
    146157elif name == "triaxial_ellipsoid": 
    147158    model = Model(kernel, 
    148         scale=0.08, req_minor=15, req_major=20, rpolar=500, 
     159        scale=0.08, background=35, 
     160        radius_equat_minor=15, radius_equat_major=20, radius_polar=500, 
    149161        sld=7.105, solvent_sld=.291, 
    150         background=5, theta=0, phi=phi, psi=0, 
     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, 
    151166        theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3, 
    152167        phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
    153168        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, 
    157169        ) 
    158170 
    159171    # SET THE FITTING PARAMETERS 
    160     model.req_minor.range(15, 1000) 
    161     model.req_major.range(15, 1000) 
    162     #model.rpolar.range(15, 1000) 
     172    model.radius_equat_minor.range(15, 1000) 
     173    model.radius_equat_major.range(15, 1000) 
     174    #model.radius_polar.range(15, 1000) 
    163175    #model.background.range(0,1000) 
    164176    #model.theta_pd.range(0, 360) 
    … …  
    173185M = Experiment(data=data, model=model) 
    174186if section == "both": 
    175    tan_model = Model(model.kernel, **model.parameters()) 
     187   tan_model = Model(model.sasmodel, **model.parameters()) 
    176188   tan_model.phi = model.phi - 90 
    177189   tan_model.cutoff = cutoff 

  • sasmodels/models/rpa.c

    rfa4a994 r20c856a  
    3939  double S31,S32,S33,S34,S41,S42,S43,S44; 
    4040  double Lad,Lbd,Lcd,Nav,Intg; 
    41  
     41   
     42  // Set values for non existent parameters (eg. no A or B in case 0 and 1 etc) 
    4243  //icase was shifted to N-1 from the original code 
    4344  if (icase <= 1){ 
    … …  
    5758    Kab = Kac = Kad = -0.0004; 
    5859  } 
    59  
     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) 
    6065  Nab=sqrt(N[0]*N[1]); 
    6166  Nac=sqrt(N[0]*N[2]); 
    … …  
    7984  Phicd=sqrt(Phi[2]*Phi[3]); 
    8085 
     86  // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk  
    8187  Xa=q*q*b[0]*b[0]*N[0]/6.0; 
    8288  Xb=q*q*b[1]*b[1]*N[1]/6.0; 
    … …  
    8490  Xd=q*q*b[3]*b[3]*N[3]/6.0; 
    8591 
    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); 
     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 
    8996  S0ab=(Phiab*vab*Nab)*Pab; 
    90   Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); 
     97  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); //ABC triblock AC partial form factor 
    9198  S0ac=(Phiac*vac*Nac)*Pac; 
    92   Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); 
     99  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); //ABCD four block 
    93100  S0ad=(Phiad*vad*Nad)*Pad; 
    94101 
    95102  S0ba=S0ab; 
    96   Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); 
     103  Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); // free B chain 
    97104  S0bb=N[1]*Phi[1]*v[1]*Pbb; 
    98   Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); 
     105  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); // BC diblock 
    99106  S0bc=(Phibc*vbc*Nbc)*Pbc; 
    100   Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); 
     107  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); // BCD triblock 
    101108  S0bd=(Phibd*vbd*Nbd)*Pbd; 
    102109 
    103110  S0ca=S0ac; 
    104111  S0cb=S0bc; 
    105   Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); 
     112  Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); // Free C chain 
    106113  S0cc=N[2]*Phi[2]*v[2]*Pcc; 
    107   Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); 
     114  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); // CD diblock 
    108115  S0cd=(Phicd*vcd*Ncd)*Pcd; 
    109116 
    … …  
    111118  S0db=S0bd; 
    112119  S0dc=S0cd; 
    113   Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); 
     120  Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); // free D chain 
    114121  S0dd=N[3]*Phi[3]*v[3]*Pdd; 
     122 
     123  // Reset all unused partial structure factors to 0 (depends on case) 
    115124  //icase was shifted to N-1 from the original code 
    116125  switch(icase){ 
    … …  
    193202  S0dc=S0cd; 
    194203 
     204  // self chi parameter is 0 ... of course 
    195205  Kaa=0.0; 
    196206  Kbb=0.0; 
    … …  
    243253  ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd); 
    244254 
     255  // D is considered the matrix or background component so enters here 
    245256  m=1.0/(S0dd-ZZ); 
    246257 
    … …  
    297308  S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23; 
    298309 
     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.  
    299314  Nav=6.022045e+23; 
    300315  Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav); 
    … …  
    303318 
    304319  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;     
    305323 
    306324  return Intg; 

  • sasmodels/models/rpa.py

    rfa4a994 r20c856a  
    11r""" 
     2Definition 
     3---------- 
     4 
    25Calculates the macroscopic scattering intensity for a multi-component 
    36homogeneous mixture of polymers using the Random Phase Approximation. 
    … …  
    2427Case 9: A-B-C-D tetra-block copolymer 
    2528 
    26 **NB: these case numbers are different from those in the NIST SANS package!** 
     29.. note:: 
     30    These case numbers are different from those in the NIST SANS package! 
    2731 
    28 Only one case can be used at any one time. 
     32USAGE NOTES: 
    2933 
    30 The RPA (mean field) formalism only applies only when the multicomponent 
    31 polymer mixture is in the homogeneous mixed-phase region. 
    32  
    33 **Component D is assumed to be the "background" component (ie, all contrasts 
    34 are calculated with respect to component D).** So the scattering contrast 
    35 for a C/D blend = [SLD(component C) - SLD(component D)]\ :sup:`2`. 
    36  
    37 Depending on which case is being used, the number of fitting parameters - the 
    38 segment lengths (ba, bb, etc) and $\chi$ parameters (Kab, Kac, etc) - vary. 
    39 The *scale* parameter should be held equal to unity. 
    40  
    41 The input parameters are the degrees of polymerization, the volume fractions, 
    42 the specific volumes, and the neutron scattering length densities for each 
    43 component. 
     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. 
    4446 
    4547 
    … …  
    4749---------- 
    4850 
    49 A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136 
     51.. [#] A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136 
    5052""" 
    5153 
    … …  
    5355 
    5456name = "rpa" 
    55 title = "Random Phase Approximation - unfinished work in progress" 
     57title = "Random Phase Approximation" 
    5658description = """ 
    5759This formalism applies to multicomponent polymer mixtures in the 
    … …  
    9092    ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    9193    ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 
    92     ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
     94    ["v[4]", "mL/mol", 100.0, [0, inf], "", "molar volume"], 
    9395    ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    9496    ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], 
    … …  
    107109 
    108110control = "case_num" 
    109 HIDE_NONE = set() 
    110 HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()) 
     111HIDE_ALL = set("Phi4".split()) 
     112HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()).union(HIDE_ALL) 
    111113HIDE_AB = set("N2 Phi2 v2 L2 b2 K23 K24".split()).union(HIDE_A) 
    112114def hidden(case_num): 
    … …  
    120122        return HIDE_A 
    121123    else: 
    122         return HIDE_NONE 
     124        return HIDE_ALL 
    123125 
Note: See TracChangeset for help on using the changeset viewer.