[a8b3cdb] | 1 | double form_volume(double r_minor, double r_ratio, double length); |
---|
| 2 | double Iq(double q, double r_minor, double r_ratio, double length, |
---|
| 3 | double sld, double solvent_sld); |
---|
| 4 | double Iqxy(double qx, double qy, double r_minor, double r_ratio, double length, |
---|
| 5 | double sld, double solvent_sld, double theta, double phi, double psi); |
---|
| 6 | |
---|
| 7 | |
---|
| 8 | double _elliptical_cylinder_kernel(double q, double r_minor, double r_ratio, double theta); |
---|
| 9 | |
---|
| 10 | double _elliptical_cylinder_kernel(double q, double r_minor, double r_ratio, double theta) |
---|
| 11 | { |
---|
| 12 | // This is the function LAMBDA1^2 in Feigin's notation |
---|
| 13 | // q is the q-value for the calculation (1/A) |
---|
| 14 | // r_minor is the transformed radius"a" in Feigin's notation |
---|
| 15 | // r_ratio is the ratio (major radius)/(minor radius) of the Ellipsoid [=] --- |
---|
| 16 | // theta is the dummy variable of the integration |
---|
| 17 | |
---|
| 18 | double retval,arg; |
---|
| 19 | |
---|
| 20 | arg = q*r_minor*sqrt((1.0+r_ratio*r_ratio)/2+(1.0-r_ratio*r_ratio)*cos(theta)/2); |
---|
| 21 | if (arg == 0.0){ |
---|
| 22 | retval = 1.0; |
---|
| 23 | }else{ |
---|
[43b7eea] | 24 | //retval = 2.0*NR_BessJ1(arg)/arg; |
---|
| 25 | retval = sas_J1c(arg); |
---|
[a8b3cdb] | 26 | } |
---|
| 27 | return retval*retval ; |
---|
| 28 | } |
---|
| 29 | |
---|
| 30 | |
---|
| 31 | double form_volume(double r_minor, double r_ratio, double length) |
---|
| 32 | { |
---|
| 33 | return M_PI * r_minor * r_minor * r_ratio * length; |
---|
| 34 | } |
---|
| 35 | |
---|
| 36 | double Iq(double q, double r_minor, double r_ratio, double length, |
---|
| 37 | double sld, double solvent_sld) { |
---|
| 38 | |
---|
| 39 | const int nordi=76; //order of integration |
---|
| 40 | const int nordj=20; |
---|
| 41 | double va,vb; //upper and lower integration limits |
---|
| 42 | double summ,zi,yyy,answer; //running tally of integration |
---|
| 43 | double summj,vaj,vbj,zij,arg,si; //for the inner integration |
---|
| 44 | |
---|
| 45 | // orientational average limits |
---|
| 46 | va = 0.0; |
---|
| 47 | vb = 1.0; |
---|
| 48 | // inner integral limits |
---|
| 49 | vaj=0.0; |
---|
| 50 | vbj=M_PI; |
---|
| 51 | |
---|
| 52 | //initialize integral |
---|
| 53 | summ = 0.0; |
---|
| 54 | |
---|
| 55 | const double delrho = sld - solvent_sld; |
---|
| 56 | |
---|
| 57 | for(int i=0;i<nordi;i++) { |
---|
| 58 | //setup inner integral over the ellipsoidal cross-section |
---|
| 59 | summj=0; |
---|
| 60 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy |
---|
| 61 | arg = r_minor*sqrt(1.0-zi*zi); |
---|
| 62 | for(int j=0;j<nordj;j++) { |
---|
| 63 | //20 gauss points for the inner integral |
---|
| 64 | zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy |
---|
| 65 | yyy = Gauss20Wt[j] * _elliptical_cylinder_kernel(q, arg, r_ratio, zij); |
---|
| 66 | summj += yyy; |
---|
| 67 | } |
---|
| 68 | //now calculate the value of the inner integral |
---|
| 69 | answer = (vbj-vaj)/2.0*summj; |
---|
| 70 | //divide integral by Pi |
---|
| 71 | answer /= M_PI; |
---|
| 72 | |
---|
| 73 | //now calculate outer integral |
---|
| 74 | arg = q*length*zi/2.0; |
---|
| 75 | if (arg == 0.0){ |
---|
| 76 | si = 1.0; |
---|
| 77 | }else{ |
---|
| 78 | si = sin(arg) * sin(arg) / arg / arg; |
---|
| 79 | } |
---|
| 80 | yyy = Gauss76Wt[i] * answer * si; |
---|
| 81 | summ += yyy; |
---|
| 82 | } |
---|
| 83 | |
---|
| 84 | answer = (vb-va)/2.0*summ; |
---|
| 85 | // Multiply by contrast^2 |
---|
| 86 | answer *= delrho*delrho; |
---|
| 87 | |
---|
| 88 | const double vol = form_volume(r_minor, r_ratio, length); |
---|
[abdd01c] | 89 | return answer*vol*vol*1.0e-4; |
---|
[a8b3cdb] | 90 | } |
---|
| 91 | |
---|
| 92 | |
---|
| 93 | double Iqxy(double qx, double qy, double r_minor, double r_ratio, double length, |
---|
| 94 | double sld, double solvent_sld, double theta, double phi, double psi) { |
---|
| 95 | const double _theta = theta * M_PI / 180.0; |
---|
| 96 | const double _phi = phi * M_PI / 180.0; |
---|
| 97 | const double _psi = psi * M_PI / 180.0; |
---|
| 98 | const double q = sqrt(qx*qx+qy*qy); |
---|
| 99 | const double q_x = qx/q; |
---|
| 100 | const double q_y = qy/q; |
---|
| 101 | |
---|
| 102 | //Cylinder orientation |
---|
| 103 | double cyl_x = cos(_theta) * cos(_phi); |
---|
| 104 | double cyl_y = sin(_theta); |
---|
| 105 | |
---|
| 106 | //cyl_z = -cos(_theta) * sin(_phi); |
---|
| 107 | |
---|
| 108 | // q vector |
---|
| 109 | //q_z = 0; |
---|
| 110 | |
---|
| 111 | // Note: cos(alpha) = 0 and 1 will get an |
---|
| 112 | // undefined value from CylKernel |
---|
| 113 | //alpha = acos( cos_val ); |
---|
| 114 | |
---|
| 115 | //ellipse orientation: |
---|
| 116 | // the elliptical corss section was transformed and projected |
---|
| 117 | // into the detector plane already through sin(alpha)and furthermore psi remains as same |
---|
| 118 | // on the detector plane. |
---|
| 119 | // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt |
---|
| 120 | // the wave vector q. |
---|
| 121 | |
---|
| 122 | //x- y- component on the detector plane. |
---|
| 123 | const double ella_x = -cos(_phi)*sin(_psi) * sin(_theta)+sin(_phi)*cos(_psi); |
---|
| 124 | const double ella_y = sin(_psi)*cos(_theta); |
---|
| 125 | const double ellb_x = -sin(_theta)*cos(_psi)*cos(_phi)-sin(_psi)*sin(_phi); |
---|
| 126 | const double ellb_y = cos(_theta)*cos(_psi); |
---|
| 127 | |
---|
| 128 | // Compute the angle btw vector q and the |
---|
| 129 | // axis of the cylinder |
---|
| 130 | double cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z; |
---|
| 131 | |
---|
| 132 | // calculate the axis of the ellipse wrt q-coord. |
---|
| 133 | double cos_nu = ella_x*q_x + ella_y*q_y; |
---|
| 134 | double cos_mu = ellb_x*q_x + ellb_y*q_y; |
---|
| 135 | |
---|
| 136 | // The following test should always pass |
---|
| 137 | if (fabs(cos_val)>1.0) { |
---|
| 138 | //printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n"); |
---|
| 139 | cos_val = 1.0; |
---|
| 140 | } |
---|
| 141 | if (fabs(cos_nu)>1.0) { |
---|
| 142 | //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n"); |
---|
| 143 | cos_nu = 1.0; |
---|
| 144 | } |
---|
| 145 | if (fabs(cos_mu)>1.0) { |
---|
| 146 | //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n"); |
---|
| 147 | cos_mu = 1.0; |
---|
| 148 | } |
---|
| 149 | |
---|
| 150 | const double r_major = r_ratio * r_minor; |
---|
| 151 | const double qr = q*sqrt( r_major*r_major*cos_nu*cos_nu + r_minor*r_minor*cos_mu*cos_mu ); |
---|
| 152 | const double qL = q*length*cos_val/2.0; |
---|
| 153 | |
---|
| 154 | double Be; |
---|
| 155 | if (qr==0){ |
---|
| 156 | Be = 0.5; |
---|
| 157 | }else{ |
---|
[43b7eea] | 158 | //Be = NR_BessJ1(qr)/qr; |
---|
| 159 | Be = 0.5*sas_J1c(qr); |
---|
[a8b3cdb] | 160 | } |
---|
| 161 | |
---|
| 162 | double Si; |
---|
| 163 | if (qL==0){ |
---|
| 164 | Si = 1.0; |
---|
| 165 | }else{ |
---|
| 166 | Si = sin(qL)/qL; |
---|
| 167 | } |
---|
| 168 | |
---|
| 169 | const double k = 2.0 * Be * Si; |
---|
| 170 | const double vol = form_volume(r_minor, r_ratio, length); |
---|
| 171 | return (sld - solvent_sld) * (sld - solvent_sld) * k * k *vol*vol*1.0e-4; |
---|
| 172 | } |
---|