Changeset 0ba3b08 in sasview for sansmodels/src/c_models


Ignore:
Timestamp:
Jan 5, 2012 12:16:29 PM (12 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

Location:
sansmodels/src/c_models
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/c_models/DiamCyl.cpp

    r67424cd r0ba3b08  
    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 "DiamCyl.h" 
    2827 
    2928extern "C" { 
    3029        #include "libStructureFactor.h" 
    31         #include "DiamCyl.h" 
    3230} 
    3331 
    … …  
    112110  return 0.0; 
    113111} 
    114 // Testing code 
    115  
    116 /** 
    117 int main(void) 
    118 { 
    119         DiamCylFunc c = DiamCylFunc(); 
    120  
    121         printf("I(Qx=%g,Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001)); 
    122         printf("I(Q=%g) = %g\n", 0.001, c(0.001)); 
    123         c.radius.dispersion = new GaussianDispersion(); 
    124         c.radius.dispersion->npts = 100; 
    125         c.radius.dispersion->width = 20; 
    126  
    127         //c.length.dispersion = GaussianDispersion(); 
    128         //c.length.dispersion.npts = 20; 
    129         //c.length.dispersion.width = 65; 
    130  
    131         //printf("I(Q=%g) = %g\n", 0.001, c(0.001)); 
    132         //printf("I(Q=%g) = %g\n", 0.001, c(0.001)); 
    133         //printf("I(Qx=%g, Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001)); 
    134         //printf("I(Q=%g,  Phi=%g) = %g\n", 0.00447, .7854, c.evaluate_rphi(sqrt(0.00002), .7854)); 
    135  
    136  
    137  
    138         double i_avg = c(0.01, 0.01); 
    139         double i_1d = c(sqrt(0.0002)); 
    140  
    141         printf("\nI(Qx=%g, Qy=%g) = %g\n", 0.01, 0.01, i_avg); 
    142         printf("I(Q=%g)         = %g\n", sqrt(0.0002), i_1d); 
    143         printf("ratio %g %g\n", i_avg/i_1d, i_1d/i_avg); 
    144  
    145  
    146         return 0; 
    147 } 
    148  **/ 

  • sansmodels/src/c_models/DiamEllip.cpp

    r67424cd r0ba3b08  
    2121 
    2222#include <math.h> 
    23 #include "models.hh" 
    2423#include "parameters.hh" 
    2524#include <stdio.h> 
    2625using namespace std; 
     26#include "DiamEllip.h" 
    2727 
    2828extern "C" { 
    2929        #include "libStructureFactor.h" 
    30         #include "DiamEllip.h" 
    3130} 
    3231 
    … …  
    109108  return 0.0; 
    110109} 
    111 // Testing code 
    112 /* 
    113 int main(void) 
    114 { 
    115         SquareWellModel c = SquareWellModel(); 
    116  
    117         printf("I(Qx=%g,Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001)); 
    118         printf("I(Q=%g) = %g\n", 0.001, c(0.001)); 
    119         c.radius.dispersion = new GaussianDispersion(); 
    120         c.radius.dispersion->npts = 100; 
    121         c.radius.dispersion->width = 5; 
    122  
    123         //c.length.dispersion = GaussianDispersion(); 
    124         //c.length.dispersion.npts = 20; 
    125         //c.length.dispersion.width = 65; 
    126  
    127         printf("I(Q=%g) = %g\n", 0.001, c(0.001)); 
    128         printf("I(Q=%g) = %g\n", 0.001, c(0.001)); 
    129         printf("I(Qx=%g, Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001)); 
    130         printf("I(Q=%g,  Phi=%g) = %g\n", 0.00447, .7854, c.evaluate_rphi(sqrt(0.00002), .7854)); 
    131  
    132  
    133  
    134         double i_avg = c(0.01, 0.01); 
    135         double i_1d = c(sqrt(0.0002)); 
    136  
    137         printf("\nI(Qx=%g, Qy=%g) = %g\n", 0.01, 0.01, i_avg); 
    138         printf("I(Q=%g)         = %g\n", sqrt(0.0002), i_1d); 
    139         printf("ratio %g %g\n", i_avg/i_1d, i_1d/i_avg); 
    140  
    141  
    142         return 0; 
    143 } 
    144 */ 

  • sansmodels/src/c_models/Hardsphere.cpp

    r67424cd r0ba3b08  
    2222 
    2323#include <math.h> 
    24 #include "models.hh" 
    2524#include "parameters.hh" 
    2625#include <stdio.h> 
    2726using namespace std; 
     27#include "Hardsphere.h" 
    2828 
    2929extern "C" { 
    30         #include "libStructureFactor.h" 
    31         #include "Hardsphere.h" 
     30#include "libStructureFactor.h" 
    3231} 
    3332 
    3433HardsphereStructure :: HardsphereStructure() { 
    35         effect_radius      = Parameter(50.0, true); 
    36         effect_radius.set_min(0.0); 
    37         volfraction = Parameter(0.20, true); 
    38         volfraction.set_min(0.0); 
     34  effect_radius      = Parameter(50.0, true); 
     35  effect_radius.set_min(0.0); 
     36  volfraction = Parameter(0.20, true); 
     37  volfraction.set_min(0.0); 
    3938} 
    4039 
    … …  
    4645 */ 
    4746double HardsphereStructure :: operator()(double q) { 
    48         double dp[2]; 
     47  double dp[2]; 
    4948 
    50         // Fill parameter array for IGOR library 
    51         // Add the background after averaging 
    52         dp[0] = effect_radius(); 
    53         dp[1] = volfraction(); 
     49  // Fill parameter array for IGOR library 
     50  // Add the background after averaging 
     51  dp[0] = effect_radius(); 
     52  dp[1] = volfraction(); 
    5453 
    55         // Get the dispersion points for the radius 
    56         vector<WeightPoint> weights_rad; 
    57         effect_radius.get_weights(weights_rad); 
     54  // Get the dispersion points for the radius 
     55  vector<WeightPoint> weights_rad; 
     56  effect_radius.get_weights(weights_rad); 
    5857 
    59         // Perform the computation, with all weight points 
    60         double sum = 0.0; 
    61         double norm = 0.0; 
     58  // Perform the computation, with all weight points 
     59  double sum = 0.0; 
     60  double norm = 0.0; 
    6261 
    63         // Loop over radius weight points 
    64         for(size_t i=0; i<weights_rad.size(); i++) { 
    65                 dp[0] = weights_rad[i].value; 
     62  // Loop over radius weight points 
     63  for(size_t i=0; i<weights_rad.size(); i++) { 
     64    dp[0] = weights_rad[i].value; 
    6665 
    67                 sum += weights_rad[i].weight 
    68                                 * HardSphereStruct(dp, q); 
    69                 norm += weights_rad[i].weight; 
    70         } 
    71         return sum/norm ; 
     66    sum += weights_rad[i].weight 
     67        * HardSphereStruct(dp, q); 
     68    norm += weights_rad[i].weight; 
     69  } 
     70  return sum/norm ; 
    7271} 
    7372 
    … …  
    7978 */ 
    8079double HardsphereStructure :: operator()(double qx, double qy) { 
    81         double q = sqrt(qx*qx + qy*qy); 
    82         return (*this).operator()(q); 
     80  double q = sqrt(qx*qx + qy*qy); 
     81  return (*this).operator()(q); 
    8382} 
    8483 
    … …  
    9190 */ 
    9291double HardsphereStructure :: evaluate_rphi(double q, double phi) { 
    93         double qx = q*cos(phi); 
    94         double qy = q*sin(phi); 
    95         return (*this).operator()(qx, qy); 
     92  double qx = q*cos(phi); 
     93  double qy = q*sin(phi); 
     94  return (*this).operator()(qx, qy); 
    9695} 
    9796/** 
    … …  
    10099 */ 
    101100double HardsphereStructure :: calculate_ER() { 
    102 //NOT implemented yet!!! 
     101  //NOT implemented yet!!! 
    103102  return 0.0; 
    104103} 

  • sansmodels/src/c_models/HayterMSA.cpp

    r67424cd r0ba3b08  
    2222 
    2323#include <math.h> 
    24 #include "models.hh" 
    2524#include "parameters.hh" 
    2625#include <stdio.h> 
    2726using namespace std; 
     27#include "HayterMSA.h" 
    2828 
    2929extern "C" { 
    30         #include "libStructureFactor.h" 
    31         #include "HayterMSA.h" 
     30#include "libStructureFactor.h" 
    3231} 
    3332 
    3433HayterMSAStructure :: HayterMSAStructure() { 
    35         effect_radius      = Parameter(20.75, true); 
    36         effect_radius.set_min(0.0); 
    37         charge      = Parameter(19.0, true); 
    38         volfraction = Parameter(0.0192, true); 
    39         volfraction.set_min(0.0); 
    40         temperature = Parameter(318.16, true); 
    41         temperature.set_min(0.0); 
    42         saltconc   = Parameter(0.0); 
    43         dielectconst  = Parameter(71.08); 
     34  effect_radius      = Parameter(20.75, true); 
     35  effect_radius.set_min(0.0); 
     36  charge      = Parameter(19.0, true); 
     37  volfraction = Parameter(0.0192, true); 
     38  volfraction.set_min(0.0); 
     39  temperature = Parameter(318.16, true); 
     40  temperature.set_min(0.0); 
     41  saltconc   = Parameter(0.0); 
     42  dielectconst  = Parameter(71.08); 
    4443} 
    4544 
    … …  
    5150 */ 
    5251double HayterMSAStructure :: operator()(double q) { 
    53         double dp[6]; 
     52  double dp[6]; 
    5453 
    55         // Fill parameter array for IGOR library 
    56         // Add the background after averaging 
    57         dp[0] = 2.0*effect_radius(); 
    58         dp[1] = fabs(charge()); 
    59         dp[2] = volfraction(); 
    60         dp[3] = temperature(); 
    61         dp[4] = saltconc(); 
    62         dp[5] = dielectconst(); 
     54  // Fill parameter array for IGOR library 
     55  // Add the background after averaging 
     56  dp[0] = 2.0*effect_radius(); 
     57  dp[1] = fabs(charge()); 
     58  dp[2] = volfraction(); 
     59  dp[3] = temperature(); 
     60  dp[4] = saltconc(); 
     61  dp[5] = dielectconst(); 
    6362 
    64         // Get the dispersion points for the radius 
    65         vector<WeightPoint> weights_rad; 
    66         effect_radius.get_weights(weights_rad); 
     63  // Get the dispersion points for the radius 
     64  vector<WeightPoint> weights_rad; 
     65  effect_radius.get_weights(weights_rad); 
    6766 
    68         // Perform the computation, with all weight points 
    69         double sum = 0.0; 
    70         double norm = 0.0; 
     67  // Perform the computation, with all weight points 
     68  double sum = 0.0; 
     69  double norm = 0.0; 
    7170 
    72         // Loop over radius weight points 
    73         for(size_t i=0; i<weights_rad.size(); i++) { 
    74                 dp[0] = 2.0*weights_rad[i].value; 
     71  // Loop over radius weight points 
     72  for(size_t i=0; i<weights_rad.size(); i++) { 
     73    dp[0] = 2.0*weights_rad[i].value; 
    7574 
    76                 sum += weights_rad[i].weight 
    77                                 * HayterPenfoldMSA(dp, q); 
    78                 norm += weights_rad[i].weight; 
    79         } 
    80         return sum/norm ; 
     75    sum += weights_rad[i].weight 
     76        * HayterPenfoldMSA(dp, q); 
     77    norm += weights_rad[i].weight; 
     78  } 
     79  return sum/norm ; 
    8180} 
    8281 
    … …  
    8887 */ 
    8988double HayterMSAStructure :: operator()(double qx, double qy) { 
    90         double q = sqrt(qx*qx + qy*qy); 
    91         return (*this).operator()(q); 
     89  double q = sqrt(qx*qx + qy*qy); 
     90  return (*this).operator()(q); 
    9291} 
    9392 
    … …  
    10099 */ 
    101100double HayterMSAStructure :: evaluate_rphi(double q, double phi) { 
    102         double qx = q*cos(phi); 
    103         double qy = q*sin(phi); 
    104         return (*this).operator()(qx, qy); 
     101  double qx = q*cos(phi); 
     102  double qy = q*sin(phi); 
     103  return (*this).operator()(qx, qy); 
    105104} 
    106105/** 
    … …  
    109108 */ 
    110109double HayterMSAStructure :: calculate_ER() { 
    111 //NOT implemented yet!!! 
     110  //NOT implemented yet!!! 
    112111  return 0.0; 
    113112} 
    114  
    115 // Testing code 
    116 /* 
    117 int main(void) 
    118 { 
    119         SquareWellModel c = SquareWellModel(); 
    120  
    121         printf("I(Qx=%g,Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001)); 
    122         printf("I(Q=%g) = %g\n", 0.001, c(0.001)); 
    123         c.radius.dispersion = new GaussianDispersion(); 
    124         c.radius.dispersion->npts = 100; 
    125         c.radius.dispersion->width = 5; 
    126  
    127         //c.length.dispersion = GaussianDispersion(); 
    128         //c.length.dispersion.npts = 20; 
    129         //c.length.dispersion.width = 65; 
    130  
    131         printf("I(Q=%g) = %g\n", 0.001, c(0.001)); 
    132         printf("I(Q=%g) = %g\n", 0.001, c(0.001)); 
    133         printf("I(Qx=%g, Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001)); 
    134         printf("I(Q=%g,  Phi=%g) = %g\n", 0.00447, .7854, c.evaluate_rphi(sqrt(0.00002), .7854)); 
    135  
    136  
    137  
    138         double i_avg = c(0.01, 0.01); 
    139         double i_1d = c(sqrt(0.0002)); 
    140  
    141         printf("\nI(Qx=%g, Qy=%g) = %g\n", 0.01, 0.01, i_avg); 
    142         printf("I(Q=%g)         = %g\n", sqrt(0.0002), i_1d); 
    143         printf("ratio %g %g\n", i_avg/i_1d, i_1d/i_avg); 
    144  
    145  
    146         return 0; 
    147 } 
    148 */ 

  • sansmodels/src/c_models/SquareWell.cpp

    r67424cd r0ba3b08  
    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 "SquareWell.h" 
    2827 
    2928extern "C" { 
    30         #include "libStructureFactor.h" 
    31         #include "SquareWell.h" 
     29#include "libStructureFactor.h" 
    3230} 
    3331 
    3432SquareWellStructure :: SquareWellStructure() { 
    35         effect_radius      = Parameter(50.0, true); 
    36         effect_radius.set_min(0.0); 
    37         volfraction = Parameter(0.04, true); 
    38         volfraction.set_min(0.0); 
    39         welldepth   = Parameter(1.50); 
    40         wellwidth  = Parameter(1.20); 
     33  effect_radius      = Parameter(50.0, true); 
     34  effect_radius.set_min(0.0); 
     35  volfraction = Parameter(0.04, true); 
     36  volfraction.set_min(0.0); 
     37  welldepth   = Parameter(1.50); 
     38  wellwidth  = Parameter(1.20); 
    4139} 
    4240 
    … …  
    4846 */ 
    4947double SquareWellStructure :: operator()(double q) { 
    50         double dp[4]; 
     48  double dp[4]; 
    5149 
    52         // Fill parameter array for IGOR library 
    53         // Add the background after averaging 
    54         dp[0] = effect_radius(); 
    55         dp[1] = volfraction(); 
    56         dp[2] = welldepth(); 
    57         dp[3] = wellwidth(); 
     50  // Fill parameter array for IGOR library 
     51  // Add the background after averaging 
     52  dp[0] = effect_radius(); 
     53  dp[1] = volfraction(); 
     54  dp[2] = welldepth(); 
     55  dp[3] = wellwidth(); 
    5856 
    59         // Get the dispersion points for the radius 
    60         vector<WeightPoint> weights_rad; 
    61         effect_radius.get_weights(weights_rad); 
     57  // Get the dispersion points for the radius 
     58  vector<WeightPoint> weights_rad; 
     59  effect_radius.get_weights(weights_rad); 
    6260 
    63         // Perform the computation, with all weight points 
    64         double sum = 0.0; 
    65         double norm = 0.0; 
     61  // Perform the computation, with all weight points 
     62  double sum = 0.0; 
     63  double norm = 0.0; 
    6664 
    67         // Loop over radius weight points 
    68         for(size_t i=0; i<weights_rad.size(); i++) { 
    69                 dp[0] = weights_rad[i].value; 
     65  // Loop over radius weight points 
     66  for(size_t i=0; i<weights_rad.size(); i++) { 
     67    dp[0] = weights_rad[i].value; 
    7068 
    71                 sum += weights_rad[i].weight 
    72                                 * SquareWellStruct(dp, q); 
    73                 norm += weights_rad[i].weight; 
    74         } 
    75         return sum/norm ; 
     69    sum += weights_rad[i].weight 
     70        * SquareWellStruct(dp, q); 
     71    norm += weights_rad[i].weight; 
     72  } 
     73  return sum/norm ; 
    7674} 
    7775 
    … …  
    8381 */ 
    8482double SquareWellStructure :: operator()(double qx, double qy) { 
    85         double q = sqrt(qx*qx + qy*qy); 
    86         return (*this).operator()(q); 
     83  double q = sqrt(qx*qx + qy*qy); 
     84  return (*this).operator()(q); 
    8785} 
    8886 
    … …  
    9593 */ 
    9694double SquareWellStructure :: evaluate_rphi(double q, double phi) { 
    97         double qx = q*cos(phi); 
    98         double qy = q*sin(phi); 
    99         return (*this).operator()(qx, qy); 
     95  double qx = q*cos(phi); 
     96  double qy = q*sin(phi); 
     97  return (*this).operator()(qx, qy); 
    10098} 
    10199/** 
    … …  
    104102 */ 
    105103double SquareWellStructure :: calculate_ER() { 
    106 //NOT implemented yet!!! 
     104  //NOT implemented yet!!! 
    107105  return 0.0; 
    108106} 
    109 // Testing code 
    110 /* 
    111 int main(void) 
    112 { 
    113         SquareWellModel c = SquareWellModel(); 
    114  
    115         printf("I(Qx=%g,Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001)); 
    116         printf("I(Q=%g) = %g\n", 0.001, c(0.001)); 
    117         c.radius.dispersion = new GaussianDispersion(); 
    118         c.radius.dispersion->npts = 100; 
    119         c.radius.dispersion->width = 5; 
    120  
    121         //c.length.dispersion = GaussianDispersion(); 
    122         //c.length.dispersion.npts = 20; 
    123         //c.length.dispersion.width = 65; 
    124  
    125         printf("I(Q=%g) = %g\n", 0.001, c(0.001)); 
    126         printf("I(Q=%g) = %g\n", 0.001, c(0.001)); 
    127         printf("I(Qx=%g, Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001)); 
    128         printf("I(Q=%g,  Phi=%g) = %g\n", 0.00447, .7854, c.evaluate_rphi(sqrt(0.00002), .7854)); 
    129  
    130  
    131  
    132         double i_avg = c(0.01, 0.01); 
    133         double i_1d = c(sqrt(0.0002)); 
    134  
    135         printf("\nI(Qx=%g, Qy=%g) = %g\n", 0.01, 0.01, i_avg); 
    136         printf("I(Q=%g)         = %g\n", sqrt(0.0002), i_1d); 
    137         printf("ratio %g %g\n", i_avg/i_1d, i_1d/i_avg); 
    138  
    139  
    140         return 0; 
    141 } 
    142 */ 

  • sansmodels/src/c_models/StickyHS.cpp

    r67424cd r0ba3b08  
    2222 
    2323#include <math.h> 
    24 #include "models.hh" 
    2524#include "parameters.hh" 
    2625#include <stdio.h> 
    2726using namespace std; 
     27#include "StickyHS.h" 
    2828 
    2929extern "C" { 
    3030        #include "libStructureFactor.h" 
    31         #include "StickyHS.h" 
    3231} 
    3332 

  • sansmodels/src/c_models/bcc.cpp

    r67424cd r0ba3b08  
    2121 
    2222#include <math.h> 
    23 #include "models.hh" 
    2423#include "parameters.hh" 
    2524#include <stdio.h> 
    2625using namespace std; 
     26#include "bcc.h" 
    2727 
    2828extern "C" { 
    2929        #include "libSphere.h" 
    30         #include "bcc.h" 
     30} 
     31 
     32// Convenience structure 
     33typedef struct { 
     34  double scale; 
     35  double dnn; 
     36  double d_factor; 
     37  double radius; 
     38  double sldSph; 
     39  double sldSolv; 
     40  double background; 
     41  double theta; 
     42  double phi; 
     43  double psi; 
     44 
     45} BCParameters; 
     46 
     47/** 
     48 * Function to evaluate 2D scattering function 
     49 * @param pars: parameters of the BCCCrystalModel 
     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 bc_analytical_2D_scaled(BCParameters *pars, double q, double q_x, double q_y) { 
     56  double b3_x, b3_y, b3_z, b1_x, b1_y; 
     57  double q_z; 
     58  double alpha, cos_val_b3, cos_val_b2, cos_val_b1; 
     59  double a1_dot_q, a2_dot_q,a3_dot_q; 
     60  double answer; 
     61  double Pi = 4.0*atan(1.0); 
     62  double aa, Da, qDa_2, latticeScale, Zq, Fkq, Fkq_2; 
     63 
     64  //convert angle degree to radian 
     65  double pi = 4.0*atan(1.0); 
     66  double theta = pars->theta * pi/180.0; 
     67  double phi = pars->phi * pi/180.0; 
     68  double psi = pars->psi * pi/180.0; 
     69 
     70  double dp[5]; 
     71  dp[0] = 1.0; 
     72  dp[1] = pars->radius; 
     73  dp[2] = pars->sldSph; 
     74  dp[3] = pars->sldSolv; 
     75  dp[4] = 0.0; 
     76 
     77  aa = pars->dnn; 
     78  Da = pars->d_factor*aa; 
     79  qDa_2 = pow(q*Da,2.0); 
     80 
     81  //the occupied volume of the lattice 
     82  latticeScale = 2.0*(4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa/sqrt(3.0/4.0),3.0); 
     83  // q vector 
     84  q_z = 0.0; // for SANS; assuming qz is negligible 
     85  /// Angles here are respect to detector coordinate 
     86  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854)) 
     87    // b3 axis orientation 
     88    b3_x = sin(theta) * cos(phi);//negative sign here??? 
     89    b3_y = sin(theta) * sin(phi); 
     90    b3_z = cos(theta); 
     91    cos_val_b3 =  b3_x*q_x + b3_y*q_y + b3_z*q_z; 
     92 
     93    alpha = acos(cos_val_b3); 
     94    // b1 axis orientation 
     95    b1_x = sin(psi); 
     96    b1_y = cos(psi); 
     97  cos_val_b1 = (b1_x*q_x + b1_y*q_y); 
     98    // b2 axis orientation 
     99  cos_val_b2 = sin(acos(cos_val_b1)); 
     100  // alpha corrections 
     101  cos_val_b2 *= sin(alpha); 
     102  cos_val_b1 *= sin(alpha); 
     103 
     104    // Compute the angle btw vector q and the a3 axis 
     105    a3_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b1-cos_val_b3); 
     106 
     107    // a1 axis 
     108    a1_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b2-cos_val_b1); 
     109 
     110    // a2 axis 
     111    a2_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b1-cos_val_b2); 
     112 
     113    // The following test should always pass 
     114    if (fabs(cos_val_b3)>1.0) { 
     115      printf("bcc_ana_2D: Unexpected error: cos(alpha)>1\n"); 
     116      return 0; 
     117    } 
     118    // Get Fkq and Fkq_2 
     119    Fkq = exp(-0.5*pow(Da/aa,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q)); 
     120    Fkq_2 = Fkq*Fkq; 
     121    // Call Zq=Z1*Z2*Z3 
     122    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2); 
     123    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2); 
     124    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2); 
     125 
     126  // Use SphereForm directly from libigor 
     127  answer = SphereForm(dp,q)*Zq; 
     128 
     129  //consider scales 
     130  answer *= latticeScale * pars->scale; 
     131 
     132  // This FIXES a singualrity the kernel in libigor. 
     133  if ( answer == INFINITY || answer == NAN){ 
     134    answer = 0.0; 
     135  } 
     136 
     137  // add background 
     138  answer += pars->background; 
     139 
     140  return answer; 
     141} 
     142 
     143/** 
     144 * Function to evaluate 2D scattering function 
     145 * @param pars: parameters of the BCC_ParaCrystal 
     146 * @param q: q-value 
     147 * @return: function value 
     148 */ 
     149static double bc_analytical_2DXY(BCParameters *pars, double qx, double qy){ 
     150  double q; 
     151  q = sqrt(qx*qx+qy*qy); 
     152  return bc_analytical_2D_scaled(pars, q, qx/q, qy/q); 
    31153} 
    32154 

  • sansmodels/src/c_models/fcc.cpp

    r67424cd r0ba3b08  
    2121 
    2222#include <math.h> 
    23 #include "models.hh" 
    2423#include "parameters.hh" 
    2524#include <stdio.h> 
    2625using namespace std; 
     26#include "fcc.h" 
    2727 
    2828extern "C" { 
    29         #include "libSphere.h" 
    30         #include "fcc.h" 
     29#include "libSphere.h" 
     30} 
     31 
     32// Convenience parameter structure 
     33typedef struct { 
     34  double scale; 
     35  double dnn; 
     36  double d_factor; 
     37  double radius; 
     38  double sldSph; 
     39  double sldSolv; 
     40  double background; 
     41  double theta; 
     42  double phi; 
     43  double psi; 
     44} FCParameters; 
     45 
     46/** 
     47 * Function to evaluate 2D scattering function 
     48 * @param pars: parameters of the FCCCrystalModel 
     49 * @param q: q-value 
     50 * @param q_x: q_x / q 
     51 * @param q_y: q_y / q 
     52 * @return: function value 
     53 */ 
     54static double fc_analytical_2D_scaled(FCParameters *pars, double q, double q_x, double q_y) { 
     55  double b3_x, b3_y, b3_z, b1_x, b1_y; 
     56  double q_z; 
     57  double alpha, cos_val_b3, cos_val_b2, cos_val_b1; 
     58  double a1_dot_q, a2_dot_q,a3_dot_q; 
     59  double answer; 
     60  double Pi = 4.0*atan(1.0); 
     61  double aa, Da, qDa_2, latticeScale, Zq, Fkq, Fkq_2; 
     62 
     63  double dp[5]; 
     64  //convert angle degree to radian 
     65  double theta = pars->theta * Pi/180.0; 
     66  double phi = pars->phi * Pi/180.0; 
     67  double psi = pars->psi * Pi/180.0; 
     68  dp[0] = 1.0; 
     69  dp[1] = pars->radius; 
     70  dp[2] = pars->sldSph; 
     71  dp[3] = pars->sldSolv; 
     72  dp[4] = 0.0; 
     73 
     74  aa = pars->dnn; 
     75  Da = pars->d_factor*aa; 
     76  qDa_2 = pow(q*Da,2.0); 
     77  //contrast = pars->sldSph - pars->sldSolv; 
     78 
     79  latticeScale = 4.0*(4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa*sqrt(2.0),3.0); 
     80  // q vector 
     81  q_z = 0.0; // for SANS; assuming qz is negligible 
     82  /// Angles here are respect to detector coordinate 
     83  ///  instead of against q coordinate in PRB 36(46), 3(6), 1754(3854) 
     84  // b3 axis orientation 
     85  b3_x = sin(theta) * cos(phi);//negative sign here??? 
     86  b3_y = sin(theta) * sin(phi); 
     87  b3_z = cos(theta); 
     88  cos_val_b3 =  b3_x*q_x + b3_y*q_y + b3_z*q_z; 
     89  alpha = acos(cos_val_b3); 
     90  // b1 axis orientation 
     91  b1_x = sin(psi); 
     92  b1_y = cos(psi); 
     93  cos_val_b1 = (b1_x*q_x + b1_y*q_y); 
     94  // b2 axis orientation 
     95  cos_val_b2 = sin(acos(cos_val_b1)); 
     96  // alpha correction 
     97  cos_val_b2 *= sin(alpha); 
     98  cos_val_b1 *= sin(alpha); 
     99 
     100  // Compute the angle btw vector q and the a3 axis 
     101  a3_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b1); 
     102 
     103  // a1 axis 
     104  a1_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b3); 
     105 
     106  // a2 axis 
     107  a2_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b1); 
     108 
     109  // The following test should always pass 
     110  if (fabs(cos_val_b3)>1.0) { 
     111    printf("fcc_ana_2D: Unexpected error: cos(alpha)>1\n"); 
     112    return 0; 
     113  } 
     114  // Get Fkq and Fkq_2 
     115  Fkq = exp(-0.5*pow(Da/aa,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q)); 
     116  Fkq_2 = Fkq*Fkq; 
     117  // Call Zq=Z1*Z2*Z3 
     118  Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2); 
     119  Zq = Zq * (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2); 
     120  Zq = Zq * (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2); 
     121 
     122  // Use SphereForm directly from libigor 
     123  answer = SphereForm(dp,q)*Zq; 
     124 
     125  //consider scales 
     126  answer *= latticeScale * pars->scale; 
     127 
     128  // This FIXES a singualrity the kernel in libigor. 
     129  if ( answer == INFINITY || answer == NAN){ 
     130    answer = 0.0; 
     131  } 
     132 
     133  // add background 
     134  answer += pars->background; 
     135 
     136  return answer; 
     137} 
     138 
     139/** 
     140 * Function to evaluate 2D scattering function 
     141 * @param pars: parameters of the FCC_ParaCrystal 
     142 * @param q: q-value 
     143 * @return: function value 
     144 */ 
     145static double fc_analytical_2DXY(FCParameters *pars, double qx, double qy){ 
     146  double q; 
     147  q = sqrt(qx*qx+qy*qy); 
     148  return fc_analytical_2D_scaled(pars, q, qx/q, qy/q); 
    31149} 
    32150 
    33151FCCrystalModel :: FCCrystalModel() { 
    34         scale      = Parameter(1.0); 
    35         dnn             = Parameter(220.0); 
    36         d_factor = Parameter(0.06); 
    37         radius     = Parameter(40.0, true); 
    38         radius.set_min(0.0); 
    39         sldSph   = Parameter(3.0e-6); 
    40         sldSolv   = Parameter(6.3e-6); 
    41         background = Parameter(0.0); 
    42         theta  = Parameter(0.0, true); 
    43         phi    = Parameter(0.0, true); 
    44         psi    = Parameter(0.0, true); 
     152  scale      = Parameter(1.0); 
     153  dnn           = Parameter(220.0); 
     154  d_factor = Parameter(0.06); 
     155  radius     = Parameter(40.0, true); 
     156  radius.set_min(0.0); 
     157  sldSph   = Parameter(3.0e-6); 
     158  sldSolv   = Parameter(6.3e-6); 
     159  background = Parameter(0.0); 
     160  theta  = Parameter(0.0, true); 
     161  phi    = Parameter(0.0, true); 
     162  psi    = Parameter(0.0, true); 
    45163} 
    46164 
    … …  
    52170 */ 
    53171double FCCrystalModel :: operator()(double q) { 
    54         double dp[7]; 
    55  
    56         // Fill parameter array for IGOR library 
    57         // Add the background after averaging 
    58         dp[0] = scale(); 
    59         dp[1] = dnn(); 
    60         dp[2] = d_factor(); 
    61         dp[3] = radius(); 
    62         dp[4] = sldSph(); 
    63         dp[5] = sldSolv(); 
    64         dp[6] = 0.0; 
    65  
    66         // Get the dispersion points for the radius 
    67         vector<WeightPoint> weights_rad; 
    68         radius.get_weights(weights_rad); 
    69  
    70         // Perform the computation, with all weight points 
    71         double sum = 0.0; 
    72         double norm = 0.0; 
    73         double vol = 0.0; 
    74         double result; 
    75  
    76         // Loop over radius weight points 
    77         for(size_t i=0; i<weights_rad.size(); i++) { 
    78                 dp[3] = weights_rad[i].value; 
    79  
    80                 //Un-normalize SphereForm by volume 
    81                 result = FCC_ParaCrystal(dp, q) * pow(weights_rad[i].value,3); 
    82                 // This FIXES a singualrity the kernel in libigor. 
    83                 if ( result == INFINITY || result == NAN){ 
    84                         result = 0.0; 
    85                 } 
    86                 sum += weights_rad[i].weight 
    87                         * result; 
    88                 //Find average volume 
    89                 vol += weights_rad[i].weight 
    90                         * pow(weights_rad[i].value,3); 
    91  
    92                 norm += weights_rad[i].weight; 
    93         } 
    94  
    95         if (vol != 0.0 && norm != 0.0) { 
    96                 //Re-normalize by avg volume 
    97                 sum = sum/(vol/norm);} 
    98         return sum/norm + background(); 
     172  double dp[7]; 
     173 
     174  // Fill parameter array for IGOR library 
     175  // Add the background after averaging 
     176  dp[0] = scale(); 
     177  dp[1] = dnn(); 
     178  dp[2] = d_factor(); 
     179  dp[3] = radius(); 
     180  dp[4] = sldSph(); 
     181  dp[5] = sldSolv(); 
     182  dp[6] = 0.0; 
     183 
     184  // Get the dispersion points for the radius 
     185  vector<WeightPoint> weights_rad; 
     186  radius.get_weights(weights_rad); 
     187 
     188  // Perform the computation, with all weight points 
     189  double sum = 0.0; 
     190  double norm = 0.0; 
     191  double vol = 0.0; 
     192  double result; 
     193 
     194  // Loop over radius weight points 
     195  for(size_t i=0; i<weights_rad.size(); i++) { 
     196    dp[3] = weights_rad[i].value; 
     197 
     198    //Un-normalize SphereForm by volume 
     199    result = FCC_ParaCrystal(dp, q) * pow(weights_rad[i].value,3); 
     200    // This FIXES a singualrity the kernel in libigor. 
     201    if ( result == INFINITY || result == NAN){ 
     202      result = 0.0; 
     203    } 
     204    sum += weights_rad[i].weight 
     205        * result; 
     206    //Find average volume 
     207    vol += weights_rad[i].weight 
     208        * pow(weights_rad[i].value,3); 
     209 
     210    norm += weights_rad[i].weight; 
     211  } 
     212 
     213  if (vol != 0.0 && norm != 0.0) { 
     214    //Re-normalize by avg volume 
     215    sum = sum/(vol/norm);} 
     216  return sum/norm + background(); 
    99217} 
    100218 
    … …  
    106224 */ 
    107225double FCCrystalModel :: operator()(double qx, double qy) { 
    108         FCParameters dp; 
    109         dp.scale      = scale(); 
    110         dp.dnn   = dnn(); 
    111         dp.d_factor   = d_factor(); 
    112         dp.radius  = radius(); 
    113         dp.sldSph   = sldSph(); 
    114         dp.sldSolv   = sldSolv(); 
    115         dp.background = 0.0; 
    116         dp.theta  = theta(); 
    117         dp.phi    = phi(); 
    118         dp.psi    = psi(); 
    119  
    120         // Get the dispersion points for the radius 
    121         vector<WeightPoint> weights_rad; 
    122         radius.get_weights(weights_rad); 
    123  
    124         // Get angular averaging for theta 
    125         vector<WeightPoint> weights_theta; 
    126         theta.get_weights(weights_theta); 
    127  
    128         // Get angular averaging for phi 
    129         vector<WeightPoint> weights_phi; 
    130         phi.get_weights(weights_phi); 
    131  
    132         // Get angular averaging for psi 
    133         vector<WeightPoint> weights_psi; 
    134         psi.get_weights(weights_psi); 
    135  
    136         // Perform the computation, with all weight points 
    137         double sum = 0.0; 
    138         double norm = 0.0; 
    139         double norm_vol = 0.0; 
    140         double vol = 0.0; 
    141         double pi = 4.0*atan(1.0); 
    142         // Loop over radius weight points 
    143         for(size_t i=0; i<weights_rad.size(); i++) { 
    144                 dp.radius = weights_rad[i].value; 
    145                 // Average over theta distribution 
    146                 for(size_t j=0; j< weights_theta.size(); j++) { 
    147                         dp.theta = weights_theta[j].value; 
    148                         // Average over phi distribution 
    149                         for(size_t k=0; k< weights_phi.size(); k++) { 
    150                                 dp.phi = weights_phi[k].value; 
    151                                 // Average over phi distribution 
    152                                 for(size_t l=0; l< weights_psi.size(); l++) { 
    153                                         dp.psi = weights_psi[l].value; 
    154                                         //Un-normalize SphereForm by volume 
    155                                         double _ptvalue = weights_rad[i].weight 
    156                                                                                 * weights_theta[j].weight 
    157                                                                                 * weights_phi[k].weight 
    158                                                                                 * weights_psi[l].weight 
    159                                                                                 * fc_analytical_2DXY(&dp, qx, qy); 
    160                                                                                 //* pow(weights_rad[i].value,3.0); 
    161                                         // Consider when there is infinte or nan. 
    162                                         if ( _ptvalue == INFINITY || _ptvalue == NAN){ 
    163                                                 _ptvalue = 0.0; 
    164                                         } 
    165                                         if (weights_theta.size()>1) { 
    166                                                 _ptvalue *= fabs(sin(weights_theta[j].value*pi/180.0)); 
    167                                         } 
    168                                         sum += _ptvalue; 
    169                                         // This model dose not need the volume of spheres correction!!! 
    170                                         norm += weights_rad[i].weight 
    171                                                         * weights_theta[j].weight 
    172                                                         * weights_phi[k].weight 
    173                                                         * weights_psi[l].weight; 
    174                                 } 
    175                         } 
    176                 } 
    177         } 
    178         // Averaging in theta needs an extra normalization 
    179         // factor to account for the sin(theta) term in the 
    180         // integration (see documentation). 
    181         if (weights_theta.size()>1) norm = norm / asin(1.0); 
    182  
    183         if (vol != 0.0 && norm_vol != 0.0) { 
    184                 //Re-normalize by avg volume 
    185                 sum = sum/(vol/norm_vol);} 
    186  
    187         return sum/norm + background(); 
     226  FCParameters dp; 
     227  dp.scale      = scale(); 
     228  dp.dnn   = dnn(); 
     229  dp.d_factor   = d_factor(); 
     230  dp.radius  = radius(); 
     231  dp.sldSph   = sldSph(); 
     232  dp.sldSolv   = sldSolv(); 
     233  dp.background = 0.0; 
     234  dp.theta  = theta(); 
     235  dp.phi    = phi(); 
     236  dp.psi    = psi(); 
     237 
     238  // Get the dispersion points for the radius 
     239  vector<WeightPoint> weights_rad; 
     240  radius.get_weights(weights_rad); 
     241 
     242  // Get angular averaging for theta 
     243  vector<WeightPoint> weights_theta; 
     244  theta.get_weights(weights_theta); 
     245 
     246  // Get angular averaging for phi 
     247  vector<WeightPoint> weights_phi; 
     248  phi.get_weights(weights_phi); 
     249 
     250  // Get angular averaging for psi 
     251  vector<WeightPoint> weights_psi; 
     252  psi.get_weights(weights_psi); 
     253 
     254  // Perform the computation, with all weight points 
     255  double sum = 0.0; 
     256  double norm = 0.0; 
     257  double norm_vol = 0.0; 
     258  double vol = 0.0; 
     259  double pi = 4.0*atan(1.0); 
     260  // Loop over radius weight points 
     261  for(size_t i=0; i<weights_rad.size(); i++) { 
     262    dp.radius = weights_rad[i].value; 
     263    // Average over theta distribution 
     264    for(size_t j=0; j< weights_theta.size(); j++) { 
     265      dp.theta = weights_theta[j].value; 
     266      // Average over phi distribution 
     267      for(size_t k=0; k< weights_phi.size(); k++) { 
     268        dp.phi = weights_phi[k].value; 
     269        // Average over phi distribution 
     270        for(size_t l=0; l< weights_psi.size(); l++) { 
     271          dp.psi = weights_psi[l].value; 
     272          //Un-normalize SphereForm by volume 
     273          double _ptvalue = weights_rad[i].weight 
     274              * weights_theta[j].weight 
     275              * weights_phi[k].weight 
     276              * weights_psi[l].weight 
     277              * fc_analytical_2DXY(&dp, qx, qy); 
     278          //* pow(weights_rad[i].value,3.0); 
     279          // Consider when there is infinte or nan. 
     280          if ( _ptvalue == INFINITY || _ptvalue == NAN){ 
     281            _ptvalue = 0.0; 
     282          } 
     283          if (weights_theta.size()>1) { 
     284            _ptvalue *= fabs(sin(weights_theta[j].value*pi/180.0)); 
     285          } 
     286          sum += _ptvalue; 
     287          // This model dose not need the volume of spheres correction!!! 
     288          norm += weights_rad[i].weight 
     289              * weights_theta[j].weight 
     290              * weights_phi[k].weight 
     291              * weights_psi[l].weight; 
     292        } 
     293      } 
     294    } 
     295  } 
     296  // Averaging in theta needs an extra normalization 
     297  // factor to account for the sin(theta) term in the 
     298  // integration (see documentation). 
     299  if (weights_theta.size()>1) norm = norm / asin(1.0); 
     300 
     301  if (vol != 0.0 && norm_vol != 0.0) { 
     302    //Re-normalize by avg volume 
     303    sum = sum/(vol/norm_vol);} 
     304 
     305  return sum/norm + background(); 
    188306} 
    189307 
    … …  
    196314 */ 
    197315double FCCrystalModel :: evaluate_rphi(double q, double phi) { 
    198         return (*this).operator()(q); 
     316  return (*this).operator()(q); 
    199317} 
    200318 
    … …  
    204322 */ 
    205323double FCCrystalModel :: calculate_ER() { 
    206         //NOT implemented yet!!! 
    207         return 0.0; 
    208 } 
     324  //NOT implemented yet!!! 
     325  return 0.0; 
     326} 

  • sansmodels/src/c_models/fuzzysphere.cpp

    r67424cd r0ba3b08  
    1515 
    1616#include <math.h> 
    17 #include "models.hh" 
    1817#include "parameters.hh" 
    1918#include <stdio.h> 
     19#include <stdlib.h> 
    2020using namespace std; 
     21#include "fuzzysphere.h" 
    2122 
    2223extern "C" { 
    23         #include "fuzzysphere.h" 
     24#include "libSphere.h" 
     25} 
     26 
     27// scattering from a uniform sphere w/ fuzzy surface 
     28// Modified from FuzzySpheres in libigor/libSphere.c without polydispersion: JHC 
     29static double fuzzysphere_kernel(double dp[], double q){ 
     30  double pi,x,xr; 
     31  double radius,sldSph,sldSolv,scale,bkg,delrho,fuzziness,f2,bes,vol,f;   //my local names 
     32 
     33  pi = 4.0*atan(1.0); 
     34  x= q; 
     35  scale = dp[0]; 
     36  radius = dp[1]; 
     37  fuzziness = dp[2]; 
     38  sldSph = dp[3]; 
     39  sldSolv = dp[4]; 
     40  bkg = dp[5]; 
     41  delrho=sldSph-sldSolv; 
     42 
     43  xr = x*radius; 
     44  //handle xr==0 separately 
     45  if(xr == 0.0){ 
     46    bes = 1.0; 
     47  }else{ 
     48    bes = 3.0*(sin(xr)-xr*cos(xr))/(xr*xr*xr); 
     49  } 
     50  vol = 4.0*pi/3.0*radius*radius*radius; 
     51  f = vol*bes*delrho;   // [=] A 
     52  f *= exp(-0.5*fuzziness*fuzziness*x*x); 
     53  // normalize to single particle volume, convert to 1/cm 
     54  f2 = f * f / vol * 1.0e8;   // [=] 1/cm 
     55 
     56  f2 *= scale; 
     57  f2 += bkg; 
     58 
     59    return(f2); //scale, and add in the background 
    2460} 
    2561 

  • sansmodels/src/c_models/models.hh

    r6e10cff r0ba3b08  
    2727using namespace std; 
    2828 
    29  
    30  
    31 class PearlNecklaceModel{ 
    32 public: 
    33         // Model parameters 
    34         Parameter scale; 
    35         Parameter radius; 
    36         Parameter edge_separation; 
    37         Parameter thick_string; 
    38         Parameter num_pearls; 
    39         Parameter sld_pearl; 
    40         Parameter sld_string; 
    41         Parameter sld_solv; 
    42         Parameter background; 
    43  
    44         // Constructor 
    45         PearlNecklaceModel(); 
    46  
    47         // Operators to get I(Q) 
    48         double operator()(double q); 
    49         double operator()(double qx, double qy); 
    50         double calculate_ER(); 
    51         double evaluate_rphi(double q, double phi); 
    52 }; 
    53  
    54 class FCCrystalModel{ 
    55 public: 
    56         // Model parameters 
    57         Parameter scale; 
    58         Parameter dnn; 
    59         Parameter d_factor; 
    60         Parameter radius; 
    61         Parameter sldSph; 
    62         Parameter sldSolv; 
    63         Parameter background; 
    64         Parameter theta; 
    65         Parameter phi; 
    66         Parameter psi; 
    67  
    68         // Constructor 
    69         FCCrystalModel(); 
    70  
    71         // Operators to get I(Q) 
    72         double operator()(double q); 
    73         double operator()(double qx, double qy); 
    74         double calculate_ER(); 
    75         double evaluate_rphi(double q, double phi); 
    76 }; 
    77  
    78  
    79 class BCCrystalModel{ 
    80 public: 
    81         // Model parameters 
    82         Parameter scale; 
    83         Parameter dnn; 
    84         Parameter d_factor; 
    85         Parameter radius; 
    86         Parameter sldSph; 
    87         Parameter sldSolv; 
    88         Parameter background; 
    89         Parameter theta; 
    90         Parameter phi; 
    91         Parameter psi; 
    92  
    93         // Constructor 
    94         BCCrystalModel(); 
    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  
    104 class FuzzySphereModel{ 
    105 public: 
    106         // Model parameters 
    107         Parameter radius; 
    108         Parameter scale; 
    109         Parameter fuzziness; 
    110         Parameter sldSph; 
    111         Parameter sldSolv; 
    112         Parameter background; 
    113  
    114         // Constructor 
    115         FuzzySphereModel(); 
    116  
    117         // Operators to get I(Q) 
    118         double operator()(double q); 
    119         double operator()(double qx, double qy); 
    120         double calculate_ER(); 
    121         double evaluate_rphi(double q, double phi); 
    122 }; 
    123  
    124 class HardsphereStructure{ 
    125 public: 
    126         // Model parameters 
    127         Parameter effect_radius; 
    128         Parameter volfraction; 
    129  
    130         // Constructor 
    131         HardsphereStructure(); 
    132  
    133         // Operators to get I(Q) 
    134         double operator()(double q); 
    135         double operator()(double qx, double qy); 
    136         double calculate_ER(); 
    137         double evaluate_rphi(double q, double phi); 
    138 }; 
    139  
    140 class StickyHSStructure{ 
    141 public: 
    142         // Model parameters 
    143         Parameter effect_radius; 
    144         Parameter volfraction; 
    145         Parameter perturb; 
    146         Parameter stickiness; 
    147  
    148         // Constructor 
    149         StickyHSStructure(); 
    150  
    151         // Operators to get I(Q) 
    152         double operator()(double q); 
    153         double operator()(double qx, double qy); 
    154         double calculate_ER(); 
    155         double evaluate_rphi(double q, double phi); 
    156 }; 
    157  
    158 class SquareWellStructure{ 
    159 public: 
    160         // Model parameters 
    161         Parameter effect_radius; 
    162         Parameter volfraction; 
    163         Parameter welldepth; 
    164         Parameter wellwidth; 
    165  
    166         // Constructor 
    167         SquareWellStructure(); 
    168  
    169         // Operators to get I(Q) 
    170         double operator()(double q); 
    171         double operator()(double qx, double qy); 
    172         double calculate_ER(); 
    173         double evaluate_rphi(double q, double phi); 
    174 }; 
    175  
    176 class HayterMSAStructure{ 
    177 public: 
    178         // Model parameters 
    179         Parameter effect_radius; 
    180         Parameter charge; 
    181         Parameter volfraction; 
    182         Parameter temperature; 
    183         Parameter saltconc; 
    184         Parameter dielectconst; 
    185  
    186         // Constructor 
    187         HayterMSAStructure(); 
    188  
    189         // Operators to get I(Q) 
    190         double operator()(double q); 
    191         double operator()(double qx, double qy); 
    192         double calculate_ER(); 
    193         double evaluate_rphi(double q, double phi); 
    194 }; 
    195  
    196 class DiamEllipFunc{ 
    197 public: 
    198         // Model parameters 
    199         Parameter radius_a; 
    200         Parameter radius_b; 
    201  
    202         // Constructor 
    203         DiamEllipFunc(); 
    204  
    205         // Operators to get I(Q) 
    206         double operator()(double q); 
    207         double operator()(double qx, double qy); 
    208         double calculate_ER(); 
    209         double evaluate_rphi(double q, double phi); 
    210 }; 
    211  
    212 class DiamCylFunc{ 
    213 public: 
    214         // Model parameters 
    215         Parameter radius; 
    216         Parameter length; 
    217  
    218         // Constructor 
    219         DiamCylFunc(); 
    220  
    221         // Operators to get I(Q) 
    222         double operator()(double q); 
    223         double operator()(double qx, double qy); 
    224         double calculate_ER(); 
    225         double evaluate_rphi(double q, double phi); 
    226 }; 
    227  
    228  
    229 class SLDCalFunc{ 
    230 public: 
    231         // Model parameters 
    232         Parameter fun_type; 
    233         Parameter npts_inter; 
    234         Parameter shell_num; 
    235         Parameter nu_inter; 
    236         Parameter sld_left; 
    237         Parameter sld_right; 
    238  
    239         // Constructor 
    240         SLDCalFunc(); 
    241  
    242         // Operators to get SLD 
    243         double operator()(double q); 
    244         double operator()(double qx, double qy); 
    245         double calculate_ER(); 
    246         double evaluate_rphi(double q, double phi); 
    247 }; 
    24829 
    24930 

  • 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} 

  • sansmodels/src/c_models/sld_cal.cpp

    r67424cd r0ba3b08  
    66 
    77#include <math.h> 
    8 #include "models.hh" 
    98#include "parameters.hh" 
    109#include <stdio.h> 
    1110using namespace std; 
     11#include "sld_cal.h" 
    1212 
    1313extern "C" { 
    14         #include "sld_cal.h" 
     14#include "libmultifunc/librefl.h" 
     15} 
     16 
     17// Convenience structure 
     18typedef struct { 
     19    double fun_type; 
     20    double npts_inter; 
     21    double shell_num; 
     22    double nu_inter; 
     23    double sld_left; 
     24    double sld_right; 
     25} SLDCalParameters; 
     26 
     27/** 
     28 * Function to calculate sld 
     29 * @param pars: parameters 
     30 * @param x: independent param-value 
     31 * @return: sld value 
     32 */ 
     33static double sld_cal_analytical_1D(SLDCalParameters *pars, double x) { 
     34  double fun, nsl, n_s, fun_coef, sld_l, sld_r, sld_out; 
     35  int fun_type; 
     36 
     37  fun = pars->fun_type; 
     38  nsl = pars->npts_inter; 
     39  n_s = pars->shell_num; 
     40  fun_coef = pars->nu_inter; 
     41  sld_l = pars-> sld_left; 
     42  sld_r = pars-> sld_right; 
     43 
     44  fun_type = floor(fun); 
     45 
     46  sld_out = intersldfunc(fun_type, nsl, n_s, fun_coef, sld_l, sld_r); 
     47 
     48  return sld_out; 
    1549} 
    1650 
Note: See TracChangeset for help on using the changeset viewer.