Changeset 011e0e4 in sasview for sansmodels/src/c_models


Ignore:
Timestamp:
Jan 5, 2012 2:23:15 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:
98fdccd
Parents:
0ba3b08
Message:

core-shell + ellipsoid refactor

Location:
sansmodels/src/c_models
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/c_models/corefourshell.cpp

    r67424cd r011e0e4  
    2121 
    2222#include <math.h> 
    23 #include "models.hh" 
    2423#include "parameters.hh" 
    2524#include <stdio.h> 
    2625using namespace std; 
     26#include "corefourshell.h" 
    2727 
    2828extern "C" { 
    2929        #include "libSphere.h" 
    30         #include "corefourshell.h" 
    31 } 
     30} 
     31 
     32typedef struct { 
     33  double scale; 
     34  double rad_core0; 
     35  double sld_core0; 
     36  double thick_shell1; 
     37  double sld_shell1; 
     38  double thick_shell2; 
     39  double sld_shell2; 
     40  double thick_shell3; 
     41  double sld_shell3; 
     42  double thick_shell4; 
     43  double sld_shell4; 
     44  double sld_solv; 
     45  double background; 
     46} CoreFourShellParameters; 
    3247 
    3348CoreFourShellModel :: CoreFourShellModel() { 

  • sansmodels/src/c_models/coreshellcylinder.cpp

    r67424cd r011e0e4  
    1818 *   sansmodels/src/libigor 
    1919 * 
    20  *      TODO: refactor so that we pull in the old sansmodels.c_extensions 
    2120 */ 
    2221 
    2322#include <math.h> 
    24 #include "models.hh" 
    2523#include "parameters.hh" 
    2624#include <stdio.h> 
    2725using namespace std; 
     26#include "core_shell_cylinder.h" 
    2827 
    2928extern "C" { 
    3029        #include "libCylinder.h" 
    3130        #include "libStructureFactor.h" 
    32         #include "core_shell_cylinder.h" 
    33 } 
     31} 
     32 
     33typedef struct { 
     34  double scale; 
     35  double radius; 
     36  double thickness; 
     37  double length; 
     38  double core_sld; 
     39  double shell_sld; 
     40  double solvent_sld; 
     41  double background; 
     42  double axis_theta; 
     43  double axis_phi; 
     44} CoreShellCylinderParameters; 
     45 
     46 
     47/** 
     48 * Function to evaluate 2D scattering function 
     49 * @param pars: parameters of the core-shell cylinder 
     50 * @param q: q-value 
     51 * @param q_x: q_x / q 
     52 * @param q_y: q_y / q 
     53 * @return: function value 
     54 */ 
     55static double core_shell_cylinder_analytical_2D_scaled(CoreShellCylinderParameters *pars, double q, double q_x, double q_y) { 
     56  double cyl_x, cyl_y, cyl_z; 
     57  double q_z; 
     58  double alpha, vol, cos_val; 
     59  double answer; 
     60  //convert angle degree to radian 
     61  double pi = 4.0*atan(1.0); 
     62  double theta = pars->axis_theta * pi/180.0; 
     63  double phi = pars->axis_phi * pi/180.0; 
     64 
     65    // Cylinder orientation 
     66    cyl_x = sin(theta) * cos(phi); 
     67    cyl_y = sin(theta) * sin(phi); 
     68    cyl_z = cos(theta); 
     69 
     70    // q vector 
     71    q_z = 0; 
     72 
     73    // Compute the angle btw vector q and the 
     74    // axis of the cylinder 
     75    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z; 
     76 
     77    // The following test should always pass 
     78    if (fabs(cos_val)>1.0) { 
     79      printf("core_shell_cylinder_analytical_2D: Unexpected error: cos(alpha)=%g\n", cos_val); 
     80      return 0; 
     81    } 
     82 
     83  alpha = acos( cos_val ); 
     84 
     85  // Call the IGOR library function to get the kernel 
     86  answer = CoreShellCylKernel(q, pars->radius, pars->thickness, 
     87                pars->core_sld,pars->shell_sld, 
     88                pars->solvent_sld, pars->length/2.0, alpha) / fabs(sin(alpha)); 
     89 
     90  //normalize by cylinder volume 
     91  vol=pi*(pars->radius+pars->thickness) 
     92      *(pars->radius+pars->thickness) 
     93      *(pars->length+2.0*pars->thickness); 
     94  answer /= vol; 
     95 
     96  //convert to [cm-1] 
     97  answer *= 1.0e8; 
     98 
     99  //Scale 
     100  answer *= pars->scale; 
     101 
     102  // add in the background 
     103  answer += pars->background; 
     104 
     105  return answer; 
     106} 
     107 
     108/** 
     109 * Function to evaluate 2D scattering function 
     110 * @param pars: parameters of the core-shell cylinder 
     111 * @param q: q-value 
     112 * @return: function value 
     113 */ 
     114static double core_shell_cylinder_analytical_2DXY(CoreShellCylinderParameters *pars, double qx, double qy) { 
     115  double q; 
     116  q = sqrt(qx*qx+qy*qy); 
     117    return core_shell_cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q); 
     118} 
     119 
    34120 
    35121CoreShellCylinderModel :: CoreShellCylinderModel() { 

  • sansmodels/src/c_models/coreshellsphere.cpp

    r67424cd r011e0e4  
    2121 
    2222#include <math.h> 
    23 #include "models.hh" 
    2423#include "parameters.hh" 
    2524#include <stdio.h> 
    2625using namespace std; 
     26#include "core_shell.h" 
    2727 
    2828extern "C" { 
    29         #include "libSphere.h" 
    30         #include "core_shell.h" 
     29#include "libSphere.h" 
    3130} 
    3231 
     32typedef struct { 
     33  double scale; 
     34  double radius; 
     35  double thickness; 
     36  double core_sld; 
     37  double shell_sld; 
     38  double solvent_sld; 
     39  double background; 
     40} CoreShellParameters; 
     41 
    3342CoreShellModel :: CoreShellModel() { 
    34         scale      = Parameter(1.0); 
    35         radius     = Parameter(60.0, true); 
    36         radius.set_min(0.0); 
    37         thickness  = Parameter(10.0, true); 
    38         thickness.set_min(0.0); 
    39         core_sld   = Parameter(1.e-6); 
    40         shell_sld  = Parameter(2.e-6); 
    41         solvent_sld = Parameter(3.e-6); 
    42         background = Parameter(0.0); 
     43  scale      = Parameter(1.0); 
     44  radius     = Parameter(60.0, true); 
     45  radius.set_min(0.0); 
     46  thickness  = Parameter(10.0, true); 
     47  thickness.set_min(0.0); 
     48  core_sld   = Parameter(1.e-6); 
     49  shell_sld  = Parameter(2.e-6); 
     50  solvent_sld = Parameter(3.e-6); 
     51  background = Parameter(0.0); 
    4352} 
    4453 
    … …  
    5059 */ 
    5160double CoreShellModel :: operator()(double q) { 
    52         double dp[7]; 
     61  double dp[7]; 
    5362 
    54         // Fill parameter array for IGOR library 
    55         // Add the background after averaging 
     63  // Fill parameter array for IGOR library 
     64  // Add the background after averaging 
    5665 
    57         dp[0] = scale(); 
    58         dp[1] = radius(); 
    59         dp[2] = thickness(); 
    60         dp[3] = core_sld(); 
    61         dp[4] = shell_sld(); 
    62         dp[5] = solvent_sld(); 
    63         dp[6] = 0.0; 
     66  dp[0] = scale(); 
     67  dp[1] = radius(); 
     68  dp[2] = thickness(); 
     69  dp[3] = core_sld(); 
     70  dp[4] = shell_sld(); 
     71  dp[5] = solvent_sld(); 
     72  dp[6] = 0.0; 
    6473 
    6574 
    66         // Get the dispersion points for the radius 
    67         vector<WeightPoint> weights_rad; 
    68         radius.get_weights(weights_rad); 
     75  // Get the dispersion points for the radius 
     76  vector<WeightPoint> weights_rad; 
     77  radius.get_weights(weights_rad); 
    6978 
    70         // Get the dispersion points for the thickness 
    71         vector<WeightPoint> weights_thick; 
    72         thickness.get_weights(weights_thick); 
     79  // Get the dispersion points for the thickness 
     80  vector<WeightPoint> weights_thick; 
     81  thickness.get_weights(weights_thick); 
    7382 
    74         // Perform the computation, with all weight points 
    75         double sum = 0.0; 
    76         double norm = 0.0; 
    77         double vol = 0.0; 
     83  // Perform the computation, with all weight points 
     84  double sum = 0.0; 
     85  double norm = 0.0; 
     86  double vol = 0.0; 
    7887 
    79         // Loop over radius weight points 
    80         for(size_t i=0; i<weights_rad.size(); i++) { 
    81                 dp[1] = weights_rad[i].value; 
     88  // Loop over radius weight points 
     89  for(size_t i=0; i<weights_rad.size(); i++) { 
     90    dp[1] = weights_rad[i].value; 
    8291 
    83                 // Loop over thickness weight points 
    84                 for(size_t j=0; j<weights_thick.size(); j++) { 
    85                         dp[2] = weights_thick[j].value; 
    86                         //Un-normalize SphereForm by volume 
    87                         sum += weights_rad[i].weight 
    88                                 * weights_thick[j].weight * CoreShellForm(dp, q)* pow(weights_rad[i].value+weights_thick[j].value,3); 
     92    // Loop over thickness weight points 
     93    for(size_t j=0; j<weights_thick.size(); j++) { 
     94      dp[2] = weights_thick[j].value; 
     95      //Un-normalize SphereForm by volume 
     96      sum += weights_rad[i].weight 
     97          * weights_thick[j].weight * CoreShellForm(dp, q)* pow(weights_rad[i].value+weights_thick[j].value,3); 
    8998 
    90                         //Find average volume 
    91                         vol += weights_rad[i].weight * weights_thick[j].weight 
    92                                 * pow(weights_rad[i].value+weights_thick[j].value,3); 
    93                         norm += weights_rad[i].weight 
    94                                 * weights_thick[j].weight; 
    95                 } 
    96         } 
     99      //Find average volume 
     100      vol += weights_rad[i].weight * weights_thick[j].weight 
     101          * pow(weights_rad[i].value+weights_thick[j].value,3); 
     102      norm += weights_rad[i].weight 
     103          * weights_thick[j].weight; 
     104    } 
     105  } 
    97106 
    98         if (vol != 0.0 && norm != 0.0) { 
    99                 //Re-normalize by avg volume 
    100                 sum = sum/(vol/norm);} 
     107  if (vol != 0.0 && norm != 0.0) { 
     108    //Re-normalize by avg volume 
     109    sum = sum/(vol/norm);} 
    101110 
    102         return sum/norm + background(); 
     111  return sum/norm + background(); 
    103112} 
    104113 
    … …  
    110119 */ 
    111120double CoreShellModel :: operator()(double qx, double qy) { 
    112         double q = sqrt(qx*qx + qy*qy); 
    113         return (*this).operator()(q); 
     121  double q = sqrt(qx*qx + qy*qy); 
     122  return (*this).operator()(q); 
    114123} 
    115124 
    … …  
    122131 */ 
    123132double CoreShellModel :: evaluate_rphi(double q, double phi) { 
    124         return (*this).operator()(q); 
     133  return (*this).operator()(q); 
    125134} 
    126135/** 
    … …  
    129138 */ 
    130139double CoreShellModel :: calculate_ER() { 
    131         CoreShellParameters dp; 
     140  CoreShellParameters dp; 
    132141 
    133         dp.radius     = radius(); 
    134         dp.thickness  = thickness(); 
     142  dp.radius     = radius(); 
     143  dp.thickness  = thickness(); 
    135144 
    136         double rad_out = 0.0; 
     145  double rad_out = 0.0; 
    137146 
    138         // Perform the computation, with all weight points 
    139         double sum = 0.0; 
    140         double norm = 0.0; 
     147  // Perform the computation, with all weight points 
     148  double sum = 0.0; 
     149  double norm = 0.0; 
    141150 
    142151 
    143         // Get the dispersion points for the major shell 
    144         vector<WeightPoint> weights_thickness; 
    145         thickness.get_weights(weights_thickness); 
     152  // Get the dispersion points for the major shell 
     153  vector<WeightPoint> weights_thickness; 
     154  thickness.get_weights(weights_thickness); 
    146155 
    147         // Get the dispersion points for the minor shell 
    148         vector<WeightPoint> weights_radius ; 
    149         radius.get_weights(weights_radius); 
     156  // Get the dispersion points for the minor shell 
     157  vector<WeightPoint> weights_radius ; 
     158  radius.get_weights(weights_radius); 
    150159 
    151         // Loop over major shell weight points 
    152         for(int j=0; j< (int)weights_thickness.size(); j++) { 
    153                 dp.thickness = weights_thickness[j].value; 
    154                 for(int k=0; k< (int)weights_radius.size(); k++) { 
    155                         dp.radius = weights_radius[k].value; 
    156                         sum += weights_thickness[j].weight 
    157                                 * weights_radius[k].weight*(dp.radius+dp.thickness); 
    158                         norm += weights_thickness[j].weight* weights_radius[k].weight; 
    159                 } 
    160         } 
    161         if (norm != 0){ 
    162                 //return the averaged value 
    163                 rad_out =  sum/norm;} 
    164         else{ 
    165                 //return normal value 
    166                 rad_out = (dp.radius+dp.thickness);} 
     160  // Loop over major shell weight points 
     161  for(int j=0; j< (int)weights_thickness.size(); j++) { 
     162    dp.thickness = weights_thickness[j].value; 
     163    for(int k=0; k< (int)weights_radius.size(); k++) { 
     164      dp.radius = weights_radius[k].value; 
     165      sum += weights_thickness[j].weight 
     166          * weights_radius[k].weight*(dp.radius+dp.thickness); 
     167      norm += weights_thickness[j].weight* weights_radius[k].weight; 
     168    } 
     169  } 
     170  if (norm != 0){ 
     171    //return the averaged value 
     172    rad_out =  sum/norm;} 
     173  else{ 
     174    //return normal value 
     175    rad_out = (dp.radius+dp.thickness);} 
    167176 
    168         return rad_out; 
     177  return rad_out; 
    169178} 

  • sansmodels/src/c_models/ellipsoid.cpp

    r67424cd r011e0e4  
    1818 *   sansmodels/src/libigor 
    1919 * 
    20  *      TODO: refactor so that we pull in the old sansmodels.c_extensions 
    2120 */ 
    2221 
    2322#include <math.h> 
    24 #include "models.hh" 
    2523#include "parameters.hh" 
    2624#include <stdio.h> 
    2725using namespace std; 
     26#include "ellipsoid.h" 
    2827 
    2928extern "C" { 
    3029        #include "libCylinder.h" 
    3130        #include "libStructureFactor.h" 
    32         #include "ellipsoid.h" 
     31} 
     32 
     33typedef struct { 
     34  double scale; 
     35  double radius_a; 
     36  double radius_b; 
     37  double sldEll; 
     38  double sldSolv; 
     39  double background; 
     40  double axis_theta; 
     41  double axis_phi; 
     42} EllipsoidParameters; 
     43 
     44/** 
     45 * Function to evaluate 2D scattering function 
     46 * @param pars: parameters of the ellipsoid 
     47 * @param q: q-value 
     48 * @param q_x: q_x / q 
     49 * @param q_y: q_y / q 
     50 * @return: function value 
     51 */ 
     52static double ellipsoid_analytical_2D_scaled(EllipsoidParameters *pars, double q, double q_x, double q_y) { 
     53  double cyl_x, cyl_y, cyl_z; 
     54  double q_z; 
     55  double alpha, vol, cos_val; 
     56  double answer; 
     57  //convert angle degree to radian 
     58  double pi = 4.0*atan(1.0); 
     59  double theta = pars->axis_theta * pi/180.0; 
     60  double phi = pars->axis_phi * pi/180.0; 
     61 
     62    // Ellipsoid orientation 
     63    cyl_x = sin(theta) * cos(phi); 
     64    cyl_y = sin(theta) * sin(phi); 
     65    cyl_z = cos(theta); 
     66 
     67    // q vector 
     68    q_z = 0.0; 
     69 
     70    // Compute the angle btw vector q and the 
     71    // axis of the cylinder 
     72    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z; 
     73 
     74    // The following test should always pass 
     75    if (fabs(cos_val)>1.0) { 
     76      printf("ellipsoid_ana_2D: Unexpected error: cos(alpha)>1\n"); 
     77      return 0; 
     78    } 
     79 
     80    // Angle between rotation axis and q vector 
     81  alpha = acos( cos_val ); 
     82 
     83  // Call the IGOR library function to get the kernel 
     84  answer = EllipsoidKernel(q, pars->radius_b, pars->radius_a, cos_val); 
     85 
     86  // Multiply by contrast^2 
     87  answer *= (pars->sldEll - pars->sldSolv) * (pars->sldEll - pars->sldSolv); 
     88 
     89  //normalize by cylinder volume 
     90    vol = 4.0/3.0 * acos(-1.0) * pars->radius_b * pars->radius_b * pars->radius_a; 
     91  answer *= vol; 
     92 
     93  //convert to [cm-1] 
     94  answer *= 1.0e8; 
     95 
     96  //Scale 
     97  answer *= pars->scale; 
     98 
     99  // add in the background 
     100  answer += pars->background; 
     101 
     102  return answer; 
     103} 
     104 
     105/** 
     106 * Function to evaluate 2D scattering function 
     107 * @param pars: parameters of the ellipsoid 
     108 * @param q: q-value 
     109 * @return: function value 
     110 */ 
     111static double ellipsoid_analytical_2DXY(EllipsoidParameters *pars, double qx, double qy) { 
     112  double q; 
     113  q = sqrt(qx*qx+qy*qy); 
     114    return ellipsoid_analytical_2D_scaled(pars, q, qx/q, qy/q); 
    33115} 
    34116 

  • sansmodels/src/c_models/models.hh

    r0ba3b08 r011e0e4  
    2727using namespace std; 
    2828 
    29  
    30  
    31 class CoreShellModel{ 
    32 public: 
    33         // Model parameters 
    34         Parameter radius; 
    35         Parameter scale; 
    36         Parameter thickness; 
    37         Parameter core_sld; 
    38         Parameter shell_sld; 
    39         Parameter solvent_sld; 
    40         Parameter background; 
    41  
    42         // Constructor 
    43         CoreShellModel(); 
    44  
    45         // Operators to get I(Q) 
    46         double operator()(double q); 
    47         double operator()(double qx, double qy); 
    48         double calculate_ER(); 
    49         double evaluate_rphi(double q, double phi); 
    50 }; 
    51  
    52 class CoreFourShellModel{ 
    53 public: 
    54         // Model parameters 
    55         Parameter scale; 
    56         Parameter rad_core0; 
    57         Parameter sld_core0; 
    58         Parameter thick_shell1; 
    59         Parameter sld_shell1; 
    60         Parameter thick_shell2; 
    61         Parameter sld_shell2; 
    62         Parameter thick_shell3; 
    63         Parameter sld_shell3; 
    64         Parameter thick_shell4; 
    65         Parameter sld_shell4; 
    66         Parameter sld_solv; 
    67         Parameter background; 
    68  
    69         // Constructor 
    70         CoreFourShellModel(); 
    71  
    72         // Operators to get I(Q) 
    73         double operator()(double q); 
    74         double operator()(double qx, double qy); 
    75         double calculate_ER(); 
    76         double evaluate_rphi(double q, double phi); 
    77 }; 
    78  
    79 class CoreShellCylinderModel{ 
    80 public: 
    81         // Model parameters 
    82         Parameter radius; 
    83         Parameter scale; 
    84         Parameter thickness; 
    85         Parameter length; 
    86         Parameter core_sld; 
    87         Parameter shell_sld; 
    88         Parameter solvent_sld; 
    89         Parameter background; 
    90         Parameter axis_theta; 
    91         Parameter axis_phi; 
    92  
    93         // Constructor 
    94         CoreShellCylinderModel(); 
    95  
    96         // Operators to get I(Q) 
    97         double operator()(double q); 
    98         double operator()(double qx, double qy); 
    99         double calculate_ER(); 
    100         double evaluate_rphi(double q, double phi); 
    101 }; 
    102  
    103 class EllipsoidModel{ 
    104 public: 
    105         // Model parameters 
    106         Parameter radius_a; 
    107         Parameter scale; 
    108         Parameter radius_b; 
    109         Parameter sldEll; 
    110         Parameter sldSolv; 
    111         Parameter background; 
    112         Parameter axis_theta; 
    113         Parameter axis_phi; 
    114  
    115         // Constructor 
    116         EllipsoidModel(); 
    117  
    118         // Operators to get I(Q) 
    119         double operator()(double q); 
    120         double operator()(double qx, double qy); 
    121         double calculate_ER(); 
    122         double evaluate_rphi(double q, double phi); 
    123 }; 
    124  
    12529class EllipticalCylinderModel{ 
    12630public: 
Note: See TracChangeset for help on using the changeset viewer.