Changeset c5ac2b2 in sasmodels for sasmodels/models/core_multi_shell.c


Ignore:
Timestamp:
Jul 18, 2016 2:05:12 AM (8 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
4a21b121
Parents:
256dfe1
Message:

working core multi shell

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/core_multi_shell.c

    rabdd01c rc5ac2b2  
    1 FourShell(double dp[], double q) 
    2 { 
    3     // variables are: 
    4     //[0] scale factor 
    5     //[1] radius of core [ï¿œ] 
    6     //[2] SLD of the core   [ï¿œ-2] 
    7     //[3] thickness of shell 1 [ï¿œ] 
    8     //[4] SLD of shell 1 
    9     //[5] thickness of shell 2 [ï¿œ] 
    10     //[6] SLD of shell 2 
    11     //[7] thickness of shell 3 
    12     //[8] SLD of shell 3 
    13     //[9] thickness of shell 3 
    14     //[10] SLD of shell 3 
    15     //[11] SLD of solvent 
    16     //[12] background   [cm-1] 
    17      
    18     double x,pi; 
    19     double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;     //my local names 
    20     double bes,f,vol,qr,contr,f2; 
    21     double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4; 
    22      
    23     pi = 4.0*atan(1.0); 
    24     x=q; 
    25      
    26     scale = dp[0]; 
    27     rcore = dp[1]; 
    28     rhocore = dp[2]; 
    29     thick1 = dp[3]; 
    30     rhoshel1 = dp[4]; 
    31     thick2 = dp[5]; 
    32     rhoshel2 = dp[6];    
    33     thick3 = dp[7]; 
    34     rhoshel3 = dp[8]; 
    35     thick4 = dp[9]; 
    36     rhoshel4 = dp[10];   
    37     rhosolv = dp[11]; 
    38     bkg = dp[12]; 
    39      
    40         // core first, then add in shells 
    41     qr=x*rcore; 
    42     contr = rhocore-rhoshel1; 
    43     if(qr == 0){ 
    44         bes = 1.0; 
    45     }else{ 
    46         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
    47     } 
    48     vol = 4.0*pi/3.0*rcore*rcore*rcore; 
    49     f = vol*bes*contr; 
    50     //now the shell (1) 
    51     qr=x*(rcore+thick1); 
    52     contr = rhoshel1-rhoshel2; 
    53     if(qr == 0){ 
    54         bes = 1.0; 
    55     }else{ 
    56         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
    57     } 
    58     vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 
    59     f += vol*bes*contr; 
    60     //now the shell (2) 
    61     qr=x*(rcore+thick1+thick2); 
    62     contr = rhoshel2-rhoshel3; 
    63     if(qr == 0){ 
    64         bes = 1.0; 
    65     }else{ 
    66         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
    67     } 
    68     vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 
    69     f += vol*bes*contr; 
    70     //now the shell (3) 
    71     qr=x*(rcore+thick1+thick2+thick3); 
    72     contr = rhoshel3-rhoshel4; 
    73     if(qr == 0){ 
    74         bes = 1.0; 
    75     }else{ 
    76         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
    77     } 
    78     vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); 
    79     f += vol*bes*contr; 
    80     //now the shell (4) 
    81     qr=x*(rcore+thick1+thick2+thick3+thick4); 
    82     contr = rhoshel4-rhosolv; 
    83     if(qr == 0){ 
    84         bes = 1.0; 
    85     }else{ 
    86         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
    87     } 
    88     vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4); 
    89     f += vol*bes*contr; 
    90      
    91          
    92     // normalize to particle volume and rescale from [ï¿œ-1] to [cm-1] 
    93     f2 = f*f/vol*1.0e8; 
    94      
    95     //scale if desired 
    96     f2 *= scale; 
    97     // then add in the background 
    98     f2 += bkg; 
    99      
    100     return(f2); 
    101 } 
    102  
    103 // above is cut and past with no changes from libigor four shell 
    104 static double 
    105 f_exp(double q, double r, double sld_in, double sld_out, 
    106     double thickness, double A) 
    107 { 
    108   const double vol = 4.0/3.0 * M_PI * r * r * r; 
    109   const double qr = q * r; 
    110   const double alpha = A * r/thickness; 
    111   const double bes = sph_j1c(qr); 
    112   const double B = (sld_out - sld_in)/expm1(A); 
    113   const double C = sld_in - B; 
    114   double fun; 
    115   if (qr == 0.0) { 
    116     fun = 1.0; 
    117   } else if (fabs(A) > 0.0) { 
    118     const double qrsq = qr*qr; 
    119     const double alphasq = alpha*alpha; 
    120     const double sumsq = alphasq + qrsq; 
    121     double sinqr, cosqr; 
    122     SINCOS(qr, sinqr, cosqr); 
    123     fun = -3.0( 
    124             ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq 
    125                 - (alpha*sinqr/qr - cosqr) 
    126         ) / sumsq; 
    127   } else { 
    128     fun = bes; 
    129   } 
    130   return vol * (B*fun + C*bes); 
    131 } 
    132  
    133 static double 
    134 f_linear(double q, double r, double sld, double slope) 
    135 { 
    136   const double vol = 4.0/3.0 * M_PI * r * r * r; 
    137   const double qr = q * r; 
    138   const double bes = sph_j1c(qr); 
    139   double fun = 0.0; 
    140   if (qr > 0.0) { 
    141     const double qrsq = qr*qr; 
    142     double sinqr, cosqr; 
    143     SINCOS(qr, sinqr, cosqr); 
    144     // Jae-He's code seems to simplify to this 
    145     //     fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq); 
    146     // Rederiving the math, we get the following instead: 
    147     fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq); 
    148   } 
    149   return vol * (sld*bes + fun); 
    150 } 
    1511 
    1522static double 
    … …  
    1544{ 
    1555  const double bes = sph_j1c(q * r); 
    156   const double vol = 4.0/3.0 * M_PI * r * r * r; 
     6  const double vol = M_4PI_3 * cube(r); 
    1577  return sld * vol * bes; 
    1588} 
    1599 
    160 static double 
     10double 
     11form_volume(double core_radius, double n, double thickness[]); 
     12double 
    16113form_volume(double core_radius, double n, double thickness[]) 
    16214{ 
    163   int i; 
    16415  double r = core_radius; 
    165   for (i=0; i < n; i++) { 
     16  for (int i=0; i < n; i++) { 
    16617    r += thickness[i]; 
    16718  } 
    168   return 4.0/3.0 * M_PI * r * r * r; 
     19  return M_4PI_3 * cube(r); 
    16920} 
    17021 
    171 static double 
    172 Iq(double q, double core_sld, double core_radius, double solvent_sld, 
    173     double n, double in_sld[], double out_sld[], double thickness[], 
    174     double A[]) 
     22double 
     23Iq(double q, double core_sld, double core_radius, 
     24   double solvent_sld, double num_shells, double sld[], double thickness[]); 
     25double 
     26Iq(double q, double core_sld, double core_radius, 
     27   double solvent_sld, double num_shells, double sld[], double thickness[]) 
    17528{ 
    176   int i; 
     29  const int n = (int)ceil(num_shells); 
     30  double f, r, last_sld; 
    17731  r = core_radius; 
    178   f = f_constant(q, r, core_sld); 
    179   for (i=0; i<n; i++){ 
    180     const double r0 = r; 
     32  last_sld = core_sld; 
     33  f = 0.; 
     34  for (int i=0; i<n; i++) { 
     35    f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * sph_j1c(q*r); 
     36    last_sld = sld[i]; 
    18137    r += thickness[i]; 
    182     if (r == r0) { 
    183       // no thickness, so nothing to add 
    184     } else if (fabs(A[i]) < 1.0e-16 || sld_out[i] == sld_in[i]) { 
    185       f -= f_constant(q, r0, sld_in[i]); 
    186       f += f_constant(q, r, sld_in[i]); 
    187     } else if (fabs(A[i]) < 1.0e-4) { 
    188       const double slope = (sld_out[i] - sld_in[i])/thickness[i]; 
    189       f -= f_linear(q, r0, sld_in[i], slope); 
    190       f += f_linear(q, r, sld_out[i], slope); 
    191     } else { 
    192       f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]); 
    193       f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]); 
    194     } 
    19538  } 
    196   f -= f_constant(q, r, solvent_sld); 
    197   f2 = f * f * 1.0e-4; 
    198  
    199   return f2; 
     39  f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sph_j1c(q*r); 
     40  return f * f * 1.0e-4; 
    20041} 
Note: See TracChangeset for help on using the changeset viewer.