[f94d8a2] | 1 | /* |
---|
| 2 | Functions for WRC implementation of flexible cylinders |
---|
| 3 | */ |
---|
[ba32cdd] | 4 | |
---|
[f94d8a2] | 5 | static double |
---|
[6e72989] | 6 | Rgsquare(double L, double b) |
---|
[f94d8a2] | 7 | { |
---|
[6e72989] | 8 | const double x = L/b; |
---|
[237c9cf] | 9 | // Use Horner's method to evaluate Pedersen eq 15: |
---|
| 10 | // alpha^2 = [1.0 + (x/3.12)^2 + (x/8.67)^3] ^ (0.176/3) |
---|
[6e72989] | 11 | const double alphasq = |
---|
[237c9cf] | 12 | pow(1.0 + x*x*(1.027284681130835e-01 + 1.534414548417740e-03*x), |
---|
[6e72989] | 13 | 5.866666666666667e-02); |
---|
| 14 | return alphasq*L*b/6.0; |
---|
[f94d8a2] | 15 | } |
---|
| 16 | |
---|
| 17 | static double |
---|
[6e72989] | 18 | Rgsquareshort(double L, double b) |
---|
[f94d8a2] | 19 | { |
---|
[237c9cf] | 20 | const double r = b/L; // = 1/n_b in Pedersen ref. |
---|
[6e72989] | 21 | return Rgsquare(L, b) * (1.0 + r*(-1.5 + r*(1.5 + r*0.75*expm1(-2.0/r)))); |
---|
[f94d8a2] | 22 | } |
---|
| 23 | |
---|
| 24 | static double |
---|
[237c9cf] | 25 | w_WR(double x) |
---|
| 26 | { |
---|
| 27 | // Pedersen eq. 16: |
---|
[5772737] | 28 | // w = [1 + tanh((x-C4)/C5)]/2 |
---|
[237c9cf] | 29 | const double C4 = 1.523; |
---|
| 30 | const double C5 = 0.1477; |
---|
| 31 | return 0.5 + 0.5*tanh((x - C4)/C5); |
---|
| 32 | } |
---|
| 33 | |
---|
| 34 | static double |
---|
| 35 | Sdebye(double qsq) |
---|
| 36 | { |
---|
| 37 | #if FLOAT_SIZE>4 |
---|
| 38 | # define DEBYE_CUTOFF 0.25 // 1e-15 error |
---|
| 39 | #else |
---|
| 40 | # define DEBYE_CUTOFF 1.0 // 4e-7 error |
---|
| 41 | #endif |
---|
| 42 | |
---|
| 43 | if (qsq < DEBYE_CUTOFF) { |
---|
| 44 | const double x = qsq; |
---|
| 45 | // mathematica: PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] |
---|
| 46 | const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; |
---|
| 47 | const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.; |
---|
| 48 | return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.) |
---|
| 49 | / ((((B4*x + B3)*x + B2)*x + B1)*x + 1.); |
---|
| 50 | } else { |
---|
| 51 | return 2.*(expm1(-qsq) + qsq)/(qsq*qsq); |
---|
| 52 | } |
---|
| 53 | } |
---|
| 54 | |
---|
| 55 | static double |
---|
[a8631ca] | 56 | a_long(double q, double L, double b/*, double p1, double p2, double q0*/) |
---|
[f94d8a2] | 57 | { |
---|
[6e72989] | 58 | const double p1 = 4.12; |
---|
| 59 | const double p2 = 4.42; |
---|
| 60 | const double q0 = 3.1; |
---|
[f94d8a2] | 61 | |
---|
[237c9cf] | 62 | // Constants C1, ..., C5 derived from least squares fit (Pedersen, eq 13,16) |
---|
[e7678b2] | 63 | const double C1 = 1.22; |
---|
| 64 | const double C2 = 0.4288; |
---|
| 65 | const double C3 = -1.651; |
---|
| 66 | const double C4 = 1.523; |
---|
| 67 | const double C5 = 0.1477; |
---|
| 68 | const double miu = 0.585; |
---|
| 69 | |
---|
[6e72989] | 70 | const double C = (L/b>10.0 ? 3.06*pow(L/b, -0.44) : 1.0); |
---|
[a8631ca] | 71 | //printf("branch B-%d q=%g L=%g b=%g\n", C==1.0, q, L, b); |
---|
[6e72989] | 72 | const double r2 = Rgsquare(L,b); |
---|
| 73 | const double r = sqrt(r2); |
---|
| 74 | const double qr_b = q0*r/b; |
---|
| 75 | const double qr_b_sq = qr_b*qr_b; |
---|
| 76 | const double qr_b_4 = qr_b_sq*qr_b_sq; |
---|
| 77 | const double qr_b_miu = pow(qr_b, -1.0/miu); |
---|
| 78 | const double em1_qr_b_sq = expm1(-qr_b_sq); |
---|
| 79 | const double sech2 = 1.0/square(cosh((qr_b-C4)/C5)); |
---|
[237c9cf] | 80 | const double w = w_WR(qr_b); |
---|
[6e72989] | 81 | |
---|
| 82 | const double t1 = pow(q0, 1.0 + p1 + p2)/(b*(p1-p2)); |
---|
| 83 | const double t2 = C/(15.0*L) * ( |
---|
| 84 | + 14.0*b*b*em1_qr_b_sq/(q0*qr_b_sq) |
---|
| 85 | + 2.0*q0*r2*exp(-qr_b_sq)*(11.0 + 7.0/qr_b_sq)); |
---|
| 86 | const double t11 = ((C3*qr_b_miu + C2)*qr_b_miu + C1)*qr_b_miu; |
---|
| 87 | const double t3 = r*sech2/(2.*C5)*t11; |
---|
| 88 | const double t4 = r*(em1_qr_b_sq + qr_b_sq)*sech2 / (C5*qr_b_4); |
---|
[237c9cf] | 89 | const double t5 = -4.0*r*qr_b*em1_qr_b_sq/qr_b_4 * (1.0 - w); |
---|
| 90 | const double t10 = 2.0*(em1_qr_b_sq + qr_b_sq)/qr_b_4 * (1.0 - w); //=Sdebye*(1-w) |
---|
[6e72989] | 91 | const double t6 = 4.0*b/q0 * t10; |
---|
| 92 | const double t7 = r*((-3.0*C3*qr_b_miu -2.0*C2)*qr_b_miu -1.0*C1)*qr_b_miu/(miu*qr_b); |
---|
[237c9cf] | 93 | const double t9 = C*b/L * (4.0 - exp(-qr_b_sq) * (11.0 + 7.0/qr_b_sq) + 7.0/qr_b_sq)/15.0; |
---|
| 94 | const double t12 = b*b*M_PI/(L*q0*q0) + t2 + t3 - t4 + t5 - t6 + t7*w; |
---|
| 95 | const double t13 = -b*M_PI/(L*q0) + t9 + t10 + t11*w; |
---|
[6e72989] | 96 | |
---|
| 97 | const double a1 = pow(q0,p1)*t13 - t1*pow(q0,-p2)*(t12 + b*p1/q0*t13); |
---|
| 98 | const double a2 = t1*pow(q0,-p1)*(t12 + b*p1/q0*t13); |
---|
| 99 | |
---|
[a8631ca] | 100 | const double ans = a1*pow(q*b, -p1) + a2*pow(q*b, -p2) + M_PI/(q*L); |
---|
[6e72989] | 101 | return ans; |
---|
| 102 | } |
---|
| 103 | |
---|
| 104 | static double |
---|
| 105 | _short(double r2, double exp_qr_b, double L, double b, double p1short, |
---|
| 106 | double p2short, double q0) |
---|
| 107 | { |
---|
| 108 | const double qr2 = q0*q0 * r2; |
---|
[e7678b2] | 109 | const double b3 = b*b*b; |
---|
[6e72989] | 110 | const double q0p = pow(q0, -4.0 + p1short); |
---|
[f94d8a2] | 111 | |
---|
[6e72989] | 112 | double yy = 1.0/(L*r2*r2) * b/exp_qr_b*q0p |
---|
| 113 | * (8.0*b3*L |
---|
| 114 | - 8.0*b3*exp_qr_b*L |
---|
| 115 | + 2.0*b3*exp_qr_b*L*p2short |
---|
| 116 | - 2.0*b*exp_qr_b*L*p2short*qr2 |
---|
| 117 | + 4.0*b*exp_qr_b*L*qr2 |
---|
| 118 | - 2.0*b3*L*p2short |
---|
| 119 | + 4.0*b*L*qr2 |
---|
| 120 | - M_PI*exp_qr_b*qr2*q0*r2 |
---|
| 121 | + M_PI*exp_qr_b*p2short*qr2*q0*r2); |
---|
[f94d8a2] | 122 | |
---|
[6e72989] | 123 | return yy; |
---|
[f94d8a2] | 124 | } |
---|
| 125 | static double |
---|
[6e72989] | 126 | a_short(double qp, double L, double b |
---|
| 127 | /*double p1short, double p2short*/, double q0) |
---|
[f94d8a2] | 128 | { |
---|
[6e72989] | 129 | const double p1short = 5.36; |
---|
| 130 | const double p2short = 5.62; |
---|
[f94d8a2] | 131 | |
---|
[6e72989] | 132 | const double r2 = Rgsquareshort(L,b); |
---|
| 133 | const double exp_qr_b = exp(r2*square(q0/b)); |
---|
| 134 | const double pdiff = p1short - p2short; |
---|
| 135 | const double a1 = _short(r2,exp_qr_b,L,b,p1short,p2short,q0)/pdiff; |
---|
| 136 | const double a2= -_short(r2,exp_qr_b,L,b,p2short,p1short,q0)/pdiff; |
---|
| 137 | const double ans = a1*pow(qp*b, -p1short) + a2*pow(qp*b, -p2short) + M_PI/(qp*L); |
---|
| 138 | return ans; |
---|
[f94d8a2] | 139 | } |
---|
| 140 | |
---|
| 141 | |
---|
| 142 | static double |
---|
| 143 | Sexv(double q, double L, double b) |
---|
| 144 | { |
---|
[237c9cf] | 145 | // Pedersen eq 13, corrected by Chen eq A.5, swapping w and 1-w |
---|
[e7678b2] | 146 | const double C1=1.22; |
---|
| 147 | const double C2=0.4288; |
---|
| 148 | const double C3=-1.651; |
---|
| 149 | const double miu = 0.585; |
---|
[6e72989] | 150 | const double qr = q*sqrt(Rgsquare(L,b)); |
---|
[237c9cf] | 151 | const double qr_miu = pow(qr, -1.0/miu); |
---|
[6e72989] | 152 | const double w = w_WR(qr); |
---|
[5772737] | 153 | const double t10 = Sdebye(qr*qr)*(1.0 - w); |
---|
[237c9cf] | 154 | const double t11 = ((C3*qr_miu + C2)*qr_miu + C1)*qr_miu; |
---|
[f94d8a2] | 155 | |
---|
[237c9cf] | 156 | return t10 + w*t11; |
---|
[f94d8a2] | 157 | } |
---|
| 158 | |
---|
[237c9cf] | 159 | // Modified by Yun on Oct. 15, |
---|
[f94d8a2] | 160 | static double |
---|
[6e72989] | 161 | Sexv_new(double q, double L, double b) |
---|
[f94d8a2] | 162 | { |
---|
[6e72989] | 163 | const double qr = q*sqrt(Rgsquare(L,b)); |
---|
| 164 | const double qr2 = qr*qr; |
---|
[237c9cf] | 165 | const double C = (L/b > 10.0) ? 3.06*pow(L/b, -0.44) : 1.0; |
---|
| 166 | const double t9 = C*b/L * (4.0 - exp(-qr2) * (11.0 + 7.0/qr2) + 7.0/qr2)/15.0; |
---|
[f94d8a2] | 167 | |
---|
[6e72989] | 168 | const double Sexv_orig = Sexv(q, L, b); |
---|
[f94d8a2] | 169 | |
---|
[6e72989] | 170 | // calculating the derivative to decide on the correction (cutoff) term? |
---|
| 171 | // Note: this is modified from WRs original code |
---|
| 172 | const double del=1.05; |
---|
| 173 | const double qdel = (Sexv(q*del,L,b) - Sexv_orig)/(q*(del - 1.0)); |
---|
[f94d8a2] | 174 | |
---|
[6e72989] | 175 | if (qdel < 0) { |
---|
[a8631ca] | 176 | //printf("branch A1-%d q=%g L=%g b=%g\n", C==1.0, q, L, b); |
---|
[237c9cf] | 177 | return t9 + Sexv_orig; |
---|
[6e72989] | 178 | } else { |
---|
[a8631ca] | 179 | //printf("branch A2-%d q=%g L=%g b=%g\n", C==1.0, q, L, b); |
---|
[6e72989] | 180 | const double w = w_WR(qr); |
---|
[237c9cf] | 181 | const double t10 = Sdebye(qr*qr)*(1.0 - w); |
---|
| 182 | return t9 + t10; |
---|
[6e72989] | 183 | } |
---|
[f94d8a2] | 184 | } |
---|
| 185 | |
---|
| 186 | |
---|
[6e72989] | 187 | static double |
---|
| 188 | Sk_WR(double q, double L, double b) |
---|
| 189 | { |
---|
| 190 | const double Rg_short = sqrt(Rgsquareshort(L, b)); |
---|
| 191 | double q0short = fmax(1.9/Rg_short, 3.0); |
---|
| 192 | double ans; |
---|
[f94d8a2] | 193 | |
---|
[18a2bfc] | 194 | |
---|
[6e72989] | 195 | if( L > 4*b ) { // L > 4*b : Longer Chains |
---|
| 196 | if (q*b <= 3.1) { |
---|
| 197 | ans = Sexv_new(q, L, b); |
---|
| 198 | } else { //q(i)*b > 3.1 |
---|
| 199 | ans = a_long(q, L, b /*, p1, p2, q0*/); |
---|
| 200 | } |
---|
| 201 | } else { // L <= 4*b : Shorter Chains |
---|
| 202 | if (q*b <= q0short) { // q*b <= fmax(1.9/Rg_short, 3) |
---|
[a8631ca] | 203 | //printf("branch C-%d q=%g L=%g b=%g\n", square(q*Rg_short)<DEBYE_CUTOFF, q, L, b); |
---|
[18a2bfc] | 204 | // Note that q0short is usually 3, but it will be greater than 3 |
---|
| 205 | // small enough b, depending on the L/b ratio: |
---|
| 206 | // L/b == 1 => b < 2.37 |
---|
| 207 | // L/b == 2 => b < 1.36 |
---|
| 208 | // L/b == 3 => b < 1.00 |
---|
| 209 | // L/b == 4 => b < 0.816 |
---|
| 210 | // 2017-10-01 pkienzle: moved low q approximation into Sdebye() |
---|
[237c9cf] | 211 | ans = Sdebye(square(q*Rg_short)); |
---|
[6e72989] | 212 | } else { // q*b > max(1.9/Rg_short, 3) |
---|
[a8631ca] | 213 | //printf("branch D q=%g L=%g b=%g\n", q, L, b); |
---|
[6e72989] | 214 | ans = a_short(q, L, b /*, p1short, p2short*/, q0short); |
---|
| 215 | } |
---|
[f94d8a2] | 216 | } |
---|
| 217 | |
---|
[6e72989] | 218 | return ans; |
---|
[f94d8a2] | 219 | } |
---|