Changeset 3c102d4 in sasview for sansmodels/src


Ignore:
Timestamp:
Aug 18, 2009 1:56:49 PM (15 years ago)
Author:
Jae Cho <jhjcho@…>
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:
eddff027
Parents:
7d11b81
Message:

fixed problems in 2d

Location:
sansmodels/src/sans/models
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • sansmodels/src/sans/models/TriaxialEllipsoidModel.py

    r975ec8e r3c102d4  
    3333        List of default parameters: 
    3434         scale           = 1.0  
    35          semi_axisB      = 35.0 [A] 
    36          semi_axisA      = 100.0 [A] 
     35         semi_axisA      = 35.0 [A] 
     36         semi_axisB      = 100.0 [A] 
    3737         semi_axisC      = 400.0 [A] 
    3838         contrast        = 5.3e-006 [1/A²] 
    3939         background      = 0.0 [1/cm] 
    40          axis_theta      = 0.0 [rad] 
    41          axis_phi        = 0.0 [rad] 
     40         axis_theta      = 1.0 [rad] 
     41         axis_phi        = 1.0 [rad] 
    4242         axis_psi        = 0.0 [rad] 
    4343 
     
    6161        self.details = {} 
    6262        self.details['scale'] = ['', None, None] 
     63        self.details['semi_axisA'] = ['[A]', None, None] 
    6364        self.details['semi_axisB'] = ['[A]', None, None] 
    64         self.details['semi_axisA'] = ['[A]', None, None] 
    6565        self.details['semi_axisC'] = ['[A]', None, None] 
    6666        self.details['contrast'] = ['[1/A²]', None, None] 
  • sansmodels/src/sans/models/c_extensions/parallelepiped.c

    r2cb89e7 r3c102d4  
    9898 */ 
    9999double parallelepiped_analytical_2D_scaled(ParallelepipedParameters *pars, double q, double q_x, double q_y) { 
    100         double cparallel_x, cparallel_y, cparallel_z, bparallel_x, bparallel_y, parallel_x, parallel_y, parallel_z; 
     100        double parallel_x, parallel_y, parallel_z; 
    101101        double q_z; 
    102102        double alpha, vol, cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC; 
     
    111111 
    112112    // parallelepiped c axis orientation 
    113     cparallel_x = sin(pars->parallel_theta) * cos(pars->parallel_phi); 
    114     cparallel_y = sin(pars->parallel_theta) * sin(pars->parallel_phi); 
    115     cparallel_z = cos(pars->parallel_theta); 
     113    parallel_x = sin(pars->parallel_theta) * cos(pars->parallel_phi); 
     114    parallel_y = sin(pars->parallel_theta) * sin(pars->parallel_phi); 
     115    parallel_z = cos(pars->parallel_theta); 
    116116 
    117117    // q vector 
     
    120120    // Compute the angle btw vector q and the 
    121121    // axis of the parallelepiped 
    122     cos_val_c = cparallel_x*q_x + cparallel_y*q_y + cparallel_z*q_z; 
     122    cos_val_c = parallel_x*q_x + parallel_y*q_y + parallel_z*q_z; 
    123123    alpha = acos(cos_val_c); 
    124124 
     
    126126    parallel_x = -(1-sin(pars->parallel_theta)*sin(pars->parallel_phi))*sin(pars->parallel_psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)*sin(pars->parallel_psi); 
    127127    parallel_y = (1-sin(pars->parallel_theta)*sin(pars->parallel_phi))*cos(pars->parallel_psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)*cos(pars->parallel_psi); 
     128 
     129    //parallel_x = -(1-sin(pars->parallel_theta)*sin(pars->parallel_phi))*sin(pars->parallel_psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)*sin(pars->parallel_psi); 
     130    //parallel_y = (1-sin(pars->parallel_theta)*sin(pars->parallel_phi))*cos(pars->parallel_psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)*cos(pars->parallel_psi); 
    128131    cos_val_a = (parallel_x*q_x + parallel_y*q_y); 
    129132 
     
    131134 
    132135    // parallelepiped b axis orientation 
    133     bparallel_x = (1-sin(pars->parallel_theta)*cos(pars->parallel_phi))*cos(pars->parallel_psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)* cos(pars->parallel_psi); 
    134     bparallel_y = (1-sin(pars->parallel_theta)*cos(pars->parallel_phi))*sin(pars->parallel_psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)* sin(pars->parallel_psi); 
     136    parallel_x = (1-sin(pars->parallel_theta)*cos(pars->parallel_phi))*cos(pars->parallel_psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)* cos(pars->parallel_psi); 
     137    parallel_y = (1-sin(pars->parallel_theta)*cos(pars->parallel_phi))*sin(pars->parallel_psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)* sin(pars->parallel_psi); 
    135138    // axis of the parallelepiped 
    136     cos_val_b = (bparallel_x*q_x + bparallel_y*q_y) ; 
     139    cos_val_b = (parallel_x*q_x + parallel_y*q_y) ; 
    137140 
    138141 
  • sansmodels/src/sans/models/c_extensions/stacked_disks.c

    r5068697 r3c102d4  
    1919double stacked_disks_analytical_1D(StackedDisksParameters *pars, double q) { 
    2020        double dp[10]; 
    21          
     21 
    2222        // Fill paramater array 
    2323        dp[0] = pars->scale; 
    2424        dp[1] = pars->radius; 
    25         dp[2] = pars->length; 
    26         dp[3] = pars->thickness; 
     25        dp[2] = pars->core_thick; 
     26        dp[3] = pars->layer_thick; 
    2727        dp[4] = pars->core_sld; 
    2828        dp[5] = pars->layer_sld; 
    2929        dp[6] = pars->solvent_sld; 
    30         dp[7] = pars->nlayers; 
    31         dp[8] = pars->spacing; 
     30        dp[7] = pars->n_stacking; 
     31        dp[8] = pars->sigma_d; 
    3232        dp[9] = pars->background; 
    3333 
    3434        // Call library function to evaluate model 
    35         return StackedDiscs(dp, q);      
     35        return StackedDiscs(dp, q); 
    3636} 
    3737/** 
     
    4545        q = sqrt(qx*qx+qy*qy); 
    4646    return stacked_disks_analytical_2D_scaled(pars, q, qx/q, qy/q); 
    47 }  
     47} 
    4848 
    4949 
     
    5757double stacked_disks_analytical_2D(StackedDisksParameters *pars, double q, double phi) { 
    5858    return stacked_disks_analytical_2D_scaled(pars, q, cos(phi), sin(phi)); 
    59 }  
    60          
     59} 
     60 
    6161/** 
    6262 * Function to evaluate 2D scattering function 
     
    7171        double q_z; 
    7272        double alpha, vol, cos_val; 
    73         double d, dum; 
     73        double d, dum, halfheight; 
    7474        double answer; 
    75          
    76          
     75 
     76 
    7777 
    7878    // parallelepiped orientation 
     
    8080    cyl_y = sin(pars->axis_theta) * sin(pars->axis_phi); 
    8181    cyl_z = cos(pars->axis_theta); 
    82       
     82 
    8383    // q vector 
    8484    q_z = 0; 
    85          
     85 
    8686    // Compute the angle btw vector q and the 
    8787    // axis of the parallelepiped 
    8888    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z; 
    89      
     89 
    9090    // The following test should always pass 
    9191    if (fabs(cos_val)>1.0) { 
     
    9393        return 0; 
    9494    } 
    95      
     95 
    9696    // Note: cos(alpha) = 0 and 1 will get an 
    9797    // undefined value from Stackdisc_kern 
     
    9999 
    100100        // Call the IGOR library function to get the kernel 
    101         dum =0.1; 
    102         d= 0.1; 
     101        d = 2 * pars->layer_thick + pars->core_thick; 
     102        halfheight = pars->core_thick/2.0; 
     103        dum =alpha ; 
    103104        answer = Stackdisc_kern(q, pars->radius, pars->core_sld,pars->layer_sld, 
    104                 pars->solvent_sld,pars->length,pars->thickness, dum, pars->spacing, d,pars->nlayers); 
    105          
     105                pars->solvent_sld, halfheight, pars->layer_thick, dum, pars->sigma_d, d, pars->n_stacking)/sin(alpha); 
     106 
    106107        // Multiply by contrast^2 
    107108        //answer *= pars->contrast*pars->contrast; 
    108          
     109 
    109110        //normalize by staked disks volume 
    110         //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vparallel 
    111     vol = acos(-1.0) * pars->radius * pars->radius * pars->length; 
    112         answer *= vol; 
    113          
     111    vol = acos(-1.0) * pars->radius * pars->radius * d * pars->n_stacking; 
     112        answer /= vol; 
     113 
    114114        //convert to [cm-1] 
    115115        answer *= 1.0e8; 
    116          
     116 
    117117        //Scale 
    118118        answer *= pars->scale; 
    119          
     119 
    120120        // add in the background 
    121121        answer += pars->background; 
    122          
     122 
    123123        return answer; 
    124124} 
    125      
    126125 
     126 
  • sansmodels/src/sans/models/c_extensions/triaxial_ellipsoid.c

    r975ec8e r3c102d4  
    3535        double t,a,b,c; 
    3636        double kernel; 
    37         double pi = acos(-1.0); 
     37        double pi = 4.0*atan(1.0); 
    3838 
    3939        a = pars->semi_axisA ; 
     
    4242 
    4343        t = q * sqrt(a*a*cos(nu)*cos(nu)+b*b*sin(nu)*sin(nu)*sin(alpha)*sin(alpha)+c*c*cos(alpha)*cos(alpha)); 
    44         if (t==0){ 
     44        if (t==0.0){ 
    4545                kernel  = 1.0; 
    4646        }else{ 
    47                 kernel  = 3*(sin(t)-t*cos(t))/(t*t*t); 
     47                kernel  = 3.0*(sin(t)-t*cos(t))/(t*t*t); 
    4848        } 
    4949        return kernel*kernel; 
     
    137137        //normalize by cylinder volume 
    138138        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 
    139     vol = 4/3 * pi * pars->semi_axisA * pars->semi_axisB * pars->semi_axisC; 
     139    vol = 4.0* pi/3.0 * pars->semi_axisA * pars->semi_axisB * pars->semi_axisC; 
    140140        answer *= vol; 
    141  
    142141        //convert to [cm-1] 
    143142        answer *= 1.0e8; 
    144  
    145143        //Scale 
    146144        answer *= pars->scale; 
  • sansmodels/src/sans/models/c_extensions/triaxial_ellipsoid.h

    r975ec8e r3c102d4  
    1717    //  [DEFAULT]=scale=1.0 
    1818    double scale; 
    19     /// semi -axis B of the triaxial_ellipsoid [A] 
    20     //  [DEFAULT]=semi_axisB= 35.0 [A] 
     19    /// semi -axis A of the triaxial_ellipsoid [A] 
     20    //  [DEFAULT]=semi_axisA= 35.0 [A] 
    2121    double semi_axisA; 
    2222    /// semi -axis B of the triaxial_ellipsoid [A] 
    23     //  [DEFAULT]=semi_axisA=100.0 [A] 
     23    //  [DEFAULT]=semi_axisB=100.0 [A] 
    2424    double semi_axisB; 
    2525        /// semi -axis C of the triaxial_ellipsoid [A] 
     
    3333        double background; 
    3434    /// Orientation of the triaxial_ellipsoid axis w/respect incoming beam [rad] 
    35     //  [DEFAULT]=axis_theta=0.0 [rad] 
     35    //  [DEFAULT]=axis_theta=1.0 [rad] 
    3636    double axis_theta; 
    3737    /// Orientation of the triaxial_ellipsoid in the plane of the detector [rad] 
    38     //  [DEFAULT]=axis_phi=0.0 [rad] 
     38    //  [DEFAULT]=axis_phi=1.0 [rad] 
    3939    double axis_phi; 
    4040    /// Orientation of the cross section of the triaxial_ellipsoid in the plane of the detector [rad] 
  • sansmodels/src/sans/models/c_models/CTriaxialEllipsoidModel.cpp

    r975ec8e r3c102d4  
    8888        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000)); 
    8989        PyDict_SetItemString(self->params,"axis_psi",Py_BuildValue("d",0.000000)); 
    90         PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",0.000000)); 
    91         PyDict_SetItemString(self->params,"semi_axisA",Py_BuildValue("d",100.000000)); 
    92         PyDict_SetItemString(self->params,"semi_axisB",Py_BuildValue("d",35.000000)); 
     90        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.000000)); 
     91        PyDict_SetItemString(self->params,"semi_axisA",Py_BuildValue("d",35.000000)); 
     92        PyDict_SetItemString(self->params,"semi_axisB",Py_BuildValue("d",100.000000)); 
    9393        PyDict_SetItemString(self->params,"semi_axisC",Py_BuildValue("d",400.000000)); 
    94         PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.000000)); 
     94        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",1.000000)); 
    9595        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000)); 
    9696        PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.000005)); 
  • sansmodels/src/sans/models/c_models/parallelepiped.cpp

    r975ec8e r3c102d4  
    3535ParallelepipedModel :: ParallelepipedModel() { 
    3636        scale      = Parameter(1.0); 
    37         short_edgeA     = Parameter(35.0, true); 
    38         short_edgeA.set_max(1.0); 
    39         longer_edgeB     = Parameter(75.0, true); 
    40         longer_edgeB.set_min(1.0); 
    41         longuest_edgeC     = Parameter(400.0, true); 
    42         longuest_edgeC.set_min(1.0); 
     37        short_a     = Parameter(35.0, true); 
     38        short_a.set_max(1.0); 
     39        long_b     = Parameter(75.0, true); 
     40        long_b.set_min(1.0); 
     41        longer_c     = Parameter(400.0, true); 
     42        longer_c.set_min(1.0); 
    4343        contrast   = Parameter(53.e-7); 
    4444        background = Parameter(0.0); 
     
    6060        // Add the background after averaging 
    6161        dp[0] = scale(); 
    62         dp[1] = short_edgeA(); 
    63         dp[2] = longer_edgeB(); 
    64         dp[3] = longuest_edgeC(); 
     62        dp[1] = short_a(); 
     63        dp[2] = long_b(); 
     64        dp[3] = longer_c(); 
    6565        dp[4] = contrast(); 
    6666        dp[5] = 0.0; 
    6767 
    6868        // Get the dispersion points for the short_edgeA 
    69         vector<WeightPoint> weights_short_edgeA; 
    70         short_edgeA.get_weights(weights_short_edgeA); 
     69        vector<WeightPoint> weights_short_a; 
     70        short_a.get_weights(weights_short_a); 
    7171 
    7272        // Get the dispersion points for the longer_edgeB 
    73         vector<WeightPoint> weights_longer_edgeB; 
    74         longer_edgeB.get_weights(weights_longer_edgeB); 
     73        vector<WeightPoint> weights_long_b; 
     74        long_b.get_weights(weights_long_b); 
    7575 
    7676        // Get the dispersion points for the longuest_edgeC 
    77         vector<WeightPoint> weights_longuest_edgeC; 
    78         longuest_edgeC.get_weights(weights_longuest_edgeC); 
     77        vector<WeightPoint> weights_longer_c; 
     78        longer_c.get_weights(weights_longer_c); 
    7979 
    8080 
     
    8585 
    8686        // Loop over short_edgeA weight points 
    87         for(int i=0; i< (int)weights_short_edgeA.size(); i++) { 
    88                 dp[1] = weights_short_edgeA[i].value; 
     87        for(int i=0; i< (int)weights_short_a.size(); i++) { 
     88                dp[1] = weights_short_a[i].value; 
    8989 
    9090                // Loop over longer_edgeB weight points 
    91                 for(int j=0; j< (int)weights_longer_edgeB.size(); j++) { 
    92                         dp[2] = weights_longer_edgeB[i].value; 
     91                for(int j=0; j< (int)weights_long_b.size(); j++) { 
     92                        dp[2] = weights_long_b[j].value; 
    9393 
    9494                        // Loop over longuest_edgeC weight points 
    95                         for(int k=0; k< (int)weights_longuest_edgeC.size(); k++) { 
    96                                 dp[3] = weights_longuest_edgeC[j].value; 
    97  
    98                                 sum += weights_short_edgeA[i].weight * weights_longer_edgeB[j].weight 
    99                                         * weights_longuest_edgeC[k].weight * Parallelepiped(dp, q); 
    100  
    101                                 norm += weights_short_edgeA[i].weight 
    102                                          * weights_longer_edgeB[j].weight * weights_longuest_edgeC[k].weight; 
     95                        for(int k=0; k< (int)weights_longer_c.size(); k++) { 
     96                                dp[3] = weights_longer_c[k].value; 
     97 
     98                                sum += weights_short_a[i].weight * weights_long_b[j].weight 
     99                                        * weights_longer_c[k].weight * Parallelepiped(dp, q); 
     100 
     101                                norm += weights_short_a[i].weight 
     102                                         * weights_long_b[j].weight * weights_longer_c[k].weight; 
    103103                        } 
    104104                } 
     
    116116        // Fill parameter array 
    117117        dp.scale      = scale(); 
    118         dp.short_edgeA   = short_edgeA(); 
    119         dp.longer_edgeB   = longer_edgeB(); 
    120         dp.longuest_edgeC  = longuest_edgeC(); 
     118        dp.short_a   = short_a(); 
     119        dp.long_b   = long_b(); 
     120        dp.longer_c  = longer_c(); 
    121121        dp.contrast   = contrast(); 
    122122        dp.background = 0.0; 
     
    128128 
    129129        // Get the dispersion points for the short_edgeA 
    130         vector<WeightPoint> weights_short_edgeA; 
    131         short_edgeA.get_weights(weights_short_edgeA); 
     130        vector<WeightPoint> weights_short_a; 
     131        short_a.get_weights(weights_short_a); 
    132132 
    133133        // Get the dispersion points for the longer_edgeB 
    134         vector<WeightPoint> weights_longer_edgeB; 
    135         longer_edgeB.get_weights(weights_longer_edgeB); 
     134        vector<WeightPoint> weights_long_b; 
     135        long_b.get_weights(weights_long_b); 
    136136 
    137137        // Get angular averaging for the longuest_edgeC 
    138         vector<WeightPoint> weights_longuest_edgeC; 
    139         longuest_edgeC.get_weights(weights_longuest_edgeC); 
     138        vector<WeightPoint> weights_longer_c; 
     139        longer_c.get_weights(weights_longer_c); 
    140140 
    141141        // Get angular averaging for theta 
     
    156156 
    157157        // Loop over radius weight points 
    158         for(int i=0; i< (int)weights_short_edgeA.size(); i++) { 
    159                 dp.short_edgeA = weights_short_edgeA[i].value; 
     158        for(int i=0; i< (int)weights_short_a.size(); i++) { 
     159                dp.short_a = weights_short_a[i].value; 
    160160 
    161161                // Loop over longer_edgeB weight points 
    162                 for(int j=0; j< (int)weights_longer_edgeB.size(); j++) { 
    163                         dp.longer_edgeB = weights_longer_edgeB[j].value; 
     162                for(int j=0; j< (int)weights_long_b.size(); j++) { 
     163                        dp.long_b = weights_long_b[j].value; 
    164164 
    165165                        // Average over longuest_edgeC distribution 
    166                         for(int k=0; k< (int)weights_longuest_edgeC.size(); k++) { 
    167                                 dp.longuest_edgeC = weights_longuest_edgeC[k].value; 
     166                        for(int k=0; k< (int)weights_longer_c.size(); k++) { 
     167                                dp.longer_c = weights_longer_c[k].value; 
    168168 
    169169                                // Average over theta distribution 
     
    179179                                                        dp.parallel_psi = weights_parallel_psi[n].value; 
    180180 
    181                                                         double _ptvalue = weights_short_edgeA[i].weight 
    182                                                                 * weights_longer_edgeB[j].weight 
    183                                                                 * weights_longuest_edgeC[k].weight 
     181                                                        double _ptvalue = weights_short_a[i].weight 
     182                                                                * weights_long_b[j].weight 
     183                                                                * weights_longer_c[k].weight 
    184184                                                                * weights_parallel_theta[l].weight 
    185185                                                                * weights_parallel_phi[m].weight 
     
    191191                                                        sum += _ptvalue; 
    192192 
    193                                                         norm += weights_short_edgeA[i].weight 
    194                                                                 * weights_longer_edgeB[j].weight 
    195                                                                 * weights_longuest_edgeC[k].weight 
     193                                                        norm += weights_short_a[i].weight 
     194                                                                * weights_long_b[j].weight 
     195                                                                * weights_longer_c[k].weight 
    196196                                                                * weights_parallel_theta[l].weight 
    197197                                                                * weights_parallel_phi[m].weight 
  • sansmodels/src/sans/models/c_models/triaxialellipsoid.cpp

    r975ec8e r3c102d4  
    3434TriaxialEllipsoidModel :: TriaxialEllipsoidModel() { 
    3535        scale      = Parameter(1.0); 
    36         semi_axisA     = Parameter(20.0, true); 
     36        semi_axisA     = Parameter(35.0, true); 
    3737        semi_axisA.set_min(0.0); 
    38         semi_axisB     = Parameter(20.0, true); 
     38        semi_axisB     = Parameter(100.0, true); 
    3939        semi_axisB.set_min(0.0); 
    4040        semi_axisC  = Parameter(400.0, true); 
     
    4242        contrast   = Parameter(5.3e-6); 
    4343        background = Parameter(0.0); 
    44         axis_theta  = Parameter(0.0, true); 
    45         axis_phi    = Parameter(0.0, true); 
     44        axis_theta  = Parameter(1.0, true); 
     45        axis_phi    = Parameter(1.0, true); 
    4646        axis_psi    = Parameter(0.0, true); 
    4747} 
     
    190190                                                                * weights_theta[l].weight 
    191191                                                                * weights_phi[m].weight 
    192                                                                 * weights_psi[m].weight; 
     192                                                                * weights_psi[n].weight; 
    193193                                                } 
    194194                                        } 
Note: See TracChangeset for help on using the changeset viewer.