Changeset 6e10cff in sasview


Ignore:
Timestamp:
Jan 5, 2012 11:21:32 AM (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:
642a025
Parents:
37805e9
Message:

cscrystal model refactor

Location:
sansmodels/src
Files:
2 deleted
3 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/c_extensions/sc.h

    r67424cd r6e10cff  
    11#if !defined(sc_h) 
    22#define sc_h 
     3#include "parameters.hh" 
    34 
    45/** 
    56 * Structure definition for SC_ParaCrystal parameters 
    67 */ 
    7  //[PYTHONCLASS] = SCCrystalModel 
    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 simple cubic structure. 
    17  //                             [Simple 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> 
    30  
    31 typedef struct { 
    32     /// Scale factor 
    33     //  [DEFAULT]=scale= 1.0 
    34     double scale; 
    35  
    36     /// Nearest neighbor distance [A] 
    37     //  [DEFAULT]=dnn=220.0 [A] 
    38     double dnn; 
    39  
    40     /// Paracrystal distortion factor g 
    41     //  [DEFAULT]=d_factor=0.06 
    42     double d_factor; 
    43  
    44     /// Radius of sphere [A] 
    45     //  [DEFAULT]=radius=40.0 [A] 
    46     double radius; 
    47  
    48     /// sldSph [1/A^(2)] 
    49     //  [DEFAULT]=sldSph= 3.0e-6 [1/A^(2)] 
    50     double sldSph; 
    51  
    52     /// sldSolv [1/A^(2)] 
    53     //  [DEFAULT]=sldSolv= 6.3e-6 [1/A^(2)] 
    54     double sldSolv; 
    55  
    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; 
    68  
    69 } SCParameters; 
     8//[PYTHONCLASS] = SCCrystalModel 
     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 simple cubic structure. 
     18//                              [Simple 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> 
    7031 
    7132 
    7233 
    73 /// 1D scattering function 
    74 double sc_analytical_1D(SCParameters *pars, double q); 
     34class SCCrystalModel{ 
     35public: 
     36  // Model parameters 
     37  /// Scale factor 
     38  //  [DEFAULT]=scale= 1.0 
     39  Parameter scale; 
    7540 
    76 /// 2D scattering function 
    77 double sc_analytical_2D(SCParameters *pars, double q, double phi); 
    78 double sc_analytical_2DXY(SCParameters *pars, double qx, double qy); 
    79 double sc_analytical_2D_scaled(SCParameters *pars, double q, double q_x, double q_y); 
     41  /// Nearest neighbor distance [A] 
     42  //  [DEFAULT]=dnn=220.0 [A] 
     43  Parameter dnn; 
     44 
     45  /// Paracrystal distortion factor g 
     46  //  [DEFAULT]=d_factor=0.06 
     47  Parameter d_factor; 
     48 
     49  /// Radius of sphere [A] 
     50  //  [DEFAULT]=radius=40.0 [A] 
     51  Parameter radius; 
     52 
     53  /// sldSph [1/A^(2)] 
     54  //  [DEFAULT]=sldSph= 3.0e-6 [1/A^(2)] 
     55  Parameter sldSph; 
     56 
     57  /// sldSolv [1/A^(2)] 
     58  //  [DEFAULT]=sldSolv= 6.3e-6 [1/A^(2)] 
     59  Parameter sldSolv; 
     60 
     61  /// Incoherent Background [1/cm] 
     62  //  [DEFAULT]=background=0 [1/cm] 
     63  Parameter background; 
     64  /// Orientation of the a1 axis w/respect incoming beam [deg] 
     65  //  [DEFAULT]=theta=0.0 [deg] 
     66  Parameter theta; 
     67  /// Orientation of the a2 in the plane of the detector [deg] 
     68  //  [DEFAULT]=phi=0.0 [deg] 
     69  Parameter phi; 
     70  /// Orientation of the a3 in the plane of the detector [deg] 
     71  //  [DEFAULT]=psi=0.0 [deg] 
     72  Parameter psi; 
     73 
     74  // Constructor 
     75  SCCrystalModel(); 
     76 
     77  // Operators to get I(Q) 
     78  double operator()(double q); 
     79  double operator()(double qx, double qy); 
     80  double calculate_ER(); 
     81  double evaluate_rphi(double q, double phi); 
     82}; 
    8083 
    8184#endif 

  • sansmodels/src/c_models/models.hh

    r37805e9 r6e10cff  
    2727using namespace std; 
    2828 
    29  
    30  
    31  
    32  
    33  
    34 class SCCrystalModel{ 
    35 public: 
    36         // Model parameters 
    37         Parameter scale; 
    38         Parameter dnn; 
    39         Parameter d_factor; 
    40         Parameter radius; 
    41         Parameter sldSph; 
    42         Parameter sldSolv; 
    43         Parameter background; 
    44         Parameter theta; 
    45         Parameter phi; 
    46         Parameter psi; 
    47  
    48         // Constructor 
    49         SCCrystalModel(); 
    50  
    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 }; 
    5729 
    5830 

  • sansmodels/src/c_models/sc.cpp

    r67424cd r6e10cff  
    2525#include <stdio.h> 
    2626using namespace std; 
     27#include "sc.h" 
    2728 
    2829extern "C" { 
    29         #include "libSphere.h" 
    30         #include "sc.h" 
     30#include "libSphere.h" 
     31} 
     32 
     33// Convenience structure 
     34typedef struct { 
     35  double scale; 
     36  double dnn; 
     37  double d_factor; 
     38  double radius; 
     39  double sldSph; 
     40  double sldSolv; 
     41  double background; 
     42  double theta; 
     43  double phi; 
     44  double psi; 
     45} SCParameters; 
     46 
     47/** 
     48 * Function to evaluate 2D scattering function 
     49 * @param pars: parameters of the SCCrystalModel 
     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 sc_analytical_2D_scaled(SCParameters *pars, double q, double q_x, double q_y) { 
     56  double a3_x, a3_y, a3_z, a2_x, a2_y, a1_x, a1_y; 
     57  double q_z; 
     58  double alpha, cos_val_a3, cos_val_a2, cos_val_a1; 
     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; 
     63 
     64  double dp[5]; 
     65  //convert angle degree to radian 
     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  dp[0] = 1.0; 
     70  dp[1] = pars->radius; 
     71  dp[2] = pars->sldSph; 
     72  dp[3] = pars->sldSolv; 
     73  dp[4] = 0.0; 
     74 
     75 
     76  aa = pars->dnn; 
     77  Da = pars->d_factor*aa; 
     78  qDa_2 = pow(q*Da,2.0); 
     79 
     80  latticeScale = (4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa,3.0); 
     81  /// Angles here are respect to detector coordinate instead of against q coordinate(PRB 36, 3, 1754) 
     82  // a3 axis orientation 
     83  a3_x = sin(theta) * cos(phi);//negative sign here??? 
     84  a3_y = sin(theta) * sin(phi); 
     85  a3_z = cos(theta); 
     86 
     87  // q vector 
     88  q_z = 0.0; 
     89 
     90  // Compute the angle btw vector q and the a3 axis 
     91  cos_val_a3 = a3_x*q_x + a3_y*q_y + a3_z*q_z; 
     92  alpha = acos(cos_val_a3); 
     93  //alpha = acos(cos_val_a3); 
     94  a3_dot_q = aa*q*cos_val_a3; 
     95  // a1 axis orientation 
     96  a1_x = sin(psi); 
     97  a1_y = cos(psi); 
     98 
     99  cos_val_a1 = a1_x*q_x + a1_y*q_y; 
     100  a1_dot_q = aa*q*cos_val_a1*sin(alpha); 
     101 
     102  // a2 axis orientation 
     103  a2_x = sqrt(1.0-sin(theta)*cos(phi))*cos(psi); 
     104  a2_y = sqrt(1.0-sin(theta)*cos(phi))*sin(psi); 
     105  // a2 axis 
     106  cos_val_a2 =  sin(acos(cos_val_a1));//a2_x*q_x + a2_y*q_y; 
     107  a2_dot_q = aa*q*cos_val_a2*sin(alpha); 
     108 
     109  // The following test should always pass 
     110  if (fabs(cos_val_a3)>1.0) { 
     111    printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
     112    return 0; 
     113  } 
     114  // Call Zq=Z1*Z2*Z3 
     115  Zq = (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a1_dot_q)+exp(-qDa_2)); 
     116  Zq = Zq * (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a2_dot_q)+exp(-qDa_2)); 
     117  Zq = Zq * (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a3_dot_q)+exp(-qDa_2)); 
     118 
     119  // Use SphereForm directly from libigor 
     120  answer = SphereForm(dp,q)*Zq; 
     121 
     122  //consider scales 
     123  answer *= latticeScale * pars->scale; 
     124 
     125  // This FIXES a singualrity the kernel in libigor. 
     126  if ( answer == INFINITY || answer == NAN){ 
     127    answer = 0.0; 
     128  } 
     129 
     130  // add background 
     131  answer += pars->background; 
     132 
     133  return answer; 
     134} 
     135 
     136/** 
     137 * Function to evaluate 2D scattering function 
     138 * @param pars: parameters of the SC_ParaCrystal 
     139 * @param q: q-value 
     140 * @return: function value 
     141 */ 
     142static double sc_analytical_2DXY(SCParameters *pars, double qx, double qy){ 
     143  double q; 
     144  q = sqrt(qx*qx+qy*qy); 
     145  return sc_analytical_2D_scaled(pars, q, qx/q, qy/q); 
    31146} 
    32147 
    33148SCCrystalModel :: SCCrystalModel() { 
    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); 
     149  scale      = Parameter(1.0); 
     150  dnn           = Parameter(220.0); 
     151  d_factor = Parameter(0.06); 
     152  radius     = Parameter(40.0, true); 
     153  radius.set_min(0.0); 
     154  sldSph   = Parameter(3.0e-6); 
     155  sldSolv   = Parameter(6.3e-6); 
     156  background = Parameter(0.0); 
     157  theta  = Parameter(0.0, true); 
     158  phi    = Parameter(0.0, true); 
     159  psi    = Parameter(0.0, true); 
    45160} 
    46161 
    … …  
    52167 */ 
    53168double SCCrystalModel :: 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 = SC_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(); 
     169  double dp[7]; 
     170 
     171  // Fill parameter array for IGOR library 
     172  // Add the background after averaging 
     173  dp[0] = scale(); 
     174  dp[1] = dnn(); 
     175  dp[2] = d_factor(); 
     176  dp[3] = radius(); 
     177  dp[4] = sldSph(); 
     178  dp[5] = sldSolv(); 
     179  dp[6] = 0.0; 
     180 
     181  // Get the dispersion points for the radius 
     182  vector<WeightPoint> weights_rad; 
     183  radius.get_weights(weights_rad); 
     184 
     185  // Perform the computation, with all weight points 
     186  double sum = 0.0; 
     187  double norm = 0.0; 
     188  double vol = 0.0; 
     189  double result; 
     190 
     191  // Loop over radius weight points 
     192  for(size_t i=0; i<weights_rad.size(); i++) { 
     193    dp[3] = weights_rad[i].value; 
     194 
     195    //Un-normalize SphereForm by volume 
     196    result = SC_ParaCrystal(dp, q) * pow(weights_rad[i].value,3); 
     197    // This FIXES a singualrity the kernel in libigor. 
     198    if ( result == INFINITY || result == NAN){ 
     199      result = 0.0; 
     200    } 
     201    sum += weights_rad[i].weight 
     202        * result; 
     203    //Find average volume 
     204    vol += weights_rad[i].weight 
     205        * pow(weights_rad[i].value,3); 
     206 
     207    norm += weights_rad[i].weight; 
     208  } 
     209 
     210  if (vol != 0.0 && norm != 0.0) { 
     211    //Re-normalize by avg volume 
     212    sum = sum/(vol/norm);} 
     213  return sum/norm + background(); 
    99214} 
    100215 
    … …  
    106221 */ 
    107222double SCCrystalModel :: operator()(double qx, double qy) { 
    108         SCParameters 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                                                                                 * sc_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                                         //Find average volume 
    171                                         //vol += weights_rad[i].weight 
    172                                         //      * pow(weights_rad[i].value,3); 
    173                                         //Find norm for volume 
    174                                         //norm_vol += weights_rad[i].weight; 
    175  
    176                                         norm += weights_rad[i].weight 
    177                                                         * weights_theta[j].weight 
    178                                                         * weights_phi[k].weight 
    179                                                         * weights_psi[l].weight; 
    180                                 } 
    181                         } 
    182                 } 
    183         } 
    184         // Averaging in theta needs an extra normalization 
    185         // factor to account for the sin(theta) term in the 
    186         // integration (see documentation). 
    187         if (weights_theta.size()>1) norm = norm / asin(1.0); 
    188  
    189         if (vol != 0.0 && norm_vol != 0.0) { 
    190                 //Re-normalize by avg volume 
    191                 sum = sum/(vol/norm_vol);} 
    192  
    193         return sum/norm + background(); 
     223  SCParameters dp; 
     224  dp.scale      = scale(); 
     225  dp.dnn   = dnn(); 
     226  dp.d_factor   = d_factor(); 
     227  dp.radius  = radius(); 
     228  dp.sldSph   = sldSph(); 
     229  dp.sldSolv   = sldSolv(); 
     230  dp.background = 0.0; 
     231  dp.theta  = theta(); 
     232  dp.phi    = phi(); 
     233  dp.psi    = psi(); 
     234 
     235  // Get the dispersion points for the radius 
     236  vector<WeightPoint> weights_rad; 
     237  radius.get_weights(weights_rad); 
     238 
     239  // Get angular averaging for theta 
     240  vector<WeightPoint> weights_theta; 
     241  theta.get_weights(weights_theta); 
     242 
     243  // Get angular averaging for phi 
     244  vector<WeightPoint> weights_phi; 
     245  phi.get_weights(weights_phi); 
     246 
     247  // Get angular averaging for psi 
     248  vector<WeightPoint> weights_psi; 
     249  psi.get_weights(weights_psi); 
     250 
     251  // Perform the computation, with all weight points 
     252  double sum = 0.0; 
     253  double norm = 0.0; 
     254  double norm_vol = 0.0; 
     255  double vol = 0.0; 
     256  double pi = 4.0*atan(1.0); 
     257  // Loop over radius weight points 
     258  for(size_t i=0; i<weights_rad.size(); i++) { 
     259    dp.radius = weights_rad[i].value; 
     260    // Average over theta distribution 
     261    for(size_t j=0; j< weights_theta.size(); j++) { 
     262      dp.theta = weights_theta[j].value; 
     263      // Average over phi distribution 
     264      for(size_t k=0; k< weights_phi.size(); k++) { 
     265        dp.phi = weights_phi[k].value; 
     266        // Average over phi distribution 
     267        for(size_t l=0; l< weights_psi.size(); l++) { 
     268          dp.psi = weights_psi[l].value; 
     269          //Un-normalize SphereForm by volume 
     270          double _ptvalue = weights_rad[i].weight 
     271              * weights_theta[j].weight 
     272              * weights_phi[k].weight 
     273              * weights_psi[l].weight 
     274              * sc_analytical_2DXY(&dp, qx, qy); 
     275          //* pow(weights_rad[i].value,3.0); 
     276          // Consider when there is infinte or nan. 
     277          if ( _ptvalue == INFINITY || _ptvalue == NAN){ 
     278            _ptvalue = 0.0; 
     279          } 
     280          if (weights_theta.size()>1) { 
     281            _ptvalue *= fabs(sin(weights_theta[j].value*pi/180.0)); 
     282          } 
     283          sum += _ptvalue; 
     284          // This model dose not need the volume of spheres correction!!! 
     285          //Find average volume 
     286          //vol += weights_rad[i].weight 
     287          //    * pow(weights_rad[i].value,3); 
     288          //Find norm for volume 
     289          //norm_vol += weights_rad[i].weight; 
     290 
     291          norm += weights_rad[i].weight 
     292              * weights_theta[j].weight 
     293              * weights_phi[k].weight 
     294              * weights_psi[l].weight; 
     295        } 
     296      } 
     297    } 
     298  } 
     299  // Averaging in theta needs an extra normalization 
     300  // factor to account for the sin(theta) term in the 
     301  // integration (see documentation). 
     302  if (weights_theta.size()>1) norm = norm / asin(1.0); 
     303 
     304  if (vol != 0.0 && norm_vol != 0.0) { 
     305    //Re-normalize by avg volume 
     306    sum = sum/(vol/norm_vol);} 
     307 
     308  return sum/norm + background(); 
    194309} 
    195310 
    … …  
    202317 */ 
    203318double SCCrystalModel :: evaluate_rphi(double q, double phi) { 
    204         return (*this).operator()(q); 
     319  return (*this).operator()(q); 
    205320} 
    206321 
    … …  
    210325 */ 
    211326double SCCrystalModel :: calculate_ER() { 
    212         //NOT implemented yet!!! 
    213         return 0.0; 
    214 } 
     327  //NOT implemented yet!!! 
     328  return 0.0; 
     329} 
Note: See TracChangeset for help on using the changeset viewer.