Ignore:
Timestamp:
Jan 5, 2012 12:16:29 PM (13 years ago)
Author:
Mathieu Doucet <doucetm@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
011e0e4
Parents:
bbbed8c
Message:

refactored bunch of models

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/c_models/pearlnecklace.cpp

    r67424cd r0ba3b08  
    11 
    22#include <math.h> 
    3 #include "models.hh" 
    43#include "parameters.hh" 
    54#include <stdio.h> 
    65using namespace std; 
     6#include "pearlnecklace.h" 
    77 
    88extern "C" { 
    9         #include "pearlnecklace.h" 
     9#include "libmultifunc/libfunc.h" 
     10} 
     11 
     12static double pearl_necklace_kernel(double dp[], double q) { 
     13  // fit parameters 
     14  double scale = dp[0]; 
     15  double radius = dp[1]; 
     16  double edge_separation = dp[2]; 
     17  double thick_string = dp[3]; 
     18  double num_pearls = dp[4]; 
     19  double sld_pearl = dp[5]; 
     20  double sld_string = dp[6]; 
     21  double sld_solv = dp[7]; 
     22  double background = dp[8]; 
     23 
     24  //relative slds 
     25  double contrast_pearl = sld_pearl - sld_solv; 
     26  double contrast_string = sld_string - sld_solv; 
     27 
     28  // number of string segments 
     29  double num_strings = num_pearls - 1.0; 
     30 
     31  //Pi 
     32  double pi = 4.0*atan(1.0); 
     33 
     34  // each volumes 
     35  double string_vol = edge_separation * pi * pow((thick_string / 2.0), 2); 
     36  double pearl_vol = 4.0 /3.0 * pi * pow(radius, 3); 
     37 
     38  //total volume 
     39  double tot_vol; 
     40  //each masses 
     41  double m_r= contrast_string * string_vol; 
     42  double m_s= contrast_pearl * pearl_vol; 
     43  double psi; 
     44  double gamma; 
     45  double beta; 
     46  //form factors 
     47  double sss; //pearls 
     48  double srr; //strings 
     49  double srs; //cross 
     50  double A_s; 
     51  double srr_1; 
     52  double srr_2; 
     53  double srr_3; 
     54  double form_factor; 
     55 
     56  tot_vol = num_strings * string_vol; 
     57  tot_vol += num_pearls * pearl_vol; 
     58 
     59  //sine functions of a pearl 
     60  psi = sin(q*radius); 
     61  psi -= q * radius * cos(q * radius); 
     62  psi /= pow((q * radius), 3); 
     63  psi *= 3.0; 
     64 
     65  // Note take only 20 terms in Si series: 10 terms may be enough though. 
     66  gamma = Si(q* edge_separation); 
     67  gamma /= (q* edge_separation); 
     68  beta = Si(q * (edge_separation + radius)); 
     69  beta -= Si(q * radius); 
     70  beta /= (q* edge_separation); 
     71 
     72  // center to center distance between the neighboring pearls 
     73  A_s = edge_separation + 2.0 * radius; 
     74 
     75  // form factor for num_pearls 
     76  sss = 1.0 - pow(sinc(q*A_s), num_pearls ); 
     77  sss /= pow((1.0-sinc(q*A_s)), 2); 
     78  sss *= -sinc(q*A_s); 
     79  sss -= num_pearls/2.0; 
     80  sss += num_pearls/(1.0-sinc(q*A_s)); 
     81  sss *= 2.0 * pow((m_s*psi), 2); 
     82 
     83  // form factor for num_strings (like thin rods) 
     84  srr_1 = -pow(sinc(q*edge_separation/2.0), 2); 
     85 
     86  srr_1 += 2.0 * gamma; 
     87  srr_1 *= num_strings; 
     88  srr_2 = 2.0/(1.0-sinc(q*A_s)); 
     89  srr_2 *= num_strings; 
     90  srr_2 *= pow(beta, 2); 
     91  srr_3 = 1.0 - pow(sinc(q*A_s), num_strings); 
     92  srr_3 /= pow((1.0-sinc(q*A_s)), 2); 
     93  srr_3 *= pow(beta, 2); 
     94  srr_3 *= -2.0; 
     95 
     96  // total srr 
     97  srr = srr_1 + srr_2 + srr_3; 
     98  srr *= pow(m_r, 2); 
     99 
     100  // form factor for correlations 
     101  srs = 1.0; 
     102  srs -= pow(sinc(q*A_s), num_strings); 
     103  srs /= pow((1.0-sinc(q*A_s)), 2); 
     104  srs *= -sinc(q*A_s); 
     105  srs += (num_strings/(1.0-sinc(q*A_s))); 
     106  srs *= 4.0; 
     107  srs *= (m_r * m_s * beta * psi); 
     108 
     109  form_factor = sss + srr + srs; 
     110  form_factor /= (tot_vol * 1.0e-8); // norm by volume and A^-1 to cm^-1 
     111 
     112  // scale and background 
     113  form_factor *= scale; 
     114  form_factor += background; 
     115  return (form_factor); 
    10116} 
    11117 
    12118PearlNecklaceModel :: PearlNecklaceModel() { 
    13         scale = Parameter(1.0); 
    14         radius = Parameter(80.0, true); 
    15         radius.set_min(0.0); 
    16         edge_separation = Parameter(350.0, true); 
    17         edge_separation.set_min(0.0); 
    18         thick_string = Parameter(2.5, true); 
    19         thick_string.set_min(0.0); 
    20         num_pearls = Parameter(3); 
    21         num_pearls.set_min(0.0); 
    22         sld_pearl = Parameter(1.0e-06); 
    23         sld_string = Parameter(1.0e-06); 
    24         sld_solv = Parameter(6.3e-06); 
    25     background = Parameter(0.0); 
     119  scale = Parameter(1.0); 
     120  radius = Parameter(80.0, true); 
     121  radius.set_min(0.0); 
     122  edge_separation = Parameter(350.0, true); 
     123  edge_separation.set_min(0.0); 
     124  thick_string = Parameter(2.5, true); 
     125  thick_string.set_min(0.0); 
     126  num_pearls = Parameter(3); 
     127  num_pearls.set_min(0.0); 
     128  sld_pearl = Parameter(1.0e-06); 
     129  sld_string = Parameter(1.0e-06); 
     130  sld_solv = Parameter(6.3e-06); 
     131  background = Parameter(0.0); 
    26132 
    27133} 
    … …  
    33139 */ 
    34140double PearlNecklaceModel :: operator()(double q) { 
    35         double dp[9]; 
    36         // Fill parameter array for IGOR library 
    37         // Add the background after averaging 
    38         dp[0] = scale(); 
    39         dp[1] = radius(); 
    40         dp[2] = edge_separation(); 
    41         dp[3] = thick_string(); 
    42         dp[4] = num_pearls(); 
    43         dp[5] = sld_pearl(); 
    44         dp[6] = sld_string(); 
    45         dp[7] = sld_solv(); 
    46         dp[8] = 0.0; 
    47         double pi = 4.0*atan(1.0); 
    48         // No polydispersion supported in this model. 
    49         // Get the dispersion points for the radius 
    50         vector<WeightPoint> weights_radius; 
    51         radius.get_weights(weights_radius); 
    52         vector<WeightPoint> weights_edge_separation; 
    53         edge_separation.get_weights(weights_edge_separation); 
    54         // Perform the computation, with all weight points 
    55         double sum = 0.0; 
    56         double norm = 0.0; 
    57         double vol = 0.0; 
    58         double string_vol = 0.0; 
    59         double pearl_vol = 0.0; 
    60         double tot_vol = 0.0; 
    61         // Loop over core weight points 
    62         for(size_t i=0; i<weights_radius.size(); i++) { 
    63                 dp[1] = weights_radius[i].value; 
    64                 // Loop over thick_inter0 weight points 
    65                 for(size_t j=0; j<weights_edge_separation.size(); j++) { 
    66                         dp[2] = weights_edge_separation[j].value; 
    67                         pearl_vol = 4.0 /3.0 * pi * pow(dp[1], 3); 
    68                         string_vol =dp[2] * pi * pow((dp[3] / 2.0), 2); 
    69                         tot_vol = (dp[4] - 1.0) * string_vol; 
    70                         tot_vol += dp[4] * pearl_vol; 
    71                         //Un-normalize Sphere by volume 
    72                         sum += weights_radius[i].weight * weights_edge_separation[j].weight 
    73                                 * pearl_necklace_kernel(dp,q) * tot_vol; 
    74                         //Find average volume 
    75                         vol += weights_radius[i].weight * weights_edge_separation[j].weight 
    76                                 * tot_vol; 
    77                         norm += weights_radius[i].weight * weights_edge_separation[j].weight; 
    78                 } 
    79         } 
    80  
    81         if (vol != 0.0 && norm != 0.0) { 
    82                 //Re-normalize by avg volume 
    83                 sum = sum/(vol/norm);} 
    84  
    85         return sum/norm + background(); 
     141  double dp[9]; 
     142  // Fill parameter array for IGOR library 
     143  // Add the background after averaging 
     144  dp[0] = scale(); 
     145  dp[1] = radius(); 
     146  dp[2] = edge_separation(); 
     147  dp[3] = thick_string(); 
     148  dp[4] = num_pearls(); 
     149  dp[5] = sld_pearl(); 
     150  dp[6] = sld_string(); 
     151  dp[7] = sld_solv(); 
     152  dp[8] = 0.0; 
     153  double pi = 4.0*atan(1.0); 
     154  // No polydispersion supported in this model. 
     155  // Get the dispersion points for the radius 
     156  vector<WeightPoint> weights_radius; 
     157  radius.get_weights(weights_radius); 
     158  vector<WeightPoint> weights_edge_separation; 
     159  edge_separation.get_weights(weights_edge_separation); 
     160  // Perform the computation, with all weight points 
     161  double sum = 0.0; 
     162  double norm = 0.0; 
     163  double vol = 0.0; 
     164  double string_vol = 0.0; 
     165  double pearl_vol = 0.0; 
     166  double tot_vol = 0.0; 
     167  // Loop over core weight points 
     168  for(size_t i=0; i<weights_radius.size(); i++) { 
     169    dp[1] = weights_radius[i].value; 
     170    // Loop over thick_inter0 weight points 
     171    for(size_t j=0; j<weights_edge_separation.size(); j++) { 
     172      dp[2] = weights_edge_separation[j].value; 
     173      pearl_vol = 4.0 /3.0 * pi * pow(dp[1], 3); 
     174      string_vol =dp[2] * pi * pow((dp[3] / 2.0), 2); 
     175      tot_vol = (dp[4] - 1.0) * string_vol; 
     176      tot_vol += dp[4] * pearl_vol; 
     177      //Un-normalize Sphere by volume 
     178      sum += weights_radius[i].weight * weights_edge_separation[j].weight 
     179          * pearl_necklace_kernel(dp,q) * tot_vol; 
     180      //Find average volume 
     181      vol += weights_radius[i].weight * weights_edge_separation[j].weight 
     182          * tot_vol; 
     183      norm += weights_radius[i].weight * weights_edge_separation[j].weight; 
     184    } 
     185  } 
     186 
     187  if (vol != 0.0 && norm != 0.0) { 
     188    //Re-normalize by avg volume 
     189    sum = sum/(vol/norm);} 
     190 
     191  return sum/norm + background(); 
    86192} 
    87193 
    … …  
    93199 */ 
    94200double PearlNecklaceModel :: operator()(double qx, double qy) { 
    95         double q = sqrt(qx*qx + qy*qy); 
    96         return (*this).operator()(q); 
     201  double q = sqrt(qx*qx + qy*qy); 
     202  return (*this).operator()(q); 
    97203} 
    98204 
    … …  
    105211 */ 
    106212double PearlNecklaceModel :: evaluate_rphi(double q, double phi) { 
    107         return (*this).operator()(q); 
     213  return (*this).operator()(q); 
    108214} 
    109215 
    … …  
    118224// sld of solvent. 
    119225double PearlNecklaceModel :: calculate_ER() { 
    120         PeralNecklaceParameters dp; 
    121  
    122         dp.scale = scale(); 
    123         dp.radius = radius(); 
    124         dp.edge_separation = edge_separation(); 
    125         dp.thick_string = thick_string(); 
    126         dp.num_pearls = num_pearls(); 
    127  
    128         double rad_out = 0.0; 
    129         // Perform the computation, with all weight points 
    130         double sum = 0.0; 
    131         double norm = 0.0; 
    132         double pi = 4.0*atan(1.0); 
    133         // No polydispersion supported in this model. 
    134         // Get the dispersion points for the radius 
    135         vector<WeightPoint> weights_radius; 
    136         radius.get_weights(weights_radius); 
    137         vector<WeightPoint> weights_edge_separation; 
    138         edge_separation.get_weights(weights_edge_separation); 
    139         // Perform the computation, with all weight points 
    140         double string_vol = 0.0; 
    141         double pearl_vol = 0.0; 
    142         double tot_vol = 0.0; 
    143         // Loop over core weight points 
    144         for(size_t i=0; i<weights_radius.size(); i++) { 
    145                 dp.radius = weights_radius[i].value; 
    146                 // Loop over thick_inter0 weight points 
    147                 for(size_t j=0; j<weights_edge_separation.size(); j++) { 
    148                         dp.edge_separation = weights_edge_separation[j].value; 
    149                         pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3); 
    150                         string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2); 
    151                         tot_vol = (dp.num_pearls - 1.0) * string_vol; 
    152                         tot_vol += dp.num_pearls * pearl_vol; 
    153                         //Find  volume 
    154                         // This may be a too much approximation 
    155                         //Todo: decided whether or not we keep this calculation 
    156                         sum += weights_radius[i].weight * weights_edge_separation[j].weight 
    157                                 * pow(3.0*tot_vol/4.0/pi,0.333333); 
    158                         norm += weights_radius[i].weight * weights_edge_separation[j].weight; 
    159                 } 
    160         } 
    161  
    162         if (norm != 0){ 
    163                 //return the averaged value 
    164                 rad_out =  sum/norm;} 
    165         else{ 
    166                 //return normal value 
    167                 pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3); 
    168                 string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2); 
    169                 tot_vol = (dp.num_pearls - 1.0) * string_vol; 
    170                 tot_vol += dp.num_pearls * pearl_vol; 
    171  
    172                 rad_out =  pow((3.0*tot_vol/4.0/pi), 0.33333); 
    173                 } 
    174  
    175         return rad_out; 
    176  
    177 } 
     226  PeralNecklaceParameters dp; 
     227 
     228  dp.scale = scale(); 
     229  dp.radius = radius(); 
     230  dp.edge_separation = edge_separation(); 
     231  dp.thick_string = thick_string(); 
     232  dp.num_pearls = num_pearls(); 
     233 
     234  double rad_out = 0.0; 
     235  // Perform the computation, with all weight points 
     236  double sum = 0.0; 
     237  double norm = 0.0; 
     238  double pi = 4.0*atan(1.0); 
     239  // No polydispersion supported in this model. 
     240  // Get the dispersion points for the radius 
     241  vector<WeightPoint> weights_radius; 
     242  radius.get_weights(weights_radius); 
     243  vector<WeightPoint> weights_edge_separation; 
     244  edge_separation.get_weights(weights_edge_separation); 
     245  // Perform the computation, with all weight points 
     246  double string_vol = 0.0; 
     247  double pearl_vol = 0.0; 
     248  double tot_vol = 0.0; 
     249  // Loop over core weight points 
     250  for(size_t i=0; i<weights_radius.size(); i++) { 
     251    dp.radius = weights_radius[i].value; 
     252    // Loop over thick_inter0 weight points 
     253    for(size_t j=0; j<weights_edge_separation.size(); j++) { 
     254      dp.edge_separation = weights_edge_separation[j].value; 
     255      pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3); 
     256      string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2); 
     257      tot_vol = (dp.num_pearls - 1.0) * string_vol; 
     258      tot_vol += dp.num_pearls * pearl_vol; 
     259      //Find  volume 
     260      // This may be a too much approximation 
     261      //Todo: decided whether or not we keep this calculation 
     262      sum += weights_radius[i].weight * weights_edge_separation[j].weight 
     263          * pow(3.0*tot_vol/4.0/pi,0.333333); 
     264      norm += weights_radius[i].weight * weights_edge_separation[j].weight; 
     265    } 
     266  } 
     267 
     268  if (norm != 0){ 
     269    //return the averaged value 
     270    rad_out =  sum/norm;} 
     271  else{ 
     272    //return normal value 
     273    pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3); 
     274    string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2); 
     275    tot_vol = (dp.num_pearls - 1.0) * string_vol; 
     276    tot_vol += dp.num_pearls * pearl_vol; 
     277 
     278    rad_out =  pow((3.0*tot_vol/4.0/pi), 0.33333); 
     279  } 
     280 
     281  return rad_out; 
     282 
     283} 
Note: See TracChangeset for help on using the changeset viewer.