Changes in / [c68c5e9:b216880] in sasmodels
- Files:
-
- 28 edited
Legend:
- Unmodified
- Added
- Removed
-
example/fit.py
rf4b36fa r1182da5 24 24 % section) 25 25 data = radial_data if section != "tangential" else tan_data 26 theta = 89.9 if section != "tangential" else 0 27 phi = 90 26 phi = 0 if section != "tangential" else 90 28 27 kernel = load_model(name, dtype="single") 29 28 cutoff = 1e-3 … … 31 30 if name == "ellipsoid": 32 31 model = Model(kernel, 33 scale=0.08, background=35,34 r adius_polar=15, radius_equatorial=800,32 scale=0.08, 33 r_polar=15, r_equatorial=800, 35 34 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, 38 40 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,41 41 ) 42 42 43 43 44 # 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) 50 48 model.background.range(0,1000) 51 49 model.scale.range(0, 10) 52 50 53 51 52 54 53 elif name == "lamellar": 55 54 model = Model(kernel, 56 scale=0.08, background=0.003,55 scale=0.08, 57 56 thickness=19.2946, 58 57 sld=5.38,sld_sol=7.105, 58 background=0.003, 59 59 thickness_pd= 0.37765, thickness_pd_n=10, thickness_pd_nsigma=3, 60 60 ) 61 61 62 62 63 # SET THE FITTING PARAMETERS … … 76 77 radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3, 77 78 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,) 79 80 """ 80 81 pars = dict( 81 82 scale=.01, background=35, 82 83 sld=.291, sld_solvent=5.77, 83 radius=250, length=178, 84 radius=250, length=178, 85 theta=90, phi=phi, 84 86 radius_pd=0.1, radius_pd_n=5, radius_pd_nsigma=3, 85 87 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) 89 90 model = Model(kernel, **pars) 90 91 … … 92 93 model.radius.range(1, 500) 93 94 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 97 98 model.radius_pd.range(0, 1) 98 model.length_pd.range(0, 1)99 model.length_pd.range(0, 2) 99 100 model.scale.range(0, 10) 100 101 model.background.range(0, 100) … … 103 104 elif name == "core_shell_cylinder": 104 105 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 108 110 radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3, 109 111 length_pd=0.26, length_pd_n=10, length_pd_nsigma=3, 110 112 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, 114 115 ) 115 116 116 117 # 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) 119 120 #model.thickness.range(18, 38) 120 121 #model.thickness_pd.range(0, 1) 121 122 #model.phi.range(0, 90) 122 model.phi_pd.range(0,20)123 123 #model.radius_pd.range(0, 1) 124 124 #model.length_pd.range(0, 1) … … 131 131 elif name == "capped_cylinder": 132 132 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, 135 134 sld=1, sld_solvent=6.3, 135 background=0, theta=0, phi=phi, 136 136 radius_pd=.1, radius_pd_n=5, radius_pd_nsigma=3, 137 137 cap_radius_pd=.1, cap_radius_pd_n=5, cap_radius_pd_nsigma=3, 138 138 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, 142 141 ) 143 142 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)154 143 model.scale.range(0, 1) 155 144 … … 157 146 elif name == "triaxial_ellipsoid": 158 147 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, 161 149 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, 166 151 theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3, 167 152 phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 168 153 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, 169 157 ) 170 158 171 159 # SET THE FITTING PARAMETERS 172 model.r adius_equat_minor.range(15, 1000)173 model.r adius_equat_major.range(15, 1000)174 #model.r adius_polar.range(15, 1000)160 model.req_minor.range(15, 1000) 161 model.req_major.range(15, 1000) 162 #model.rpolar.range(15, 1000) 175 163 #model.background.range(0,1000) 176 164 #model.theta_pd.range(0, 360) … … 185 173 M = Experiment(data=data, model=model) 186 174 if section == "both": 187 tan_model = Model(model. sasmodel, **model.parameters())175 tan_model = Model(model.kernel, **model.parameters()) 188 176 tan_model.phi = model.phi - 90 189 177 tan_model.cutoff = cutoff -
sasmodels/compare.py
rd504bcd rfe25eda 340 340 if pars['radius'] < pars['thick_string']: 341 341 pars['radius'], pars['thick_string'] = pars['thick_string'], pars['radius'] 342 pars['num_pearls'] = math.ceil(pars['num_pearls']) 342 343 pass 343 344 … … 352 353 for c in '1234': 353 354 pars['Phi'+c] /= total 355 356 elif name == 'stacked_disks': 357 pars['n_stacking'] = math.ceil(pars['n_stacking']) 354 358 355 359 def parlist(model_info, pars, is2d): -
sasmodels/models/core_multi_shell.c
r925ad6e r925ad6e 8 8 } 9 9 10 static double 11 form_volume(double core_radius, double fp_n, double thickness[]) 10 double 11 form_volume(double core_radius, double n, double thickness[]); 12 double 13 form_volume(double core_radius, double n, double thickness[]) 12 14 { 13 15 double r = core_radius; 14 int n = (int)(fp_n+0.5);15 16 for (int i=0; i < n; i++) { 16 17 r += thickness[i]; … … 19 20 } 20 21 21 staticdouble22 double 22 23 Iq(double q, double core_sld, double core_radius, 23 double solvent_sld, double fp_n, double sld[], double thickness[]) 24 double solvent_sld, double num_shells, double sld[], double thickness[]); 25 double 26 Iq(double q, double core_sld, double core_radius, 27 double solvent_sld, double num_shells, double sld[], double thickness[]) 24 28 { 25 const int n = (int) (fp_n+0.5);29 const int n = (int)ceil(num_shells); 26 30 double f, r, last_sld; 27 31 r = core_radius; -
sasmodels/models/core_multi_shell.py
r5a0b3d7 r925ad6e 107 107 Returns the SLD profile *r* (Ang), and *rho* (1e-6/Ang^2). 108 108 """ 109 n = int(n+0.5)110 109 z = [] 111 110 rho = [] … … 134 133 def ER(radius, n, thickness): 135 134 """Effective radius""" 136 n = int(n[0] +0.5) # n is a control parameter and is notpolydisperse135 n = int(n[0]) # n cannot be polydisperse 137 136 return np.sum(thickness[:n], axis=0) + radius 138 137 -
sasmodels/models/flexible_cylinder.c
r592343f r592343f 1 1 static double 2 form_volume( double length, double kuhn_length, doubleradius)2 form_volume(length, kuhn_length, radius) 3 3 { 4 4 return 1.0; -
sasmodels/models/lamellar_hg_stack_caille.c
r1c6286d ra807206 3 3 */ 4 4 5 static double 6 Iq(double qval, 7 double length_tail, 8 double length_head, 9 double fp_Nlayers, 10 double dd, 11 double Cp, 12 double tail_sld, 13 double head_sld, 14 double solvent_sld) 5 double Iq(double qval, 6 double length_tail, 7 double length_head, 8 double Nlayers, 9 double dd, 10 double Cp, 11 double tail_sld, 12 double head_sld, 13 double solvent_sld); 14 15 double Iq(double qval, 16 double length_tail, 17 double length_head, 18 double Nlayers, 19 double dd, 20 double Cp, 21 double tail_sld, 22 double head_sld, 23 double solvent_sld) 15 24 { 16 int Nlayers = (int)(fp_Nlayers+0.5); //cast to an integer for the loop25 double NN; //local variables of coefficient wave 17 26 double inten,Pq,Sq,alpha,temp,t2; 18 27 //double dQ, dQDefault, t1, t3; 28 int ii,NNint; 19 29 // from wikipedia 0.577215664901532860606512090082402431042159335 20 30 const double Euler = 0.577215664901533; // Euler's constant, increased sig figs for new models Feb 2015 … … 22 32 //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015 23 33 24 Pq = (head_sld-solvent_sld)*(sin(qval*(length_head+length_tail))-sin(qval*length_tail)) 25 + (tail_sld-solvent_sld)*sin(qval*length_tail); 34 NN = trunc(Nlayers); //be sure that NN is an integer 35 36 Pq = (head_sld-solvent_sld)*(sin(qval*(length_head+length_tail))-sin(qval*length_tail)) + 37 (tail_sld-solvent_sld)*sin(qval*length_tail); 26 38 Pq *= Pq; 27 39 Pq *= 4.0/(qval*qval); 28 40 41 NNint = (int)NN; //cast to an integer for the loop 42 ii=0; 29 43 Sq = 0.0; 30 for(int ii=1; ii < Nlayers; ii++) { 44 for(ii=1;ii<=(NNint-1);ii+=1) { 45 46 //fii = (double)ii; //do I really need to do this? - unused variable, removed 18Feb2015 47 31 48 temp = 0.0; 32 49 alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler); … … 35 52 //t3 = dQ*dQ*dd*dd*ii*ii; 36 53 37 temp = 1.0- (double)ii/(double)Nlayers;54 temp = 1.0-ii/NN; 38 55 //temp *= cos(dd*qval*ii/(1.0+t1)); 39 56 temp *= cos(dd*qval*ii); 40 57 //if (temp < 0) printf("q=%g: ii=%d, cos(dd*q*ii)=cos(%g) < 0\n",qval,ii,dd*qval*ii); 41 58 //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 42 temp *= exp(-t2/2.0 );59 temp *= exp(-t2/2.0 ); 43 60 //temp /= sqrt(1.0+t1); 44 61 … … 54 71 55 72 inten *= 1.0e-04; // 1/A to 1/cm 56 return inten;73 return(inten); 57 74 } 58 75 -
sasmodels/models/lamellar_hg_stack_caille.py
ra57b31d r7c57861 98 98 ["length_head", "Ang", 2, [0, inf], "volume", 99 99 "head thickness"], 100 ["Nlayers", "", 30, [ 1, inf], "",100 ["Nlayers", "", 30, [0, inf], "", 101 101 "Number of layers"], 102 102 ["d_spacing", "Ang", 40., [0.0, inf], "volume", -
sasmodels/models/lamellar_stack_caille.c
r1c6286d r0bef47b 3 3 */ 4 4 5 static double 6 Iq(double qval, 7 double del, 8 double fp_Nlayers, 9 double dd, 10 double Cp, 11 double sld, 12 double solvent_sld) 5 double Iq(double qval, 6 double del, 7 double Nlayers, 8 double dd, 9 double Cp, 10 double sld, 11 double solvent_sld); 12 13 double Iq(double qval, 14 double del, 15 double Nlayers, 16 double dd, 17 double Cp, 18 double sld, 19 double solvent_sld) 13 20 { 14 int Nlayers = (int)(fp_Nlayers+0.5); //cast to an integer for the loop 15 double contr; //local variables of coefficient wave 21 double contr,NN; //local variables of coefficient wave 16 22 double inten,Pq,Sq,alpha,temp,t2; 17 23 //double dQ, dQDefault, t1, t3; 24 int ii,NNint; 18 25 // from wikipedia 0.577215664901532860606512090082402431042159335 19 26 const double Euler = 0.577215664901533; // Euler's constant, increased sig figs for new models Feb 2015 … … 21 28 //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015 22 29 30 NN = trunc(Nlayers); //be sure that NN is an integer 31 23 32 contr = sld - solvent_sld; 24 33 25 34 Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)); 26 35 36 NNint = (int)NN; //cast to an integer for the loop 37 ii=0; 27 38 Sq = 0.0; 28 for (int ii=1; ii < Nlayers; ii++) { 39 // the vital "=" in ii<= added March 2015 40 for(ii=1;ii<=(NNint-1);ii+=1) { 41 42 //fii = (double)ii; //do I really need to do this? - unused variable, removed 18Feb2015 43 29 44 temp = 0.0; 30 45 alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler); … … 33 48 //t3 = dQ*dQ*dd*dd*ii*ii; 34 49 35 temp = 1.0 - (double)ii / (double)Nlayers;50 temp = 1.0-ii/NN; 36 51 //temp *= cos(dd*qval*ii/(1.0+t1)); 37 52 temp *= cos(dd*qval*ii); 38 53 //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 39 temp *= exp(-t2/2.0 );54 temp *= exp(-t2/2.0 ); 40 55 //temp /= sqrt(1.0+t1); 41 56 -
sasmodels/models/lamellar_stack_caille.py
ra57b31d r7c57861 88 88 parameters = [ 89 89 ["thickness", "Ang", 30.0, [0, inf], "volume", "sheet thickness"], 90 ["Nlayers", "", 20, [ 1, inf], "", "Number of layers"],90 ["Nlayers", "", 20, [0, inf], "", "Number of layers"], 91 91 ["d_spacing", "Ang", 400., [0.0,inf], "volume", "lamellar d-spacing of Caille S(Q)"], 92 92 ["Caille_parameter", "1/Ang^2", 0.1, [0.0,0.8], "", "Caille parameter"], -
sasmodels/models/lamellar_stack_paracrystal.c
r5467cd8 r4962519 2 2 3 3 */ 4 double paraCryst_sn(double ww, double qval, double davg, int Nlayers, double an); 5 double paraCryst_an(double ww, double qval, double davg, int Nlayers); 4 double Iq(double qval, 5 double th, 6 double Nlayers, 7 double davg, 8 double pd, 9 double sld, 10 double solvent_sld); 11 double paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an); 12 double paraCryst_an(double ww, double qval, double davg, long Nlayers); 6 13 7 static double 8 Iq(double qval, 9 double th, 10 double fp_Nlayers, 11 double davg, 12 double pd, 13 double sld, 14 double solvent_sld) 14 double Iq(double qval, 15 double th, 16 double Nlayers, 17 double davg, 18 double pd, 19 double sld, 20 double solvent_sld) 15 21 { 22 23 double inten,contr,xn; 24 double xi,ww,Pbil,Znq,Snq,an; 25 long n1,n2; 26 27 contr = sld - solvent_sld; 16 28 //get the fractional part of Nlayers, to determine the "mixing" of N's 17 int n1 = (int)(fp_Nlayers); //truncate towards zero18 int n2 = n1 + 1;19 const double xn = (double)n2 - fp_Nlayers; //fractional contribution of n120 29 21 const double ww = exp(-0.5*square(qval*pd*davg)); 30 n1 = (long)trunc(Nlayers); //rounds towards zero 31 n2 = n1 + 1; 32 xn = (double)n2 - Nlayers; //fractional contribution of n1 33 34 ww = exp(-qval*qval*pd*pd*davg*davg/2.0); 22 35 23 36 //calculate the n1 contribution 24 double Znq,Snq,an;25 37 an = paraCryst_an(ww,qval,davg,n1); 26 38 Snq = paraCryst_sn(ww,qval,davg,n1,an); 27 39 28 40 Znq = xn*Snq; 29 41 … … 40 52 // Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg)) 41 53 42 const doublexi = th/2.0; //use 1/2 the bilayer thickness43 const double Pbil = square(sas_sinx_x(qval*xi));54 xi = th/2.0; //use 1/2 the bilayer thickness 55 Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi)); 44 56 45 const double contr = sld - solvent_sld;46 const double inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval);57 inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval); 58 inten *= 1.0e-04; 47 59 //printf("q=%.7e wwm1=%g ww=%.5e an=% 12.5e Snq=% 12.5e Znq=% 12.5e Pbil=% 12.5e\n",qval,wwm1,ww,an,Snq,Znq,Pbil); 48 return 1.0e-4*inten;60 return(inten); 49 61 } 50 62 51 63 // functions for the lamellar paracrystal model 52 64 double 53 paraCryst_sn(double ww, double qval, double davg, intNlayers, double an) {65 paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an) { 54 66 55 67 double Snq; … … 57 69 Snq = an/( (double)Nlayers*square(1.0+ww*ww-2.0*ww*cos(qval*davg)) ); 58 70 59 return Snq;71 return(Snq); 60 72 } 61 73 62 74 double 63 paraCryst_an(double ww, double qval, double davg, int Nlayers) { 75 paraCryst_an(double ww, double qval, double davg, long Nlayers) { 76 64 77 double an; 65 78 … … 69 82 an += 2.0*pow(ww,(Nlayers+1))*cos((double)(Nlayers+1)*qval*davg); 70 83 71 return an;84 return(an); 72 85 } 73 86 -
sasmodels/models/lamellar_stack_paracrystal.py
ra0168e8 r7c57861 113 113 parameters = [["thickness", "Ang", 33.0, [0, inf], "volume", 114 114 "sheet thickness"], 115 ["Nlayers", "", 20, [ 1, inf], "",115 ["Nlayers", "", 20, [0, inf], "", 116 116 "Number of layers"], 117 117 ["d_spacing", "Ang", 250., [0.0, inf], "", -
sasmodels/models/linear_pearls.c
r925ad6e r925ad6e 4 4 double radius, 5 5 double edge_sep, 6 double fp_num_pearls,6 double num_pearls, 7 7 double pearl_sld, 8 8 double solvent_sld); … … 11 11 double radius, 12 12 double edge_sep, 13 intnum_pearls,13 double num_pearls, 14 14 double pearl_sld, 15 15 double solvent_sld); 16 16 17 17 18 double form_volume(double radius, double fp_num_pearls)18 double form_volume(double radius, double num_pearls) 19 19 { 20 int num_pearls = (int)(fp_num_pearls + 0.5);21 20 // Pearl volume 22 21 double pearl_vol = M_4PI_3 * cube(radius); … … 28 27 double radius, 29 28 double edge_sep, 30 intnum_pearls,29 double num_pearls, 31 30 double pearl_sld, 32 31 double solvent_sld) 33 32 { 33 double n_contrib; 34 34 //relative sld 35 35 double contrast_pearl = pearl_sld - solvent_sld; … … 46 46 double psi = sas_3j1x_x(q * radius); 47 47 48 // N pearls interaction terms 49 double structure_factor = (double)num_pearls; 50 for(int num=1; num<num_pearls; num++) { 51 structure_factor += 2.0*(num_pearls-num)*sas_sinx_x(q*separation*num); 48 // N pearls contribution 49 int n_max = num_pearls - 1; 50 n_contrib = num_pearls; 51 for(int num=1; num<=n_max; num++) { 52 n_contrib += (2.0*(num_pearls-num)*sas_sinx_x(q*separation*num)); 52 53 } 53 54 // form factor for num_pearls 54 double form_factor = 1.0e-4 * structure_factor* square(m_s*psi) / tot_vol;55 double form_factor = 1.0e-4 * n_contrib * square(m_s*psi) / tot_vol; 55 56 56 57 return form_factor; … … 60 61 double radius, 61 62 double edge_sep, 62 double fp_num_pearls,63 double num_pearls, 63 64 double pearl_sld, 64 65 double solvent_sld) 65 66 { 66 67 67 int num_pearls = (int)(fp_num_pearls + 0.5);68 68 double result = linear_pearls_kernel(q, 69 69 radius, -
sasmodels/models/linear_pearls.py
r925ad6e r925ad6e 16 16 .. math:: 17 17 18 P(Q) = \frac{ \text{scale}}{V}\left[ m_{p}^219 \left(N+2\sum_{n-1}^{N-1}(N-n)\frac{ \sin(qnl)}{qnl}\right)20 \left( 3\frac{ \sin(qR)-qR\cos(qR)}{(qr)^3}\right)^2\right]18 P(Q) = \frac{scale}{V}\left[ m_{p}^2 19 \left(N+2\sum_{n-1}^{N-1}(N-n)\frac{sin(qnl)}{qnl}\right) 20 \left( 3\frac{sin(qR)-qRcos(qR)}{(qr)^3}\right)^2\right] 21 21 22 22 where the mass $m_p$ is $(SLD_{pearl}-SLD_{solvent})*(volume\ of\ N\ pearls)$. … … 56 56 ["radius", "Ang", 80.0, [0, inf], "", "Radius of the pearls"], 57 57 ["edge_sep", "Ang", 350.0, [0, inf], "", "Length of the string segment - surface to surface"], 58 ["num_pearls", "", 3.0, [ 1, inf], "", "Number of the pearls"],58 ["num_pearls", "", 3.0, [0, inf], "", "Number of the pearls"], 59 59 ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", "SLD of the pearl spheres"], 60 60 ["sld_solvent", "1e-6/Ang^2", 6.3, [-inf, inf], "sld", "SLD of the solvent"], -
sasmodels/models/multilayer_vesicle.c
r925ad6e r925ad6e 7 7 double sld_solvent, 8 8 double sld, 9 intn_pairs)9 double n_pairs) 10 10 { 11 11 //calculate with a loop, two shells at a time … … 47 47 double sld_solvent, 48 48 double sld, 49 double fp_n_pairs)49 double n_pairs) 50 50 { 51 int n_pairs = (int)(fp_n_pairs + 0.5);52 51 return multilayer_vesicle_kernel(q, 53 52 volfraction, -
sasmodels/models/multilayer_vesicle.py
r925ad6e r925ad6e 19 19 20 20 .. math:: 21 P(q) = \text{scale} \cdot \frac{V_f}{V_t} F^2(q) + \text{background}22 21 23 for 22 P(q) = \frac{\text{scale.volfraction}}{V_t} F^2(q) + \text{background} 23 24 where 24 25 25 26 .. math:: 26 F(q) = (\rho_\text{shell}-\rho_\text{solv}) \sum_{i=1}^{n_\text{pairs}}27 \left[28 3V(R_i)\frac{\sin(qR_i)-qR_i\cos(qR_i)}{(qR_i)^3} \\29 - 3V(R_i+t_s)\frac{\sin(q(R_i+t_s))-q(R_i+t_s)\cos(q(R_i+t_s))}{(q(R_i+t_s))^3}30 \right]31 27 32 and 28 F(q) = (\rho_{shell}-\rho_{solv}) \sum_{i=1}^{n\_pairs} \left[ 29 3V(R_i)\frac{\sin(qR_i)-qR_i\cos(qR_i)}{(qR_i)^3} \\ 30 - 3V(R_i+t_s)\frac{\sin(q(R_i+t_s))-q(R_i+t_s)\cos(q(R_i+t_s))}{(q(R_i+t_s))^3} 31 \right] 33 32 34 .. math::35 R_i = r_c + (i-1)(t_s + t_w)36 33 37 where $V_f$ is the volume fraction of particles, $V_t$ is the volume of the 38 whole particle, $V(r)$ is the volume of a sphere of radius $r$, $r_c$ is the 39 radius of the core, $\rho_\text{shell}$ is the scattering length density of a 40 shell, $\rho_\text{solv}$ is the scattering length density of the solvent. 34 where $R_i = r_c + (i-1)(t_s + t_w)$ 35 36 where $V_t$ is the volume of the whole particle, $V(R)$ is the volume of a sphere 37 of radius $R$, $r_c$ is the radius of the core, $\rho_{shell}$ is the scattering length 38 density of a shell, $\rho_{solv}$ is the scattering length density of the solvent. 41 39 42 The outer most radius, $r_o = R_n + t_s$, is used for both the volume fraction43 normalization and for the effective radius for *S(Q)* when $P(Q) * S(Q)$44 is applied.45 40 46 41 The 2D scattering intensity is the same as 1D, regardless of the orientation … … 50 45 51 46 q = \sqrt{q_x^2 + q_y^2} 47 48 49 The outer most radius 50 51 $radius + n\_pairs * thick\_shell + (n\_pairs- 1) * thick\_solvent$ 52 53 is used for both the volume fraction normalization and for the 54 effective radius for *S(Q)* when $P(Q) * S(Q)$ is applied. 52 55 53 56 For information about polarised and magnetic scattering, see … … 67 70 **Author:** NIST IGOR/DANSE **on:** pre 2010 68 71 69 **Last Modified by:** Piotr Rozyczko 72 **Last Modified by:** Piotr Rozyczko**on:** Feb 24, 2016 70 73 71 74 **Last Reviewed by:** Paul Butler **on:** March 20, 2016 … … 106 109 source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 107 110 108 # TODO: the following line does nothing109 111 polydispersity = ["radius", "n_pairs"] 110 112 -
sasmodels/models/onion.c
r925ad6e r925ad6e 30 30 31 31 static double 32 form_volume(double radius_core, double n _shells, double thickness[])32 form_volume(double radius_core, double n, double thickness[]) 33 33 { 34 int n = (int)(n_shells+0.5);34 int i; 35 35 double r = radius_core; 36 for (i nt i=0; i < n; i++) {36 for (i=0; i < n; i++) { 37 37 r += thickness[i]; 38 38 } -
sasmodels/models/onion.py
r925ad6e r925ad6e 323 323 Returns shape profile with x=radius, y=SLD. 324 324 """ 325 n_shells = int(n_shells+0.5) 325 326 326 total_radius = 1.25*(sum(thickness[:n_shells]) + radius_core + 1) 327 327 dz = total_radius/400 # 400 points for a smooth plot … … 366 366 return np.asarray(z), np.asarray(rho) 367 367 368 def ER(radius_core, n _shells, thickness):368 def ER(radius_core, n, thickness): 369 369 """Effective radius""" 370 n = int(n_shells[0]+0.5) 371 return np.sum(thickness[:n], axis=0) + radius_core 370 return np.sum(thickness[:int(n[0])], axis=0) + radius_core 372 371 373 372 demo = { -
sasmodels/models/pearl_necklace.c
r3f853beb r4b541ac 1 1 double form_volume(double radius, double edge_sep, 2 double thick_string, double fp_num_pearls);2 double thick_string, double num_pearls); 3 3 double Iq(double q, double radius, double edge_sep, 4 double thick_string, double fp_num_pearls, double sld,4 double thick_string, double num_pearls, double sld, 5 5 double string_sld, double solvent_sld); 6 6 … … 9 9 // From Igor library 10 10 static double 11 pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string,12 intnum_pearls, double sld_pearl, double sld_string, double sld_solv)11 _pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string, 12 double num_pearls, double sld_pearl, double sld_string, double sld_solv) 13 13 { 14 14 // number of string segments 15 const int num_strings = num_pearls - 1; 15 num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 16 const double num_strings = num_pearls - 1.0; 16 17 17 18 //each masses: contrast * volume … … 68 69 69 70 double form_volume(double radius, double edge_sep, 70 double thick_string, double fp_num_pearls)71 double thick_string, double num_pearls) 71 72 { 72 const int num_pearls = (int)(fp_num_pearls + 0.5); //Force integer number of pearls 73 const int num_strings = num_pearls - 1; 73 num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 74 75 const double num_strings = num_pearls - 1.0; 74 76 const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string; 75 77 const double pearl_vol = M_4PI_3 * radius * radius * radius; … … 80 82 81 83 double Iq(double q, double radius, double edge_sep, 82 double thick_string, double fp_num_pearls, double sld,84 double thick_string, double num_pearls, double sld, 83 85 double string_sld, double solvent_sld) 84 86 { 85 const int num_pearls = (int)(fp_num_pearls + 0.5); //Force integer number of pearls 86 const double form = pearl_necklace_kernel(q, radius, edge_sep, 87 const double form = _pearl_necklace_kernel(q, radius, edge_sep, 87 88 thick_string, num_pearls, sld, string_sld, solvent_sld); 88 89 -
sasmodels/models/pearl_necklace.py
r1bd1ea2 r4b541ac 82 82 ["thick_string", "Ang", 2.5, [0, inf], "volume", 83 83 "Thickness of the chain linkage"], 84 ["num_pearls", "none", 3, [ 1, inf], "volume",84 ["num_pearls", "none", 3, [0, inf], "volume", 85 85 "Number of pearls in the necklace (must be integer)"], 86 86 ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", … … 100 100 Redundant with form_volume. 101 101 """ 102 num_pearls = int(num_pearls + 0.5)103 102 number_of_strings = num_pearls - 1.0 104 103 string_vol = edge_sep * pi * pow((thick_string / 2.0), 2.0) … … 112 111 Calculation for effective radius. 113 112 """ 114 num_pearls = int(num_pearls + 0.5)115 113 tot_vol = volume(radius, edge_sep, thick_string, num_pearls) 116 114 rad_out = (tot_vol/(4.0/3.0*pi)) ** (1./3.) -
sasmodels/models/raspberry.py
r8e68ea0 r8e68ea0 10 10 Schematic of the raspberry model 11 11 12 In order to calculate the form factor of the entire complex, the 13 self-correlation of the large droplet, the self-correlation of the particles, 14 the correlation terms between different particles and the cross terms between 15 largedroplet and small particles all need to be calculated.12 In order to calculate the form factor of the entire complex, the self- 13 correlation of the large droplet, the self-correlation of the particles, the 14 correlation terms between different particles and the cross terms between large 15 droplet and small particles all need to be calculated. 16 16 17 Consider two infinitely thin shells of radii $R_1$ and $R_2$ separated by18 distance $r$. The general structure of the equation is then the form factor 19 of the two shells multiplied by the phase factor that accounts for the 20 separation of theircenters.17 Consider two infinitely thin shells of radii R1 and R2 separated by distance r. 18 The general structure of the equation is then the form factor of the two shells 19 multiplied by the phase factor that accounts for the separation of their 20 centers. 21 21 22 22 .. math:: -
sasmodels/models/rpa.c
r20c856a r6351bfa 1 double Iq(double q, double fp_case_num,1 double Iq(double q, double case_num, 2 2 double N[], double Phi[], double v[], double L[], double b[], 3 3 double Kab, double Kac, double Kad, … … 5 5 ); 6 6 7 double Iq(double q, double fp_case_num,7 double Iq(double q, double case_num, 8 8 double N[], // DEGREE OF POLYMERIZATION 9 9 double Phi[], // VOL FRACTION … … 15 15 ) 16 16 { 17 int icase = (int) (fp_case_num+0.5);17 int icase = (int)case_num; 18 18 19 19 double Nab,Nac,Nad,Nbc,Nbd,Ncd; … … 39 39 double S31,S32,S33,S34,S41,S42,S43,S44; 40 40 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 43 42 //icase was shifted to N-1 from the original code 44 43 if (icase <= 1){ … … 58 57 Kab = Kac = Kad = -0.0004; 59 58 } 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 65 60 Nab=sqrt(N[0]*N[1]); 66 61 Nac=sqrt(N[0]*N[2]); … … 84 79 Phicd=sqrt(Phi[2]*Phi[3]); 85 80 86 // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk87 81 Xa=q*q*b[0]*b[0]*N[0]/6.0; 88 82 Xb=q*q*b[1]*b[1]*N[1]/6.0; … … 90 84 Xd=q*q*b[3]*b[3]*N[3]/6.0; 91 85 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); 96 89 S0ab=(Phiab*vab*Nab)*Pab; 97 Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); //ABC triblock AC partial form factor90 Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); 98 91 S0ac=(Phiac*vac*Nac)*Pac; 99 Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); //ABCD four block92 Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); 100 93 S0ad=(Phiad*vad*Nad)*Pad; 101 94 102 95 S0ba=S0ab; 103 Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); // free B chain96 Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); 104 97 S0bb=N[1]*Phi[1]*v[1]*Pbb; 105 Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); // BC diblock98 Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); 106 99 S0bc=(Phibc*vbc*Nbc)*Pbc; 107 Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); // BCD triblock100 Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); 108 101 S0bd=(Phibd*vbd*Nbd)*Pbd; 109 102 110 103 S0ca=S0ac; 111 104 S0cb=S0bc; 112 Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); // Free C chain105 Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); 113 106 S0cc=N[2]*Phi[2]*v[2]*Pcc; 114 Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); // CD diblock107 Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); 115 108 S0cd=(Phicd*vcd*Ncd)*Pcd; 116 109 … … 118 111 S0db=S0bd; 119 112 S0dc=S0cd; 120 Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); // free D chain113 Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); 121 114 S0dd=N[3]*Phi[3]*v[3]*Pdd; 122 123 // Reset all unused partial structure factors to 0 (depends on case)124 115 //icase was shifted to N-1 from the original code 125 116 switch(icase){ … … 202 193 S0dc=S0cd; 203 194 204 // self chi parameter is 0 ... of course205 195 Kaa=0.0; 206 196 Kbb=0.0; … … 253 243 ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd); 254 244 255 // D is considered the matrix or background component so enters here256 245 m=1.0/(S0dd-ZZ); 257 246 … … 308 297 S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23; 309 298 310 //calculate contrast where L[i] is the scattering length of i and D is the matrix311 //Note that should multiply by Nav to get units of SLD which will become312 // Nav*2 in the next line (SLD^2) but then normalization in that line would313 //need to divide by Nav leaving only Nav or sqrt(Nav) before squaring.314 299 Nav=6.022045e+23; 315 300 Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav); … … 318 303 319 304 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;323 305 324 306 return Intg; -
sasmodels/models/rpa.py
r20c856a r40a87fa 1 1 r""" 2 Definition3 ----------4 5 2 Calculates the macroscopic scattering intensity for a multi-component 6 3 homogeneous mixture of polymers using the Random Phase Approximation. … … 27 24 Case 9: A-B-C-D tetra-block copolymer 28 25 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!** 31 27 32 USAGE NOTES: 28 Only one case can be used at any one time. 33 29 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. 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. 46 44 47 45 … … 49 47 ---------- 50 48 51 .. [#]A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 413649 A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136 52 50 """ 53 51 … … 55 53 56 54 name = "rpa" 57 title = "Random Phase Approximation "55 title = "Random Phase Approximation - unfinished work in progress" 58 56 description = """ 59 57 This formalism applies to multicomponent polymer mixtures in the … … 92 90 ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 93 91 ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 94 ["v[4]", "mL/mol", 100.0, [0, inf], "", " molarvolume"],92 ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 95 93 ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 96 94 ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], … … 109 107 110 108 control = "case_num" 111 HIDE_ ALL = set("Phi4".split())112 HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()) .union(HIDE_ALL)109 HIDE_NONE = set() 110 HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()) 113 111 HIDE_AB = set("N2 Phi2 v2 L2 b2 K23 K24".split()).union(HIDE_A) 114 112 def hidden(case_num): … … 116 114 Return a list of parameters to hide depending on the multiplicity parameter. 117 115 """ 118 case_num = int(case_num+0.5)119 116 if case_num < 2: 120 117 return HIDE_AB … … 122 119 return HIDE_A 123 120 else: 124 return HIDE_ ALL121 return HIDE_NONE 125 122 -
sasmodels/models/spherical_sld.c
r925ad6e r925ad6e 1 1 static double form_volume( 2 double fp_n_shells,2 int n_shells, 3 3 double thickness[], 4 4 double interface[]) 5 5 { 6 int n_shells= (int)(fp_n_shells + 0.5);7 6 double r = 0.0; 8 7 for (int i=0; i < n_shells; i++) { … … 21 20 return pow(z, nu); 22 21 } else if (shape==2) { 23 return 1.0 - pow(1. 0- z, nu);22 return 1.0 - pow(1. - z, nu); 24 23 } else if (shape==3) { 25 24 return expm1(-nu*z)/expm1(-nu); … … 45 44 static double Iq( 46 45 double q, 47 double fp_n_shells,46 int n_shells, 48 47 double sld_solvent, 49 48 double sld[], … … 52 51 double shape[], 53 52 double nu[], 54 double fp_n_steps)53 int n_steps) 55 54 { 56 55 // iteration for # of shells + core + solvent 57 int n_shells = (int)(fp_n_shells + 0.5);58 int n_steps = (int)(fp_n_steps + 0.5);59 56 double f=0.0; 60 57 double r=0.0; -
sasmodels/models/spherical_sld.py
r925ad6e r925ad6e 233 233 """ 234 234 235 n_shells = int(n_shells + 0.5)236 n_steps = int(n_steps + 0.5)237 235 z = [] 238 236 rho = [] … … 242 240 rho.append(sld[0]) 243 241 244 for i in range(0, n_shells):242 for i in range(0, int(n_shells)): 245 243 z_next += thickness[i] 246 244 z.append(z_next) … … 263 261 def ER(n_shells, thickness, interface): 264 262 """Effective radius""" 265 n_shells = int(n_shells + 0.5)263 n_shells = int(n_shells) 266 264 total = (np.sum(thickness[:n_shells], axis=1) 267 265 + np.sum(interface[:n_shells], axis=1)) -
sasmodels/models/stacked_disks.c
r19f996b r6c3e266 1 static double stacked_disks_kernel( 2 double q, 3 double halfheight, 4 double thick_layer, 5 double radius, 6 int n_stacking, 7 double sigma_dnn, 8 double core_sld, 9 double layer_sld, 10 double solvent_sld, 11 double sin_alpha, 12 double cos_alpha, 13 double d) 1 double form_volume(double thick_core, 2 double thick_layer, 3 double radius, 4 double n_stacking); 5 6 double Iq(double q, 7 double thick_core, 8 double thick_layer, 9 double radius, 10 double n_stacking, 11 double sigma_dnn, 12 double core_sld, 13 double layer_sld, 14 double solvent_sld); 15 16 double Iqxy(double qx, double qy, 17 double thick_core, 18 double thick_layer, 19 double radius, 20 double n_stacking, 21 double sigma_dnn, 22 double core_sld, 23 double layer_sld, 24 double solvent_sld, 25 double theta, 26 double phi); 27 28 static 29 double _kernel(double q, 30 double radius, 31 double core_sld, 32 double layer_sld, 33 double solvent_sld, 34 double halfheight, 35 double thick_layer, 36 double sin_alpha, 37 double cos_alpha, 38 double sigma_dnn, 39 double d, 40 double n_stacking) 14 41 15 42 { … … 61 88 62 89 63 static double stacked_disks_1d(64 65 double thick_core,66 double thick_layer,67 double radius,68 intn_stacking,69 double sigma_dnn,70 double core_sld,71 double layer_sld,72 double solvent_sld)90 static 91 double stacked_disks_kernel(double q, 92 double thick_core, 93 double thick_layer, 94 double radius, 95 double n_stacking, 96 double sigma_dnn, 97 double core_sld, 98 double layer_sld, 99 double solvent_sld) 73 100 { 74 101 /* StackedDiscsX : calculates the form factor of a stacked "tactoid" of core shell disks … … 84 111 double sin_alpha, cos_alpha; // slots to hold sincos function output 85 112 SINCOS(zi, sin_alpha, cos_alpha); 86 double yyy = stacked_disks_kernel(q, 87 halfheight, 88 thick_layer, 113 double yyy = _kernel(q, 89 114 radius, 90 n_stacking,91 sigma_dnn,92 115 core_sld, 93 116 layer_sld, 94 117 solvent_sld, 118 halfheight, 119 thick_layer, 95 120 sin_alpha, 96 121 cos_alpha, 97 d); 122 sigma_dnn, 123 d, 124 n_stacking); 98 125 summ += Gauss76Wt[i] * yyy * sin_alpha; 99 126 } … … 105 132 } 106 133 107 static double form_volume( 108 double thick_core, 109 double thick_layer, 110 double radius, 111 double fp_n_stacking) 112 { 113 int n_stacking = (int)(fp_n_stacking + 0.5); 134 double form_volume(double thick_core, 135 double thick_layer, 136 double radius, 137 double n_stacking){ 114 138 double d = 2.0 * thick_layer + thick_core; 115 139 return M_PI * radius * radius * d * n_stacking; 116 140 } 117 141 118 static double Iq( 119 double q, 120 double thick_core, 121 double thick_layer, 122 double radius, 123 double fp_n_stacking, 124 double sigma_dnn, 125 double core_sld, 126 double layer_sld, 127 double solvent_sld) 128 { 129 int n_stacking = (int)(fp_n_stacking + 0.5); 130 return stacked_disks_1d(q, 142 double Iq(double q, 143 double thick_core, 144 double thick_layer, 145 double radius, 146 double n_stacking, 147 double sigma_dnn, 148 double core_sld, 149 double layer_sld, 150 double solvent_sld) 151 { 152 return stacked_disks_kernel(q, 131 153 thick_core, 132 154 thick_layer, … … 140 162 141 163 142 static double Iqxy(double qx, double qy, 143 double thick_core,144 double thick_layer,145 double radius,146 double fp_n_stacking,147 double sigma_dnn,148 double core_sld,149 double layer_sld,150 double solvent_sld,151 double theta,152 double phi)153 { 154 int n_stacking = (int)(fp_n_stacking + 0.5); 164 double 165 Iqxy(double qx, double qy, 166 double thick_core, 167 double thick_layer, 168 double radius, 169 double n_stacking, 170 double sigma_dnn, 171 double core_sld, 172 double layer_sld, 173 double solvent_sld, 174 double theta, 175 double phi) 176 { 155 177 double q, sin_alpha, cos_alpha; 156 178 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); … … 158 180 double d = 2.0 * thick_layer + thick_core; 159 181 double halfheight = 0.5*thick_core; 160 double answer = stacked_disks_kernel(q, 161 halfheight, 162 thick_layer, 182 double answer = _kernel(q, 163 183 radius, 164 n_stacking,165 sigma_dnn,166 184 core_sld, 167 185 layer_sld, 168 186 solvent_sld, 187 halfheight, 188 thick_layer, 169 189 sin_alpha, 170 190 cos_alpha, 171 d); 191 sigma_dnn, 192 d, 193 n_stacking); 172 194 173 195 //convert to [cm-1] -
sasmodels/models/stacked_disks.py
rb7e8b94 rb7e8b94 126 126 ["thick_layer", "Ang", 10.0, [0, inf], "volume", "Thickness of layer each side of core"], 127 127 ["radius", "Ang", 15.0, [0, inf], "volume", "Radius of the stacked disk"], 128 ["n_stacking", "", 1.0, [ 1, inf], "volume", "Number of stacked layer/core/layer disks"],128 ["n_stacking", "", 1.0, [0, inf], "volume", "Number of stacked layer/core/layer disks"], 129 129 ["sigma_d", "Ang", 0, [0, inf], "", "Sigma of nearest neighbor spacing"], 130 130 ["sld_core", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Core scattering length density"], -
sasmodels/models/star_polymer.c
r2586093f r3a48772 3 3 double Iq(double q, double radius2, double arms); 4 4 5 static double star_polymer_kernel(double q, double radius2, double arms)5 static double _mass_fractal_kernel(double q, double radius2, double arms) 6 6 { 7 7 … … 23 23 double Iq(double q, double radius2, double arms) 24 24 { 25 return star_polymer_kernel(q, radius2, arms);25 return _mass_fractal_kernel(q, radius2, arms); 26 26 } -
sasmodels/models/unified_power_Rg.py
r66ca2a6 r66ca2a6 97 97 98 98 def Iq(q, level, rg, power, B, G): 99 level = int(level + 0.5)100 if level == 0:99 ilevel = int(level) 100 if ilevel == 0: 101 101 with errstate(divide='ignore'): 102 102 return 1./q … … 104 104 with errstate(divide='ignore', invalid='ignore'): 105 105 result = np.zeros(q.shape, 'd') 106 for i in range( level):106 for i in range(ilevel): 107 107 exp_now = exp(-(q*rg[i])**2/3.) 108 108 pow_now = (erf(q*rg[i]/sqrt(6.))**3/q)**power[i] 109 if i < level-1:109 if i < ilevel-1: 110 110 exp_next = exp(-(q*rg[i+1])**2/3.) 111 111 else: … … 113 113 result += G[i]*exp_now + B[i]*exp_next*pow_now 114 114 115 result[q == 0] = np.sum(G[: level])115 result[q == 0] = np.sum(G[:ilevel]) 116 116 return result 117 117
Note: See TracChangeset
for help on using the changeset viewer.