Changes in / [aab23e3:f0ce6e2] in sasview


Ignore:
Location:
src/sas/sascalc/calculator
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • src/sas/sascalc/calculator/c_extensions/libfunc.c

    rb6c8abe r144e032a  
    9595// spintheta: angle (anti-clock-wise) between neutron spin(up) and x axis 
    9696// Note: all angles are in degrees. 
    97 polar_sld cal_msld(int isangle, double qx, double qy, double bn, 
     97void cal_msld(polar_sld *p_sld, int isangle, double qx, double qy, double bn, 
    9898                                double m01, double mtheta1, double mphi1, 
    9999                                double spinfraci, double spinfracf, double spintheta) 
     
    124124        double my = 0.0; 
    125125        double mz = 0.0; 
    126         polar_sld p_sld; 
    127         p_sld.uu = sld; 
    128         p_sld.dd = sld; 
    129         p_sld.re_ud = 0.0; 
    130         p_sld.im_ud = 0.0; 
    131         p_sld.re_du = 0.0; 
    132         p_sld.im_du = 0.0; 
     126        double uu = sld; 
     127        double dd = sld; 
     128        double re_ud = 0.0; 
     129        double im_ud = 0.0; 
     130        double re_du = 0.0; 
     131        double im_du = 0.0; 
    133132 
    134133        //No mag means no further calculation 
    135         if (isangle>0){ 
     134        if (isangle>0) { 
    136135                if (m_max < 1.0e-32){ 
    137                         p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu; 
    138                         p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd; 
    139                         return p_sld; 
    140                 } 
    141         } 
    142         else{ 
    143                 if (fabs(m_max)< 1.0e-32 && fabs(m_phi)< 1.0e-32 && fabs(m_theta)< 1.0e-32){ 
    144                         p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu; 
    145                         p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd; 
    146                         return p_sld; 
    147                 } 
    148         } 
    149  
    150         //These are needed because of the precision of inputs 
    151         if (in_spin < 0.0) in_spin = 0.0; 
    152         if (in_spin > 1.0) in_spin = 1.0; 
    153         if (out_spin < 0.0) out_spin = 0.0; 
    154         if (out_spin > 1.0) out_spin = 1.0; 
    155  
    156         if (q_x == 0.0) q_angle = pi / 2.0; 
    157         else q_angle = atan(q_y/q_x); 
    158         if (q_y < 0.0 && q_x < 0.0) q_angle -= pi; 
    159         else if (q_y > 0.0 && q_x < 0.0) q_angle += pi; 
    160  
    161         q_angle = pi/2.0 - q_angle; 
    162         if (q_angle > pi) q_angle -= 2.0 * pi; 
    163         else if (q_angle < -pi) q_angle += 2.0 * pi; 
    164  
    165         if (fabs(q_x) < 1.0e-16 && fabs(q_y) < 1.0e-16){ 
    166                 m_perp = 0.0; 
    167                 } 
    168         else { 
    169                 m_perp = m_max; 
    170                 } 
    171         if (is_angle > 0){ 
    172                 m_phi *= pi/180.0; 
    173                 m_theta *= pi/180.0; 
    174                 mx = m_perp * cos(m_theta) * cos(m_phi); 
    175                 my = m_perp * sin(m_theta); 
    176                 mz = -(m_perp * cos(m_theta) * sin(m_phi)); 
    177         } 
    178         else{ 
    179                 mx = m_perp; 
    180                 my = m_phi; 
    181                 mz = m_theta; 
    182         } 
    183         //ToDo: simplify these steps 
    184         // m_perp1 -m_perp2 
    185         m_perp_x = (mx) *  cos(q_angle); 
    186         m_perp_x -= (my) * sin(q_angle); 
    187         m_perp_y = m_perp_x; 
    188         m_perp_x *= cos(-q_angle); 
    189         m_perp_y *= sin(-q_angle); 
    190         m_perp_z = mz; 
    191  
    192         m_sigma_x = (m_perp_x * cos(-s_theta) - m_perp_y * sin(-s_theta)); 
    193         m_sigma_y = (m_perp_x * sin(-s_theta) + m_perp_y * cos(-s_theta)); 
    194         m_sigma_z = (m_perp_z); 
    195  
    196         //Find b 
    197         p_sld.uu -= m_sigma_x; 
    198         p_sld.dd += m_sigma_x; 
    199         p_sld.re_ud = m_sigma_y; 
    200         p_sld.re_du = m_sigma_y; 
    201         p_sld.im_ud = m_sigma_z; 
    202         p_sld.im_du = -m_sigma_z; 
    203  
    204         p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu; 
    205         p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd; 
    206  
    207         p_sld.re_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * p_sld.re_ud; 
    208         p_sld.im_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * p_sld.im_ud; 
    209         p_sld.re_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * p_sld.re_du; 
    210         p_sld.im_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * p_sld.im_du; 
    211  
    212         return p_sld; 
     136                        uu = sqrt(sqrt(in_spin * out_spin)) * uu; 
     137                        dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * dd; 
     138                } 
     139        } 
     140        else if (fabs(m_max)< 1.0e-32 && fabs(m_phi)< 1.0e-32 && fabs(m_theta)< 1.0e-32){ 
     141                        uu = sqrt(sqrt(in_spin * out_spin)) * uu; 
     142                        dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * dd; 
     143        } else { 
     144 
     145                //These are needed because of the precision of inputs 
     146                if (in_spin < 0.0) in_spin = 0.0; 
     147                if (in_spin > 1.0) in_spin = 1.0; 
     148                if (out_spin < 0.0) out_spin = 0.0; 
     149                if (out_spin > 1.0) out_spin = 1.0; 
     150 
     151                if (q_x == 0.0) q_angle = pi / 2.0; 
     152                else q_angle = atan(q_y/q_x); 
     153                if (q_y < 0.0 && q_x < 0.0) q_angle -= pi; 
     154                else if (q_y > 0.0 && q_x < 0.0) q_angle += pi; 
     155 
     156                q_angle = pi/2.0 - q_angle; 
     157                if (q_angle > pi) q_angle -= 2.0 * pi; 
     158                else if (q_angle < -pi) q_angle += 2.0 * pi; 
     159 
     160                if (fabs(q_x) < 1.0e-16 && fabs(q_y) < 1.0e-16){ 
     161                        m_perp = 0.0; 
     162                        } 
     163                else { 
     164                        m_perp = m_max; 
     165                        } 
     166                if (is_angle > 0){ 
     167                        m_phi *= pi/180.0; 
     168                        m_theta *= pi/180.0; 
     169                        mx = m_perp * cos(m_theta) * cos(m_phi); 
     170                        my = m_perp * sin(m_theta); 
     171                        mz = -(m_perp * cos(m_theta) * sin(m_phi)); 
     172                } 
     173                else{ 
     174                        mx = m_perp; 
     175                        my = m_phi; 
     176                        mz = m_theta; 
     177                } 
     178                //ToDo: simplify these steps 
     179                // m_perp1 -m_perp2 
     180                m_perp_x = (mx) *  cos(q_angle); 
     181                m_perp_x -= (my) * sin(q_angle); 
     182                m_perp_y = m_perp_x; 
     183                m_perp_x *= cos(-q_angle); 
     184                m_perp_y *= sin(-q_angle); 
     185                m_perp_z = mz; 
     186 
     187                m_sigma_x = (m_perp_x * cos(-s_theta) - m_perp_y * sin(-s_theta)); 
     188                m_sigma_y = (m_perp_x * sin(-s_theta) + m_perp_y * cos(-s_theta)); 
     189                m_sigma_z = (m_perp_z); 
     190 
     191                //Find b 
     192                uu -= m_sigma_x; 
     193                dd += m_sigma_x; 
     194                re_ud = m_sigma_y; 
     195                re_du = m_sigma_y; 
     196                im_ud = m_sigma_z; 
     197                im_du = -m_sigma_z; 
     198 
     199                uu = sqrt(sqrt(in_spin * out_spin)) * uu; 
     200                dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * dd; 
     201 
     202                re_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * re_ud; 
     203                im_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * im_ud; 
     204                re_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * re_du; 
     205                im_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * im_du; 
     206        } 
     207        p_sld->uu = uu; 
     208        p_sld->dd = dd; 
     209        p_sld->re_ud = re_ud; 
     210        p_sld->im_ud = im_ud; 
     211        p_sld->re_du = re_du; 
     212        p_sld->im_du = im_du; 
    213213} 
    214214 
  • src/sas/sascalc/calculator/c_extensions/libfunc.h

    rb6c8abe r144e032a  
    1818//double gamln(double x); 
    1919 
    20 polar_sld cal_msld(int isangle, double qx, double qy, double bn, double m01, double mtheta1,  
     20void cal_msld(polar_sld*, int isangle, double qx, double qy, double bn, double m01, double mtheta1,  
    2121                        double mphi1, double spinfraci, double spinfracf, double spintheta); 
    2222 
  • src/sas/sascalc/calculator/c_extensions/librefl.c

    ra1daf86 r144e032a  
    112112#endif // NEED_ERF 
    113113 
    114 complex cassign(real, imag) 
    115         double real, imag; 
    116 { 
    117         complex x; 
    118         x.re = real; 
    119         x.im = imag; 
    120         return x; 
    121 } 
    122  
    123  
    124 complex cplx_add(x,y) 
    125         complex x,y; 
    126 { 
    127         complex z; 
    128         z.re = x.re + y.re; 
    129         z.im = x.im + y.im; 
    130         return z; 
    131 } 
    132  
    133 complex rcmult(x,y) 
    134         double x; 
    135     complex y; 
    136 { 
    137         complex z; 
    138         z.re = x*y.re; 
    139         z.im = x*y.im; 
    140         return z; 
    141 } 
    142  
    143 complex cplx_sub(x,y) 
    144         complex x,y; 
    145 { 
    146         complex z; 
    147         z.re = x.re - y.re; 
    148         z.im = x.im - y.im; 
    149         return z; 
    150 } 
    151  
    152  
    153 complex cplx_mult(x,y) 
    154         complex x,y; 
    155 { 
    156         complex z; 
    157         z.re = x.re*y.re - x.im*y.im; 
    158         z.im = x.re*y.im + x.im*y.re; 
    159         return z; 
    160 } 
    161  
    162 complex cplx_div(x,y) 
    163         complex x,y; 
    164 { 
    165         complex z; 
    166         z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im); 
    167         z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im); 
    168         return z; 
    169 } 
    170  
    171 complex cplx_exp(b) 
    172         complex b; 
    173 { 
    174         complex z; 
     114void cassign(Cplx *x, double real, double imag) 
     115{ 
     116        x->re = real; 
     117        x->im = imag; 
     118} 
     119 
     120 
     121void cplx_add(Cplx *z, Cplx x, Cplx y) 
     122{ 
     123        z->re = x.re + y.re; 
     124        z->im = x.im + y.im; 
     125} 
     126 
     127void rcmult(Cplx *z, double x, Cplx y) 
     128{ 
     129        z->re = x*y.re; 
     130        z->im = x*y.im; 
     131} 
     132 
     133void cplx_sub(Cplx *z, Cplx x, Cplx y) 
     134{ 
     135        z->re = x.re - y.re; 
     136        z->im = x.im - y.im; 
     137} 
     138 
     139 
     140void cplx_mult(Cplx *z, Cplx x, Cplx y) 
     141{ 
     142        z->re = x.re*y.re - x.im*y.im; 
     143        z->im = x.re*y.im + x.im*y.re; 
     144} 
     145 
     146void cplx_div(Cplx *z, Cplx x, Cplx y) 
     147{ 
     148        z->re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im); 
     149        z->im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im); 
     150} 
     151 
     152void cplx_exp(Cplx *z, Cplx b) 
     153{ 
    175154        double br,bi; 
    176155        br=b.re; 
    177156        bi=b.im; 
    178         z.re = exp(br)*cos(bi); 
    179         z.im = exp(br)*sin(bi); 
    180         return z; 
    181 } 
    182  
    183  
    184 complex cplx_sqrt(z)    //see Schaum`s Math Handbook p. 22, 6.6 and 6.10 
    185         complex z; 
    186 { 
    187         complex c; 
     157        z->re = exp(br)*cos(bi); 
     158        z->im = exp(br)*sin(bi); 
     159} 
     160 
     161 
     162void cplx_sqrt(Cplx *c, Cplx z)    //see Schaum`s Math Handbook p. 22, 6.6 and 6.10 
     163{ 
    188164        double zr,zi,x,y,r,w; 
    189165 
     
    193169        if (zr==0.0 && zi==0.0) 
    194170        { 
    195     c.re=0.0; 
    196         c.im=0.0; 
    197     return c; 
    198         } 
    199         else 
    200         { 
     171                c->re=0.0; 
     172                c->im=0.0; 
     173        } else { 
    201174                x=fabs(zr); 
    202175                y=fabs(zi); 
     
    205178                        r=y/x; 
    206179                        w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r))); 
    207                 } 
    208                 else 
    209                 { 
     180                } else { 
    210181                        r=x/y; 
    211182                        w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r))); 
     
    213184                if (zr >=0.0) 
    214185                { 
    215                         c.re=w; 
    216                         c.im=zi/(2.0*w); 
     186                        c->re=w; 
     187                        c->im=zi/(2.0*w); 
     188                } else { 
     189                        c->im=(zi >= 0) ? w : -w; 
     190                        c->re=zi/(2.0*c->im); 
    217191                } 
    218                 else 
    219                 { 
    220                         c.im=(zi >= 0) ? w : -w; 
    221                         c.re=zi/(2.0*c.im); 
    222                 } 
    223                 return c; 
    224         } 
    225 } 
    226  
    227 complex cplx_cos(b) 
    228         complex b; 
    229 { 
    230         complex zero,two,z,i,bi,negbi; 
    231         zero = cassign(0.0,0.0); 
    232         two = cassign(2.0,0.0); 
    233         i = cassign(0.0,1.0); 
    234         bi = cplx_mult(b,i); 
    235         negbi = cplx_sub(zero,bi); 
    236         z = cplx_div(cplx_add(cplx_exp(bi),cplx_exp(negbi)),two); 
    237         return z; 
     192        } 
     193} 
     194 
     195void cplx_cos(Cplx *z, Cplx b) 
     196{ 
     197        // cos(b) = (e^bi + e^-bi)/2 
     198        //        = (e^b.im e^-i bi.re) + e^-b.im e^i b.re)/2 
     199        //        = (e^b.im cos(-b.re) + e^b.im sin(-b.re) i)/2 + (e^-b.im cos(b.re) + e^-b.im sin(b.re) i)/2 
     200        //        = e^b.im cos(b.re)/2 - e^b.im sin(b.re)/2 i + 1/e^b.im cos(b.re)/2 + 1/e^b.im sin(b.re)/2 i 
     201        //        = (e^b.im + 1/e^b.im)/2 cos(b.re) + (-e^b.im + 1/e^b.im)/2 sin(b.re) i 
     202        //        = cosh(b.im) cos(b.re) - sinh(b.im) sin(b.re) i 
     203        double exp_b_im = exp(b.im); 
     204        z->re = 0.5*(+exp_b_im + 1.0/exp_b_im) * cos(b.re); 
     205        z->im = -0.5*(exp_b_im - 1.0/exp_b_im) * sin(b.re); 
    238206} 
    239207 
  • src/sas/sascalc/calculator/c_extensions/librefl.h

    r9e531f2 r144e032a  
    55        double re; 
    66        double im; 
    7 } complex; 
     7} Cplx; 
    88 
    99typedef struct { 
    10         complex a; 
    11         complex b; 
    12         complex c; 
    13         complex d; 
     10        Cplx a; 
     11        Cplx b; 
     12        Cplx c; 
     13        Cplx d; 
    1414} matrix; 
    1515 
    16 complex cassign(double real, double imag); 
     16void cassign(Cplx*, double real, double imag); 
    1717 
    18 complex cplx_add(complex x,complex y); 
     18void cplx_add(Cplx*, Cplx x,Cplx y); 
    1919 
    20 complex rcmult(double x,complex y); 
     20void rcmult(Cplx*, double x,Cplx y); 
    2121 
    22 complex cplx_sub(complex x,complex y); 
     22void cplx_sub(Cplx*, Cplx x,Cplx y); 
    2323 
    24 complex cplx_mult(complex x,complex y); 
     24void cplx_mult(Cplx*, Cplx x,Cplx y); 
    2525 
    26 complex cplx_div(complex x,complex y); 
     26void cplx_div(Cplx*, Cplx x,Cplx y); 
    2727 
    28 complex cplx_exp(complex b); 
     28void cplx_exp(Cplx*, Cplx b); 
    2929 
    30 complex cplx_sqrt(complex z); 
     30void cplx_sqrt(Cplx*, Cplx z); 
    3131 
    32 complex cplx_cos(complex b); 
     32void cplx_cos(Cplx*, Cplx b); 
    3333 
    3434double intersldfunc(int fun_type, double n_sub, double i, double nu, double sld_l, double sld_r); 
  • src/sas/sascalc/calculator/c_extensions/sld2i.c

    ra1daf86 r144e032a  
    5454        polar_sld b_sld; 
    5555        double qr = 0.0; 
    56         complex iqr = cassign(0.0, 0.0); 
    57         complex ephase = cassign(0.0, 0.0); 
    58         complex comp_sld = cassign(0.0, 0.0); 
    59  
    60         complex sumj_uu; 
    61         complex sumj_ud; 
    62         complex sumj_du; 
    63         complex sumj_dd; 
    64         complex temp_fi; 
     56        Cplx iqr; 
     57        Cplx ephase; 
     58        Cplx comp_sld; 
     59 
     60        Cplx sumj_uu; 
     61        Cplx sumj_ud; 
     62        Cplx sumj_du; 
     63        Cplx sumj_dd; 
     64        Cplx temp_fi; 
    6565 
    6666        double count = 0.0; 
    6767        int i, j; 
     68 
     69        cassign(&iqr, 0.0, 0.0); 
     70        cassign(&ephase, 0.0, 0.0); 
     71        cassign(&comp_sld, 0.0, 0.0); 
    6872 
    6973        //Assume that pixel volumes are given in vol_pix in A^3 unit 
     
    7781        for(i=0; i<npoints; i++){ 
    7882                //I_out[i] = 0.0; 
    79                 sumj_uu = cassign(0.0, 0.0); 
    80                 sumj_ud = cassign(0.0, 0.0); 
    81                 sumj_du = cassign(0.0, 0.0); 
    82                 sumj_dd = cassign(0.0, 0.0); 
     83                cassign(&sumj_uu, 0.0, 0.0); 
     84                cassign(&sumj_ud, 0.0, 0.0); 
     85                cassign(&sumj_du, 0.0, 0.0); 
     86                cassign(&sumj_dd, 0.0, 0.0); 
    8387                //printf("i: %d\n", i); 
    8488                //q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]); 
    8589 
    8690                for(j=0; j<this->n_pix; j++){ 
    87                         //printf("j: %d\n", j); 
    8891                        if (this->sldn_val[j]!=0.0 
    8992                                ||this->mx_val[j]!=0.0 
     
    9194                                ||this->mz_val[j]!=0.0) 
    9295                        { 
     96                            // printf("i,j: %d,%d\n", i,j); 
    9397                                //anisotropic 
    94                                 temp_fi = cassign(0.0, 0.0); 
    95                                 b_sld = cal_msld(0, qx[i], qy[i], this->sldn_val[j], 
     98                                cassign(&temp_fi, 0.0, 0.0); 
     99                                cal_msld(&b_sld, 0, qx[i], qy[i], this->sldn_val[j], 
    96100                                                         this->mx_val[j], this->my_val[j], this->mz_val[j], 
    97101                                                         this->inspin, this->outspin, this->stheta); 
    98102                                qr = (qx[i]*this->x_val[j] + qy[i]*this->y_val[j]); 
    99                                 iqr = cassign(0.0, qr); 
    100                                 ephase = cplx_exp(iqr); 
     103                                cassign(&iqr, 0.0, qr); 
     104                                cplx_exp(&ephase, iqr); 
    101105 
    102106                                //Let's multiply pixel(atomic) volume here 
    103                                 ephase = rcmult(this->vol_pix[j], ephase); 
     107                                rcmult(&ephase, this->vol_pix[j], ephase); 
    104108                                //up_up 
    105109                                if (this->inspin > 0.0 && this->outspin > 0.0){ 
    106                                         comp_sld = cassign(b_sld.uu, 0.0); 
    107                                         temp_fi = cplx_mult(comp_sld, ephase); 
    108                                         sumj_uu = cplx_add(sumj_uu, temp_fi); 
     110                                        cassign(&comp_sld, b_sld.uu, 0.0); 
     111                                        cplx_mult(&temp_fi, comp_sld, ephase); 
     112                                        cplx_add(&sumj_uu, sumj_uu, temp_fi); 
    109113                                } 
    110114                                //down_down 
    111115                                if (this->inspin < 1.0 && this->outspin < 1.0){ 
    112                                         comp_sld = cassign(b_sld.dd, 0.0); 
    113                                         temp_fi = cplx_mult(comp_sld, ephase); 
    114                                         sumj_dd = cplx_add(sumj_dd, temp_fi); 
     116                                        cassign(&comp_sld, b_sld.dd, 0.0); 
     117                                        cplx_mult(&temp_fi, comp_sld, ephase); 
     118                                        cplx_add(&sumj_dd, sumj_dd, temp_fi); 
    115119                                } 
    116120                                //up_down 
    117121                                if (this->inspin > 0.0 && this->outspin < 1.0){ 
    118                                         comp_sld = cassign(b_sld.re_ud, b_sld.im_ud); 
    119                                         temp_fi = cplx_mult(comp_sld, ephase); 
    120                                         sumj_ud = cplx_add(sumj_ud, temp_fi); 
     122                                        cassign(&comp_sld, b_sld.re_ud, b_sld.im_ud); 
     123                                        cplx_mult(&temp_fi, comp_sld, ephase); 
     124                                        cplx_add(&sumj_ud, sumj_ud, temp_fi); 
    121125                                } 
    122126                                //down_up 
    123127                                if (this->inspin < 1.0 && this->outspin > 0.0){ 
    124                                         comp_sld = cassign(b_sld.re_du, b_sld.im_du); 
    125                                         temp_fi = cplx_mult(comp_sld, ephase); 
    126                                         sumj_du = cplx_add(sumj_du, temp_fi); 
    127                                 } 
    128  
     128                                        cassign(&comp_sld, b_sld.re_du, b_sld.im_du); 
     129                                        cplx_mult(&temp_fi, comp_sld, ephase); 
     130                                        cplx_add(&sumj_du, sumj_du, temp_fi); 
     131                                } 
    129132 
    130133                                if (i == 0){ 
  • src/sas/sascalc/calculator/sas_gen.py

    ra1daf86 r144e032a  
    133133        sldn = copy.deepcopy(self.data_sldn) 
    134134        sldn -= self.params['solvent_SLD'] 
    135         model = mod.new_GenI((1 if self.is_avg else 0), 
    136                              pos_x, pos_y, pos_z, 
    137                              sldn, self.data_mx, self.data_my, 
    138                              self.data_mz, self.data_vol, 
    139                              self.params['Up_frac_in'], 
    140                              self.params['Up_frac_out'], 
    141                              self.params['Up_theta']) 
     135        # **** WARNING **** new_GenI holds pointers to numpy vectors 
     136        # be sure that they are contiguous double precision arrays and make  
     137        # sure the GC doesn't eat them before genicom is called. 
     138        # TODO: rewrite so that the parameters are passed directly to genicom 
     139        args = ( 
     140            (1 if self.is_avg else 0), 
     141            pos_x, pos_y, pos_z, 
     142            sldn, self.data_mx, self.data_my, 
     143            self.data_mz, self.data_vol, 
     144            self.params['Up_frac_in'], 
     145            self.params['Up_frac_out'], 
     146            self.params['Up_theta']) 
     147        model = mod.new_GenI(*args) 
    142148        if len(qy): 
    143149            qx, qy = _vec(qx), _vec(qy) 
Note: See TracChangeset for help on using the changeset viewer.