Changeset 0ba3b08 in sasview for sansmodels


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

Location:
sansmodels
Files:
5 deleted
24 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/setup.py

    r67424cd r0ba3b08  
    5555print "Installing SANS models" 
    5656 
     57sys.path.append(wrapper_dir) 
     58from wrapping import generate_wrappers 
     59generate_wrappers(header_dir=srcdir,  
     60                  output_dir=os.path.join("src", "sans", "models"),  
     61                  c_wrapper_dir=wrapper_dir) 
    5762 
    5863IGNORED_FILES = ["a.exe", 

  • sansmodels/src/c_extensions/DiamCyl.h

    r67424cd r0ba3b08  
    11#if !defined(DiamCyl_h) 
    22#define DiamCyl_h 
     3#include "parameters.hh" 
    34 
    45/** 
    … …  
    1617**/ 
    1718 
    18 typedef struct { 
    19     /// Radius [A] 
    20     //  [DEFAULT]=radius=20.0 A 
    21     double radius; 
     19class DiamCylFunc{ 
     20public: 
     21  // Model parameters 
     22  /// Radius [A] 
     23  //  [DEFAULT]=radius=20.0 A 
     24  Parameter radius; 
     25  /// Length [A] 
     26  //  [DEFAULT]=length= 400 A 
     27  Parameter length; 
    2228 
    23     /// Length [A] 
    24     //  [DEFAULT]=length= 400 A 
    25     double length; 
     29  // Constructor 
     30  DiamCylFunc(); 
    2631 
    27 } DiamCyldParameters; 
    28  
    29  
    30  
    31 /// 1D scattering function 
    32 //double DiamCyld_analytical_1D(DiamCyldParameters *pars, double q); 
    33  
    34 /// 2D scattering function 
    35 //double DiamCyld_analytical_2D(DiamCyldParameters *pars, double q, double phi); 
    36 //double DiamCyld_analytical_2DXY(DiamCyldParameters *pars, double qx, double qy); 
     32  // Operators to get I(Q) 
     33  double operator()(double q); 
     34  double operator()(double qx, double qy); 
     35  double calculate_ER(); 
     36  double evaluate_rphi(double q, double phi); 
     37}; 
    3738 
    3839#endif 

  • sansmodels/src/c_extensions/DiamEllip.h

    r67424cd r0ba3b08  
    11#if !defined(DiamEllip_h) 
    22#define DiamEllip_h 
     3#include "parameters.hh" 
    34 
    45/** 
    … …  
    2021 
    2122 **/ 
    22 typedef struct { 
    23     /// Polar radius [A] 
    24     //  [DEFAULT]=radius_a=20.0 A 
    25     double radius_a; 
    2623 
    27     /// Equatorial radius [A] 
    28     //  [DEFAULT]=radius_b= 400 A 
    29     double radius_b; 
     24class DiamEllipFunc{ 
     25public: 
     26  // Model parameters 
     27  /// Polar radius [A] 
     28  //  [DEFAULT]=radius_a=20.0 A 
     29  Parameter radius_a; 
    3030 
    31 } DiamEllipsParameters; 
     31  /// Equatorial radius [A] 
     32  //  [DEFAULT]=radius_b= 400 A 
     33  Parameter radius_b; 
    3234 
     35  // Constructor 
     36  DiamEllipFunc(); 
    3337 
    34  
    35 /// 1D scattering function 
    36 //double DiamEllips_analytical_1D(DiamEllipsParameters *pars, double q); 
    37  
    38 /// 2D scattering function 
    39 //double DiamEllips_analytical_2D(DiamEllipsParameters *pars, double q, double phi); 
    40 //double DiamEllips_analytical_2DXY(DiamEllipsParameters *pars, double qx, double qy); 
     38  // Operators to get I(Q) 
     39  double operator()(double q); 
     40  double operator()(double qx, double qy); 
     41  double calculate_ER(); 
     42  double evaluate_rphi(double q, double phi); 
     43}; 
    4144 
    4245#endif 

  • sansmodels/src/c_extensions/Hardsphere.h

    r67424cd r0ba3b08  
    11#if !defined(Hardsphere_h) 
    22#define Hardsphere_h 
     3#include "parameters.hh" 
    34 
    45/** 
    56 * Structure definition for HardsphereStructure (factor) parameters 
    67 */ 
    7  //[PYTHONCLASS] = HardsphereStructure 
    8  //[DISP_PARAMS] = effect_radius 
    9  //[DESCRIPTION] =<text>Structure factor for interacting particles:                   . 
    10  // 
    11  //  The interparticle potential is 
    12  // 
    13  //                     U(r)= inf   , r < 2R 
    14  //                         = 0     , r >= 2R 
    15  // 
    16  //                                             R: effective radius of the Hardsphere particle 
    17  //                                             V:The volume fraction 
    18  // 
    19  //    Ref: Percus., J. K.,etc., J. Phy. 
    20  //     Rev. 1958, 110, 1. 
    21  //      </text> 
    22  //[FIXED]= effect_radius.width 
     8//[PYTHONCLASS] = HardsphereStructure 
     9//[DISP_PARAMS] = effect_radius 
     10//[DESCRIPTION] =<text>Structure factor for interacting particles:                   . 
     11// 
     12//  The interparticle potential is 
     13// 
     14//                     U(r)= inf   , r < 2R 
     15//                         = 0     , r >= 2R 
     16// 
     17//                                              R: effective radius of the Hardsphere particle 
     18//                                              V:The volume fraction 
     19// 
     20//    Ref: Percus., J. K.,etc., J. Phy. 
     21//     Rev. 1958, 110, 1. 
     22//       </text> 
     23//[FIXED]= effect_radius.width 
    2324 
     25class HardsphereStructure{ 
     26public: 
     27  // Model parameters 
     28  /// effect radius of hardsphere [A] 
     29  //  [DEFAULT]=effect_radius=50.0 [A] 
     30  Parameter effect_radius; 
    2431 
    25 typedef struct { 
    26     /// effect radius of hardsphere [A] 
    27     //  [DEFAULT]=effect_radius=50.0 [A] 
    28     double effect_radius; 
     32  /// Volume fraction 
     33  //  [DEFAULT]=volfraction= 0.2 
     34  Parameter volfraction; 
    2935 
    30     /// Volume fraction 
    31     //  [DEFAULT]=volfraction= 0.2 
    32     double volfraction; 
    33 } HardsphereParameters; 
     36  // Constructor 
     37  HardsphereStructure(); 
    3438 
    35  
    36  
    37 /// 1D scattering function 
    38 //double Hardsphere_analytical_1D(HardsphereParameters *pars, double q); 
    39  
    40 /// 2D scattering function 
    41 //double Hardsphere_analytical_2D(HardsphereParameters *pars, double q, double phi); 
    42 //double Hardsphere_analytical_2DXY(HardsphereParameters *pars, double qx, double qy); 
     39  // Operators to get I(Q) 
     40  double operator()(double q); 
     41  double operator()(double qx, double qy); 
     42  double calculate_ER(); 
     43  double evaluate_rphi(double q, double phi); 
     44}; 
    4345 
    4446#endif 

  • sansmodels/src/c_extensions/HayterMSA.h

    r67424cd r0ba3b08  
    11#if !defined(HayterMSA_h) 
    22#define HayterMSA_h 
     3#include "parameters.hh" 
    34 
    45/** 
    56 * Structure definition for screened Coulomb interaction 
    67 */ 
    7  //[PYTHONCLASS] = HayterMSAStructure 
    8  //[DISP_PARAMS] = effect_radius 
    9  //[DESCRIPTION] =<text>To calculate the structure factor (the Fourier transform of the 
    10  //                     pair correlation function g(r)) for a system of 
    11  //                     charged, spheroidal objects in a dielectric 
    12  //                                             medium. 
    13  //                     When combined with an appropriate form 
    14  //                     factor, this allows for inclusion of 
    15  //                     the interparticle interference effects 
    16  //                     due to screened coulomb repulsion between 
    17  //                     charged particles. 
    18  //                                             (Note: charge > 0 required.) 
    19  // 
    20  //                     Ref: JP Hansen and JB Hayter, Molecular 
    21  //                           Physics 46, 651-656 (1982). 
    22  // 
    23  //                             </text> 
    24  //[FIXED]= effect_radius.width 
     8//[PYTHONCLASS] = HayterMSAStructure 
     9//[DISP_PARAMS] = effect_radius 
     10//[DESCRIPTION] =<text>To calculate the structure factor (the Fourier transform of the 
     11//                     pair correlation function g(r)) for a system of 
     12//                     charged, spheroidal objects in a dielectric 
     13//                                              medium. 
     14//                     When combined with an appropriate form 
     15//                     factor, this allows for inclusion of 
     16//                     the interparticle interference effects 
     17//                     due to screened coulomb repulsion between 
     18//                     charged particles. 
     19//                                              (Note: charge > 0 required.) 
     20// 
     21//                     Ref: JP Hansen and JB Hayter, Molecular 
     22//                           Physics 46, 651-656 (1982). 
     23// 
     24//                              </text> 
     25//[FIXED]= effect_radius.width 
    2526 
    26 typedef struct { 
    27     /// effetitve radius of particle [A] 
    28     //  [DEFAULT]=effect_radius=20.75 [A] 
    29     double effect_radius; 
     27class HayterMSAStructure{ 
     28public: 
     29  // Model parameters 
     30  /// effetitve radius of particle [A] 
     31  //  [DEFAULT]=effect_radius=20.75 [A] 
     32  Parameter effect_radius; 
    3033 
    31     /// charge 
    32     //  [DEFAULT]=charge= 19 
    33     double charge; 
     34  /// charge 
     35  //  [DEFAULT]=charge= 19 
     36  Parameter charge; 
    3437 
    35     /// Volume fraction 
    36     //  [DEFAULT]=volfraction= 0.0192 
    37     double volfraction; 
     38  /// Volume fraction 
     39  //  [DEFAULT]=volfraction= 0.0192 
     40  Parameter volfraction; 
    3841 
    39         ///     Temperature [K] 
    40     //  [DEFAULT]=temperature= 318.16 [K] 
    41     double temperature; 
     42  /// Temperature [K] 
     43  //  [DEFAULT]=temperature= 318.16 [K] 
     44  Parameter temperature; 
    4245 
    43         ///     Monovalent salt concentration [M] 
    44     //  [DEFAULT]=saltconc= 0 [M] 
    45     double saltconc; 
     46  /// Monovalent salt concentration [M] 
     47  //  [DEFAULT]=saltconc= 0 [M] 
     48  Parameter saltconc; 
    4649 
    47     /// Dielectric constant of solvent 
    48     //  [DEFAULT]=dielectconst= 71.08 
    49     double dielectconst; 
    50 } HayterMSAParameters; 
     50  /// Dielectric constant of solvent 
     51  //  [DEFAULT]=dielectconst= 71.08 
     52  Parameter dielectconst; 
    5153 
    52 /// 1D scattering function 
    53 //double HayterMSA_analytical_1D(HayterMSAParameters *pars, double q); 
     54  // Constructor 
     55  HayterMSAStructure(); 
    5456 
    55 /// 2D scattering function 
    56 //double HayterMSA_analytical_2D(HayterMSAParameters *pars, double q, double phi); 
    57 //double HayterMSA_analytical_2DXY(HayterMSAParameters *pars, double qx, double qy); 
     57  // Operators to get I(Q) 
     58  double operator()(double q); 
     59  double operator()(double qx, double qy); 
     60  double calculate_ER(); 
     61  double evaluate_rphi(double q, double phi); 
     62}; 
    5863 
    5964#endif 

  • sansmodels/src/c_extensions/SquareWell.h

    r67424cd r0ba3b08  
    11#if !defined(SquareWell_h) 
    22#define SquareWell_h 
     3#include "parameters.hh" 
    34 
    45/**Structure definition for SquareWell parameters 
    5 */ 
     6 */ 
    67//   [PYTHONCLASS] = SquareWellStructure 
    78//   [DISP_PARAMS] = effect_radius 
    … …  
    2627//[ORIENTATION_PARAMS]= <text> </text> 
    2728 
     29class SquareWellStructure{ 
     30public: 
     31  // Model parameters 
     32  /// effective radius of particle [A] 
     33  //  [DEFAULT]=effect_radius=50.0 [A] 
     34  Parameter effect_radius; 
    2835 
    29 typedef struct { 
    30     /// effective radius of particle [A] 
    31     //  [DEFAULT]=effect_radius=50.0 [A] 
    32     double effect_radius; 
     36  /// Volume fraction 
     37  //  [DEFAULT]=volfraction= 0.040 
     38  Parameter volfraction; 
    3339 
    34     /// Volume fraction 
    35     //  [DEFAULT]=volfraction= 0.040 
    36     double volfraction; 
     40  /// Well depth [kT] 
     41  //  [DEFAULT]=welldepth= 1.50 [kT] 
     42  Parameter welldepth; 
    3743 
    38     /// Well depth [kT] 
    39     //  [DEFAULT]=welldepth= 1.50 [kT] 
    40     double welldepth; 
     44  /// Well width 
     45  //  [DEFAULT]=wellwidth= 1.20 
     46  Parameter wellwidth; 
    4147 
    42     /// Well width 
    43     //  [DEFAULT]=wellwidth= 1.20 
    44     double wellwidth; 
     48  // Constructor 
     49  SquareWellStructure(); 
    4550 
    46 } SquareWellParameters; 
    47  
    48  
    49  
    50 /// 1D scattering function 
    51 //double SquareWell_analytical_1D(SquareWellParameters *pars, double q); 
    52  
    53 /// 2D scattering function 
    54 //double SquareWell_analytical_2D(SquareWellParameters *pars, double q, double phi); 
    55 //double SquareWell_analytical_2DXY(SquareWellParameters *pars, double qx, double qy); 
     51  // Operators to get I(Q) 
     52  double operator()(double q); 
     53  double operator()(double qx, double qy); 
     54  double calculate_ER(); 
     55  double evaluate_rphi(double q, double phi); 
     56}; 
    5657 
    5758#endif 

  • sansmodels/src/c_extensions/StickyHS.h

    r67424cd r0ba3b08  
    11#if !defined(StickyHS_h) 
    22#define StickyHS_h 
     3#include "parameters.hh" 
    34 
    45/** 
    56 * Structure definition for StickyHS_struct (struncture factor) parameters 
    67 */ 
    7  //[PYTHONCLASS] = StickyHSStructure 
    8  //[DISP_PARAMS] = effect_radius 
    9  //[DESCRIPTION] =<text> Structure Factor for interacting particles:                               . 
    10  // 
    11  //  The interaction potential is 
    12  // 
    13  //                     U(r)= inf , r < 2R 
    14  //                         = -Uo  , 2R < r < 2R + w 
    15  //                         = 0   , r >= 2R +w 
    16  // 
    17  //                                             R: effective radius of the hardsphere 
    18  //                     stickiness = [exp(Uo/kT)]/(12*perturb) 
    19  //                     perturb = w/(w+ 2R) , 0.01 =< w <= 0.1 
    20  //                     w: The width of the square well ,w > 0 
    21  //                                             v: The volume fraction , v > 0 
    22  // 
    23  //                     Ref: Menon, S. V. G.,et.al., J. Chem. 
    24  //                      Phys., 1991, 95(12), 9186-9190. 
    25  //                             </text> 
    26  //[FIXED]= effect_radius.width 
    27 typedef struct { 
    28     /// effect radius of hardsphere [A] 
    29     //  [DEFAULT]=effect_radius=50.0 [A] 
    30     double effect_radius; 
     8//[PYTHONCLASS] = StickyHSStructure 
     9//[DISP_PARAMS] = effect_radius 
     10//[DESCRIPTION] =<text> Structure Factor for interacting particles:                               . 
     11// 
     12//  The interaction potential is 
     13// 
     14//                     U(r)= inf , r < 2R 
     15//                         = -Uo  , 2R < r < 2R + w 
     16//                         = 0   , r >= 2R +w 
     17// 
     18//                                              R: effective radius of the hardsphere 
     19//                     stickiness = [exp(Uo/kT)]/(12*perturb) 
     20//                     perturb = w/(w+ 2R) , 0.01 =< w <= 0.1 
     21//                     w: The width of the square well ,w > 0 
     22//                                              v: The volume fraction , v > 0 
     23// 
     24//                     Ref: Menon, S. V. G.,et.al., J. Chem. 
     25//                      Phys., 1991, 95(12), 9186-9190. 
     26//                              </text> 
     27//[FIXED]= effect_radius.width 
    3128 
    32     /// Volume fraction 
    33     //  [DEFAULT]=volfraction= 0.1 
    34     double volfraction; 
     29class StickyHSStructure{ 
     30public: 
     31  // Model parameters 
     32  /// effect radius of hardsphere [A] 
     33  //  [DEFAULT]=effect_radius=50.0 [A] 
     34  Parameter effect_radius; 
    3535 
    36         /// Perturbation Parameter ( between 0.01 - 0.1) 
    37     //  [DEFAULT]=perturb=0.05 
    38     double perturb; 
     36  /// Volume fraction 
     37  //  [DEFAULT]=volfraction= 0.1 
     38  Parameter volfraction; 
    3939 
    40     /// Stickiness 
    41     //  [DEFAULT]=stickiness= 0.2 
    42     double stickiness; 
    43 } StickyHSParameters; 
     40  /// Perturbation Parameter ( between 0.01 - 0.1) 
     41  //  [DEFAULT]=perturb=0.05 
     42  Parameter perturb; 
    4443 
     44  /// Stickiness 
     45  //  [DEFAULT]=stickiness= 0.2 
     46  Parameter stickiness; 
    4547 
     48  // Constructor 
     49  StickyHSStructure(); 
    4650 
    47 /// 1D scattering function 
    48 //double StickyHS_analytical_1D(StickyHSParameters *pars, double q); 
    49  
    50 /// 2D scattering function 
    51 //double StickyHS_analytical_2D(StickyHSParameters *pars, double q, double phi); 
    52 //double StickyHS_analytical_2DXY(StickyHSParameters *pars, double qx, double qy); 
     51  // Operators to get I(Q) 
     52  double operator()(double q); 
     53  double operator()(double qx, double qy); 
     54  double calculate_ER(); 
     55  double evaluate_rphi(double q, double phi); 
     56}; 
    5357 
    5458#endif 

  • sansmodels/src/c_extensions/bcc.h

    r67424cd r0ba3b08  
    11#if !defined(bcc_h) 
    22#define bcc_h 
     3#include "parameters.hh" 
    34 
    45/** 
    56 * Structure definition for BCC_ParaCrystal parameters 
    67 */ 
    7  //[PYTHONCLASS] = BCCrystalModel 
    8  //[DISP_PARAMS] = radius,phi, psi, theta 
    9  //[DESCRIPTION] =<text>P(q)=(scale/Vp)*V_lattice*P(q)*Z(q)+bkg where scale is the volume 
    10  //                                      fraction of sphere, 
    11  //                             Vp = volume of the primary particle, 
    12  //                             V_lattice = volume correction for 
    13  //                                     for the crystal structure, 
    14  //                             P(q)= form factor of the sphere (normalized), 
    15  //                             Z(q)= paracrystalline structure factor 
    16  //                                     for a face centered cubic structure. 
    17  //                             [Body Centered Cubic ParaCrystal Model] 
    18  //                             Parameters; 
    19  //                             scale: volume fraction of spheres 
    20  //                             bkg:background, R: radius of sphere 
    21  //                             dnn: Nearest neighbor distance 
    22  //                             d_factor: Paracrystal distortion factor 
    23  //                             radius: radius of the spheres 
    24  //                             sldSph: SLD of the sphere 
    25  //                             sldSolv: SLD of the solvent 
    26  // 
    27  //             </text> 
    28  //[FIXED]=  radius.width;phi.width;psi.width; theta.width 
    29  //[ORIENTATION_PARAMS]= <text> phi;psi; theta; phi.width;psi.width; theta.width</text> 
     8//[PYTHONCLASS] = BCCrystalModel 
     9//[DISP_PARAMS] = radius,phi, psi, theta 
     10//[DESCRIPTION] =<text>P(q)=(scale/Vp)*V_lattice*P(q)*Z(q)+bkg where scale is the volume 
     11//                                       fraction of sphere, 
     12//                              Vp = volume of the primary particle, 
     13//                              V_lattice = volume correction for 
     14//                                      for the crystal structure, 
     15//                              P(q)= form factor of the sphere (normalized), 
     16//                              Z(q)= paracrystalline structure factor 
     17//                                      for a face centered cubic structure. 
     18//                              [Body Centered Cubic ParaCrystal Model] 
     19//                              Parameters; 
     20//                              scale: volume fraction of spheres 
     21//                              bkg:background, R: radius of sphere 
     22//                              dnn: Nearest neighbor distance 
     23//                              d_factor: Paracrystal distortion factor 
     24//                              radius: radius of the spheres 
     25//                              sldSph: SLD of the sphere 
     26//                              sldSolv: SLD of the solvent 
     27// 
     28//              </text> 
     29//[FIXED]=  radius.width;phi.width;psi.width; theta.width 
     30//[ORIENTATION_PARAMS]= <text> phi;psi; theta; phi.width;psi.width; theta.width</text> 
    3031 
    31 typedef struct { 
    32     /// Scale factor 
    33     //  [DEFAULT]=scale= 1.0 
    34     double scale; 
     32class BCCrystalModel{ 
     33public: 
     34  // Model parameters 
     35  /// Scale factor 
     36  //  [DEFAULT]=scale= 1.0 
     37  Parameter scale; 
    3538 
    36     /// Nearest neighbor distance [A] 
    37     //  [DEFAULT]=dnn=220.0 [A] 
    38     double dnn; 
     39  /// Nearest neighbor distance [A] 
     40  //  [DEFAULT]=dnn=220.0 [A] 
     41  Parameter dnn; 
    3942 
    40     /// Paracrystal distortion factor g 
    41     //  [DEFAULT]=d_factor=0.06 
    42     double d_factor; 
     43  /// Paracrystal distortion factor g 
     44  //  [DEFAULT]=d_factor=0.06 
     45  Parameter d_factor; 
    4346 
    44     /// Radius of sphere [A] 
    45     //  [DEFAULT]=radius=40.0 [A] 
    46     double radius; 
     47  /// Radius of sphere [A] 
     48  //  [DEFAULT]=radius=40.0 [A] 
     49  Parameter radius; 
    4750 
    48     /// sldSph [1/A^(2)] 
    49     //  [DEFAULT]=sldSph= 3.0e-6 [1/A^(2)] 
    50     double sldSph; 
     51  /// sldSph [1/A^(2)] 
     52  //  [DEFAULT]=sldSph= 3.0e-6 [1/A^(2)] 
     53  Parameter sldSph; 
    5154 
    52     /// sldSolv [1/A^(2)] 
    53     //  [DEFAULT]=sldSolv= 6.3e-6 [1/A^(2)] 
    54     double sldSolv; 
     55  /// sldSolv [1/A^(2)] 
     56  //  [DEFAULT]=sldSolv= 6.3e-6 [1/A^(2)] 
     57  Parameter sldSolv; 
    5558 
    56         /// Incoherent Background [1/cm] 
    57         //  [DEFAULT]=background=0 [1/cm] 
    58         double background; 
    59     /// Orientation of the a1 axis w/respect incoming beam [deg] 
    60     //  [DEFAULT]=theta=0.0 [deg] 
    61     double theta; 
    62     /// Orientation of the a2 in the plane of the detector [deg] 
    63     //  [DEFAULT]=phi=0.0 [deg] 
    64     double phi; 
    65     /// Orientation of the a3 in the plane of the detector [deg] 
    66     //  [DEFAULT]=psi=0.0 [deg] 
    67     double psi; 
     59  /// Incoherent Background [1/cm] 
     60  //  [DEFAULT]=background=0 [1/cm] 
     61  Parameter background; 
     62  /// Orientation of the a1 axis w/respect incoming beam [deg] 
     63  //  [DEFAULT]=theta=0.0 [deg] 
     64  Parameter theta; 
     65  /// Orientation of the a2 in the plane of the detector [deg] 
     66  //  [DEFAULT]=phi=0.0 [deg] 
     67  Parameter phi; 
     68  /// Orientation of the a3 in the plane of the detector [deg] 
     69  //  [DEFAULT]=psi=0.0 [deg] 
     70  Parameter psi; 
    6871 
    69 } BCParameters; 
     72  // Constructor 
     73  BCCrystalModel(); 
     74 
     75  // Operators to get I(Q) 
     76  double operator()(double q); 
     77  double operator()(double qx, double qy); 
     78  double calculate_ER(); 
     79  double evaluate_rphi(double q, double phi); 
     80}; 
    7081 
    7182 
    72  
    73 /// 1D scattering function 
    74 double bc_analytical_1D(BCParameters *pars, double q); 
    75  
    76 /// 2D scattering function 
    77 double bc_analytical_2D(BCParameters *pars, double q, double phi); 
    78 double bc_analytical_2DXY(BCParameters *pars, double qx, double qy); 
    79 double bc_analytical_2D_scaled(BCParameters *pars, double q, double q_x, double q_y); 
    80  
    8183#endif 

  • sansmodels/src/c_extensions/fcc.h

    r67424cd r0ba3b08  
    11#if !defined(fcc_h) 
    22#define fcc_h 
     3#include "parameters.hh" 
    34 
    45/** 
    … …  
    2930 //[ORIENTATION_PARAMS]= <text> phi;psi; theta; phi.width;psi.width; theta.width</text> 
    3031 
    31 typedef struct { 
    32     /// Scale factor 
    33     //  [DEFAULT]=scale= 1.0 
    34     double scale; 
     32class FCCrystalModel{ 
     33public: 
     34  // Model parameters 
     35  /// Scale factor 
     36  //  [DEFAULT]=scale= 1.0 
     37  Parameter scale; 
    3538 
    36     /// Nearest neighbor distance [A] 
    37     //  [DEFAULT]=dnn=220.0 [A] 
    38     double dnn; 
     39  /// Nearest neighbor distance [A] 
     40  //  [DEFAULT]=dnn=220.0 [A] 
     41  Parameter dnn; 
    3942 
    40     /// Paracrystal distortion factor g 
    41     //  [DEFAULT]=d_factor=0.06 
    42     double d_factor; 
     43  /// Paracrystal distortion factor g 
     44  //  [DEFAULT]=d_factor=0.06 
     45  Parameter d_factor; 
    4346 
    44     /// Radius of sphere [A] 
    45     //  [DEFAULT]=radius=40.0 [A] 
    46     double radius; 
     47  /// Radius of sphere [A] 
     48  //  [DEFAULT]=radius=40.0 [A] 
     49  Parameter radius; 
    4750 
    48     /// sldSph [1/A^(2)] 
    49     //  [DEFAULT]=sldSph= 3.0e-6 [1/A^(2)] 
    50     double sldSph; 
     51  /// sldSph [1/A^(2)] 
     52  //  [DEFAULT]=sldSph= 3.0e-6 [1/A^(2)] 
     53  Parameter sldSph; 
    5154 
    52     /// sldSolv [1/A^(2)] 
    53     //  [DEFAULT]=sldSolv= 6.3e-6 [1/A^(2)] 
    54     double sldSolv; 
     55  /// sldSolv [1/A^(2)] 
     56  //  [DEFAULT]=sldSolv= 6.3e-6 [1/A^(2)] 
     57  Parameter sldSolv; 
    5558 
    56         /// Incoherent Background [1/cm] 
    57         //  [DEFAULT]=background=0 [1/cm] 
    58         double background; 
    59     /// Orientation of the a1 axis w/respect incoming beam [deg] 
    60     //  [DEFAULT]=theta=0.0 [deg] 
    61     double theta; 
    62     /// Orientation of the a2 in the plane of the detector [deg] 
    63     //  [DEFAULT]=phi=0.0 [deg] 
    64     double phi; 
    65     /// Orientation of the a3 in the plane of the detector [deg] 
    66     //  [DEFAULT]=psi=0.0 [deg] 
    67     double psi; 
     59/// Incoherent Background [1/cm] 
     60//  [DEFAULT]=background=0 [1/cm] 
     61Parameter background; 
     62  /// Orientation of the a1 axis w/respect incoming beam [deg] 
     63  //  [DEFAULT]=theta=0.0 [deg] 
     64  Parameter theta; 
     65  /// Orientation of the a2 in the plane of the detector [deg] 
     66  //  [DEFAULT]=phi=0.0 [deg] 
     67  Parameter phi; 
     68  /// Orientation of the a3 in the plane of the detector [deg] 
     69  //  [DEFAULT]=psi=0.0 [deg] 
     70  Parameter psi; 
    6871 
    69 } FCParameters; 
     72  // Constructor 
     73  FCCrystalModel(); 
    7074 
    71  
    72  
    73 /// 1D scattering function 
    74 double fc_analytical_1D(FCParameters *pars, double q); 
    75  
    76 /// 2D scattering function 
    77 double fc_analytical_2D(FCParameters *pars, double q, double phi); 
    78 double fc_analytical_2DXY(FCParameters *pars, double qx, double qy); 
    79 double fc_analytical_2D_scaled(FCParameters *pars, double q, double q_x, double q_y); 
     75  // Operators to get I(Q) 
     76  double operator()(double q); 
     77  double operator()(double qx, double qy); 
     78  double calculate_ER(); 
     79  double evaluate_rphi(double q, double phi); 
     80}; 
    8081 
    8182#endif 

  • sansmodels/src/c_extensions/fuzzysphere.h

    r67424cd r0ba3b08  
    11#if !defined(fuzzysphere_h) 
    22#define fuzzysphere_h 
     3#include "parameters.hh" 
    34 
    45/** 
    56 * Structure definition for FuzzySphere parameters 
    67 */ 
    7  //[PYTHONCLASS] = FuzzySphereModel 
    8  //[DISP_PARAMS] = radius, fuzziness 
    9  //[DESCRIPTION] =<text> 
    10  //                             scale: scale factor times volume fraction, 
    11  //                                     or just volume fraction for absolute scale data 
    12  //                             radius: radius of the solid sphere 
    13  //                             fuzziness = the STD of the height of fuzzy interfacial 
    14  //                              thickness (ie., so-called interfacial roughness) 
    15  //                             sldSph: the SLD of the sphere 
    16  //                             sldSolv: the SLD of the solvent 
    17  //                             background: incoherent background 
    18  //                     Note: By definition, this function works only when fuzziness << radius. 
    19  //             </text> 
    20  //[FIXED]=  radius.width; fuzziness.width 
    21  //[ORIENTATION_PARAMS]= <text> </text> 
     8//[PYTHONCLASS] = FuzzySphereModel 
     9//[DISP_PARAMS] = radius, fuzziness 
     10//[DESCRIPTION] =<text> 
     11//                              scale: scale factor times volume fraction, 
     12//                                      or just volume fraction for absolute scale data 
     13//                              radius: radius of the solid sphere 
     14//                              fuzziness = the STD of the height of fuzzy interfacial 
     15//                               thickness (ie., so-called interfacial roughness) 
     16//                              sldSph: the SLD of the sphere 
     17//                              sldSolv: the SLD of the solvent 
     18//                              background: incoherent background 
     19//                      Note: By definition, this function works only when fuzziness << radius. 
     20//              </text> 
     21//[FIXED]=  radius.width; fuzziness.width 
     22//[ORIENTATION_PARAMS]= <text> </text> 
    2223 
    23 typedef struct { 
    24     /// Scale factor 
    25     //  [DEFAULT]=scale= 0.01 
    26     double scale; 
     24class FuzzySphereModel{ 
     25public: 
     26  // Model parameters 
     27  /// Radius of sphere [A] 
     28  //  [DEFAULT]=radius=60.0 [A] 
     29  Parameter radius; 
     30  /// Scale factor 
     31  //  [DEFAULT]=scale= 0.01 
     32  Parameter scale; 
    2733 
    28     /// Radius of sphere [A] 
    29     //  [DEFAULT]=radius=60.0 [A] 
    30     double radius; 
    3134 
    32     /// surface roughness [A] 
    33     //  [DEFAULT]=fuzziness= 10.0 [A] 
    34     double fuzziness; 
     35  /// surface roughness [A] 
     36  //  [DEFAULT]=fuzziness= 10.0 [A] 
     37  Parameter fuzziness; 
    3538 
    36     /// sldSph [1/A^(2)] 
    37     //  [DEFAULT]=sldSph= 1.0e-6 [1/A^(2)] 
    38     double sldSph; 
     39  /// sldSph [1/A^(2)] 
     40  //  [DEFAULT]=sldSph= 1.0e-6 [1/A^(2)] 
     41  Parameter sldSph; 
    3942 
    40     /// sldSolv [1/A^(2)] 
    41     //  [DEFAULT]=sldSolv= 3.0e-6 [1/A^(2)] 
    42     double sldSolv; 
     43  /// sldSolv [1/A^(2)] 
     44  //  [DEFAULT]=sldSolv= 3.0e-6 [1/A^(2)] 
     45  Parameter sldSolv; 
    4346 
    44         /// Incoherent Background [1/cm] 
    45         //  [DEFAULT]=background=0.001 [1/cm] 
    46         double background; 
    47 } FuzzySphereParameters; 
     47  /// Incoherent Background [1/cm] 
     48  //  [DEFAULT]=background=0.001 [1/cm] 
     49  Parameter background; 
    4850 
    49 /// kernel 
    50 double fuzzysphere_kernel(double dp[], double q); 
     51  // Constructor 
     52  FuzzySphereModel(); 
    5153 
    52 /// 1D scattering function 
    53 double fuzzysphere_analytical_1D(FuzzySphereParameters *pars, double q); 
     54  // Operators to get I(Q) 
     55  double operator()(double q); 
     56  double operator()(double qx, double qy); 
     57  double calculate_ER(); 
     58  double evaluate_rphi(double q, double phi); 
     59}; 
    5460 
    55 /// 2D scattering function 
    56 double fuzzysphere_analytical_2D(FuzzySphereParameters *pars, double q, double phi); 
    57 double fuzzysphere_analytical_2DXY(FuzzySphereParameters *pars, double qx, double qy); 
    58 double fuzzysphere_analytical_2D_scaled(FuzzySphereParameters *pars, double q, double q_x, double q_y); 
    5961#endif 

  • sansmodels/src/c_extensions/pearlnecklace.h

    r67424cd r0ba3b08  
    11#if !defined(o_h) 
    22#define pearl_necklace_h 
     3#include "parameters.hh" 
    34 
    45/** 
    56 * Structure definition for sphere parameters 
    67 */ 
    7  //[PYTHONCLASS] = PearlNecklaceModel 
    8  //[DISP_PARAMS] = radius, edge_separation 
    9  //[DESCRIPTION] =<text>Calculate form factor for Pearl Necklace Model 
     8//[PYTHONCLASS] = PearlNecklaceModel 
     9//[DISP_PARAMS] = radius, edge_separation 
     10//[DESCRIPTION] =<text>Calculate form factor for Pearl Necklace Model 
    1011//                              [Macromol. Symp. 2004, 211, 25-42] 
    11  //                             Parameters: 
    12  //                             background:background 
    13  //                             scale: scale factor 
    14  //                             sld_pearl: the SLD of the pearl spheres 
    15  //                 sld_string: the SLD of the strings 
    16  //                             sld_solv: the SLD of the solvent 
    17  //                             num_pearls: number of the pearls 
    18  //                             radius: the radius of a pearl 
    19  //                             edge_separation: the length of string segment; surface to surface 
    20  //                             thick_string: thickness (ie, diameter) of the string 
    21  //             </text> 
    22  //[FIXED]=  <text>radius.width; edge_separation.width</text> 
    23  //[NON_FITTABLE_PARAMS]= <text></text> 
    24  //[ORIENTATION_PARAMS]= <text> </text> 
     12//                              Parameters: 
     13//                              background:background 
     14//                              scale: scale factor 
     15//                              sld_pearl: the SLD of the pearl spheres 
     16//                  sld_string: the SLD of the strings 
     17//                              sld_solv: the SLD of the solvent 
     18//                              num_pearls: number of the pearls 
     19//                              radius: the radius of a pearl 
     20//                              edge_separation: the length of string segment; surface to surface 
     21//                              thick_string: thickness (ie, diameter) of the string 
     22//              </text> 
     23//[FIXED]=  <text>radius.width; edge_separation.width</text> 
     24//[NON_FITTABLE_PARAMS]= <text></text> 
     25//[ORIENTATION_PARAMS]= <text> </text> 
    2526 
    2627typedef struct { 
    27     /// Scale factor 
    28     //  [DEFAULT]=scale= 1.0 
    29         double scale; 
    30     /// radius [A] 
    31     //  [DEFAULT]=radius=80.0 [A] 
    32         double radius; 
    33         ///     edge_separation 
    34         //  [DEFAULT]=edge_separation= 350 [A] 
    35         double edge_separation; 
    36         ///     thick_string [A] 
    37         //  [DEFAULT]=thick_string= 2.5 [A] 
    38         double thick_string; 
    39         ///     num_pearls 
    40         //  [DEFAULT]=num_pearls= 3 
    41         double num_pearls; 
    42         ///     sld_pearl 
    43         //  [DEFAULT]=sld_pearl= 1.0e-06 [1/A^(2)] 
    44         double sld_pearl; 
    45         ///     sld_string 
    46         //  [DEFAULT]=sld_string= 1.0e-06 [1/A^(2)] 
    47         double sld_string; 
    48         ///     sld_solv 
    49         //  [DEFAULT]=sld_solv= 6.3e-06 [1/A^(2)] 
    50         double sld_solv; 
    51         /// Background 
    52         //  [DEFAULT]=background=0 
    53         double background; 
     28  /// Scale factor 
     29  //  [DEFAULT]=scale= 1.0 
     30  double scale; 
     31  ///   radius [A] 
     32  //  [DEFAULT]=radius=80.0 [A] 
     33  double radius; 
     34  ///   edge_separation 
     35  //  [DEFAULT]=edge_separation= 350 [A] 
     36  double edge_separation; 
     37  ///   thick_string [A] 
     38  //  [DEFAULT]=thick_string= 2.5 [A] 
     39  double thick_string; 
     40  ///   num_pearls 
     41  //  [DEFAULT]=num_pearls= 3 
     42  double num_pearls; 
     43  ///   sld_pearl 
     44  //  [DEFAULT]=sld_pearl= 1.0e-06 [1/A^(2)] 
     45  double sld_pearl; 
     46  ///   sld_string 
     47  //  [DEFAULT]=sld_string= 1.0e-06 [1/A^(2)] 
     48  double sld_string; 
     49  ///   sld_solv 
     50  //  [DEFAULT]=sld_solv= 6.3e-06 [1/A^(2)] 
     51  double sld_solv; 
     52  /// Background 
     53  //  [DEFAULT]=background=0 
     54  double background; 
    5455 
    5556}PeralNecklaceParameters; 
    5657 
    57 double pearl_necklace_kernel(double dq[], double q); 
    5858 
    59 /// 1D scattering function 
    60 double pearl_necklace_analytical_1D(PeralNecklaceParameters *pars, double q); 
    6159 
    62 /// 2D scattering function 
    63 double pearl_necklace_analytical_2D(PeralNecklaceParameters *pars, double q, double phi); 
    64 double pearl_necklace_analytical_2DXY(PeralNecklaceParameters *pars, double qx, double qy); 
     60class PearlNecklaceModel{ 
     61public: 
     62  // Model parameters 
     63  /// Scale factor 
     64  //  [DEFAULT]=scale= 1.0 
     65  Parameter scale; 
     66  /// radius [A] 
     67  //  [DEFAULT]=radius=80.0 [A] 
     68  Parameter radius; 
     69  /// edge_separation 
     70  //  [DEFAULT]=edge_separation= 350 [A] 
     71  Parameter edge_separation; 
     72  /// thick_string [A] 
     73  //  [DEFAULT]=thick_string= 2.5 [A] 
     74  Parameter thick_string; 
     75  /// num_pearls 
     76  //  [DEFAULT]=num_pearls= 3 
     77  Parameter num_pearls; 
     78  /// sld_pearl 
     79  //  [DEFAULT]=sld_pearl= 1.0e-06 [1/A^(2)] 
     80  Parameter sld_pearl; 
     81  /// sld_string 
     82  //  [DEFAULT]=sld_string= 1.0e-06 [1/A^(2)] 
     83  Parameter sld_string; 
     84  /// sld_solv 
     85  //  [DEFAULT]=sld_solv= 6.3e-06 [1/A^(2)] 
     86  Parameter sld_solv; 
     87  /// Background 
     88  //  [DEFAULT]=background=0 
     89  Parameter background; 
     90 
     91  // Constructor 
     92  PearlNecklaceModel(); 
     93 
     94  // Operators to get I(Q) 
     95  double operator()(double q); 
     96  double operator()(double qx, double qy); 
     97  double calculate_ER(); 
     98  double evaluate_rphi(double q, double phi); 
     99}; 
    65100 
    66101#endif 

  • sansmodels/src/c_extensions/sld_cal.h

    r67424cd r0ba3b08  
    11#if !defined(sld_cal_h) 
    22#define sld_cal_h 
     3#include "parameters.hh" 
    34 
    45/** 
    … …  
    1213 
    1314 **/ 
    14 typedef struct { 
    15     /// fun_type 
    16     //  [DEFAULT]=fun_type=0 
    17     double fun_type; 
    1815 
    19     /// npts_inter 
    20     //  [DEFAULT]=npts_inter= 21 
    21     double npts_inter; 
     16class SLDCalFunc{ 
     17public: 
     18  // Model parameters 
     19  /// fun_type 
     20  //  [DEFAULT]=fun_type=0 
     21  Parameter fun_type; 
    2222 
    23     /// shell_num 
    24     //  [DEFAULT]=shell_num= 0 
    25     double shell_num; 
     23  /// npts_inter 
     24  //  [DEFAULT]=npts_inter= 21 
     25  Parameter npts_inter; 
    2626 
    27     /// nu_inter 
    28     //  [DEFAULT]=nu_inter= 2.5 
    29     double nu_inter; 
     27  /// shell_num 
     28  //  [DEFAULT]=shell_num= 0 
     29  Parameter shell_num; 
    3030 
    31     /// sld_left [1/A^(2)] 
    32     //  [DEFAULT]=sld_left= 0 [1/A^(2)] 
    33     double sld_left; 
     31  /// nu_inter 
     32  //  [DEFAULT]=nu_inter= 2.5 
     33  Parameter nu_inter; 
    3434 
    35     /// sld_right [1/A^(2)] 
    36     //  [DEFAULT]=sld_right= 0 [1/A^(2)] 
    37     double sld_right; 
     35  /// sld_left [1/A^(2)] 
     36  //  [DEFAULT]=sld_left= 0 [1/A^(2)] 
     37  Parameter sld_left; 
    3838 
    39 } SLDCalParameters; 
     39  /// sld_right [1/A^(2)] 
     40  //  [DEFAULT]=sld_right= 0 [1/A^(2)] 
     41  Parameter sld_right; 
     42 
     43  // Constructor 
     44  SLDCalFunc(); 
     45 
     46  // Operators to get SLD 
     47  double operator()(double q); 
     48  double operator()(double qx, double qy); 
     49  double calculate_ER(); 
     50  double evaluate_rphi(double q, double phi); 
     51}; 
    4052 
    4153 
    4254 
    43 /// 1D function 
    44 double sld_cal_analytical_1D(SLDCalParameters *pars, double q); 
    45  
    46 /// 2D function 
    47 double sld_cal_analytical_2D(SLDCalParameters *pars, double q, double phi); 
    48 double sld_cal_analytical_2DXY(SLDCalParameters *pars, double qx, double qy); 
    49  
    5055#endif 

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