Changeset 5bddd89 in sasmodels for sasmodels/models/hollow_cylinder.c
- Timestamp:
- Oct 14, 2016 4:05:39 PM (8 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 9068f4c
- Parents:
- 0b717c5
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/hollow_cylinder.c
r0d6e865 r5bddd89 1 1 double form_volume(double radius, double thickness, double length); 2 3 2 double Iq(double q, double radius, double thickness, double length, double sld, 4 3 double solvent_sld); … … 9 8 10 9 // From Igor library 11 static double hollow_cylinder_scaling(12 10 static double 11 _hollow_cylinder_scaling(double integrand, double delrho, double volume) 13 12 { 14 double answer; 15 // Multiply by contrast^2 16 answer = integrand*delrho*delrho; 17 18 //normalize by cylinder volume 19 answer *= volume*volume; 20 21 //convert to [cm-1] 22 answer *= 1.0e-4; 23 24 return answer; 13 return 1.0e-4 * square(volume * delrho * integrand); 25 14 } 26 15 27 16 28 static double _hollow_cylinder_kernel( 29 double q, double radius, double thickness, double length, double dum) 17 static double 18 _hollow_cylinder_kernel(double q, 19 double radius, double thickness, double length, double sin_val, double cos_val) 30 20 { 31 const double qs = q*s qrt(1.0-dum*dum);21 const double qs = q*sin_val; 32 22 const double lam1 = sas_J1c((radius+thickness)*qs); 33 23 const double lam2 = sas_J1c(radius*qs); … … 35 25 //Note: lim_{r -> r_c} psi = J0(radius_core*qs) 36 26 const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 37 const double t2 = sinc( q*length*dum/2.0);38 return square(psi*t2);27 const double t2 = sinc(0.5*q*length*cos_val); 28 return psi*t2; 39 29 } 40 30 41 42 static double hollow_cylinder_analytical_2D_scaled( 43 double q, double q_x, double q_y, double radius, double thickness, 44 double length, double sld, double solvent_sld, double theta, double phi) 45 { 46 double cyl_x, cyl_y; //, cyl_z 47 //double q_z; 48 double vol, cos_val, delrho; 49 double answer; 50 //convert angle degree to radian 51 theta = theta * M_PI_180; 52 phi = phi * M_PI_180; 53 delrho = solvent_sld - sld; 54 55 // Cylinder orientation 56 cyl_x = sin(theta) * cos(phi); 57 cyl_y = sin(theta) * sin(phi); 58 //cyl_z = -cos(theta) * sin(phi); 59 60 // q vector 61 //q_z = 0; 62 63 // Compute the angle btw vector q and the 64 // axis of the cylinder 65 cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z; 66 67 answer = _hollow_cylinder_kernel(q, radius, thickness, length, cos_val); 68 69 vol = form_volume(radius, thickness, length); 70 answer = hollow_cylinder_scaling(answer, delrho, vol); 71 72 return answer; 73 } 74 75 76 double form_volume(double radius, double thickness, double length) 31 double 32 form_volume(double radius, double thickness, double length) 77 33 { 78 34 double v_shell = M_PI*length*((radius+thickness)*(radius+thickness)-radius*radius); … … 81 37 82 38 83 double Iq(double q, double radius, double thickness, double length, 39 double 40 Iq(double q, double radius, double thickness, double length, 84 41 double sld, double solvent_sld) 85 42 { 86 int i; 87 double lower,upper,zi, inter; //upper and lower integration limits 88 double summ,answer,delrho; //running tally of integration 89 double norm,volume; //final calculation variables 43 const double lower = 0.0; 44 const double upper = 1.0; //limits of numerical integral 90 45 91 lower = 0.0; 92 upper = 1.0; //limits of numerical integral 93 94 summ = 0.0; //initialize intergral 95 for (i=0;i<76;i++) { 96 zi = ( Gauss76Z[i] * (upper-lower) + lower + upper )/2.0; 97 inter = Gauss76Wt[i] * _hollow_cylinder_kernel(q, radius, thickness, length, zi); 98 summ += inter; 46 double summ = 0.0; //initialize intergral 47 for (int i=0;i<76;i++) { 48 const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 49 const double sin_val = sqrt(1.0 - cos_val*cos_val); 50 const double inter = _hollow_cylinder_kernel(q, radius, thickness, length, 51 sin_val, cos_val); 52 summ += Gauss76Wt[i] * inter; 99 53 } 100 54 101 norm = summ*(upper-lower)/2.0;102 volume = form_volume(radius, thickness, length);103 delrho = solvent_sld - sld;104 answer = hollow_cylinder_scaling(norm, delrho, volume); 55 const double Aq = 0.5*summ*(upper-lower); 56 const double volume = form_volume(radius, thickness, length); 57 return _hollow_cylinder_scaling(Aq, solvent_sld - sld, volume); 58 } 105 59 106 return(answer); 60 double 61 Iqxy(double qx, double qy, 62 double radius, double thickness, double length, 63 double sld, double solvent_sld, double theta, double phi) 64 { 65 double q, sin_alpha, cos_alpha; 66 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 67 const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length, 68 sin_alpha, cos_alpha); 69 70 const double vol = form_volume(radius, thickness, length); 71 return _hollow_cylinder_scaling(Aq, solvent_sld-sld, vol); 107 72 } 108 73 109 74 110 double Iqxy(double qx, double qy, double radius, double thickness,111 double length, double sld, double solvent_sld, double theta, double phi)112 {113 const double q = sqrt(qx*qx+qy*qy);114 return hollow_cylinder_analytical_2D_scaled(q, qx/q, qy/q, radius, thickness, length, sld, solvent_sld, theta, phi);115 }
Note: See TracChangeset
for help on using the changeset viewer.