Changeset 9624cda in sasview


Ignore:
Timestamp:
Mar 13, 2013 9:10:51 AM (12 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:
49a8843
Parents:
5f248f6
Message:

fixed debye averging in gensans

Files:
5 edited

Legend:

Unmodified
Added
Removed

  • calculatorview/src/sans/perspectives/calculator/gen_scatter_panel.py

    r644f95b r9624cda  
    2121from data_util.calcthread import CalcThread 
    2222from sans.guiframe.local_perspectives.plotting.SimplePlot import PlotFrame 
    23 from sans.guiframe.dataFitting import Data2D 
     23from sans.guiframe.dataFitting import Data2D  
     24from sans.guiframe.dataFitting import Data1D 
    2425from sans.dataloader.data_info import Detector 
    2526from sans.dataloader.data_info import Source 
    … …  
    4647_QMAX_DEFAULT = 0.3 
    4748_NPTS_DEFAULT = 50  
     49_Q1D_MIN = 0.001 
    4850 
    4951def add_icon(parent, frame): 
    … …  
    8486                 completefn = None, 
    8587                 updatefn   = None, 
    86                  elapsed = 0, 
    87                  yieldtime  = 0.005, 
     88                 #elapsed = 0, 
     89                 yieldtime  = 0.01, 
    8890                 worktime   = 0.01): 
    8991        """ 
    … …  
    9698        self.id = id  
    9799        self.input = input  
     100        self.update_fn = updatefn 
    98101         
    99102    def compute(self): 
    … …  
    101104        excuting computation 
    102105        """ 
    103         elapsed = time.time() - self.starttime 
    104         
    105         self.complete(input=self.input, 
    106                       elapsed = elapsed) 
     106        #elapsed = time.time() - self.starttime 
     107        self.starttime = time.time() 
     108        self.complete(input=self.input, update=self.update_fn) 
    107109             
    108110class SasGenPanel(ScrolledPanel, PanelBase): 
    … …  
    140142        self.data = None 
    141143        self.scale2d = None 
     144        self.is_avg = False 
    142145        self.plot_frame = None 
    143146        self.qmax_x = _QMAX_DEFAULT 
    … …  
    299302                                      style=wx.CB_READONLY)  
    300303        orient_combo.Append('Fixed orientation') 
    301         orient_combo.Append('No orientation ') 
     304        orient_combo.Append('Debye full avg.') 
     305        #orient_combo.Append('Debye sph. sym.') 
    302306         
    303307        orient_combo.Bind(wx.EVT_COMBOBOX, self._on_orient_select) 
    … …  
    311315        event.Skip() 
    312316        cb = event.GetEventObject() 
    313         is_avg = cb.GetCurrentSelection() == 1 
    314         self.model.set_is_avg(is_avg)          
     317        if cb.GetCurrentSelection() == 2: 
     318            self.is_avg = None 
     319        else: 
     320            is_avg = cb.GetCurrentSelection() == 1 
     321            self.is_avg = is_avg 
     322        self.model.set_is_avg(self.is_avg) 
     323        self.set_est_time()     
    315324            
    316325    def _layout_qrange(self): 
    … …  
    360369            Do the layout for the button widgets 
    361370        """  
    362         self.est_time = '*Estimated Computation time : %s sec' 
    363         self.time_text = wx.StaticText(self, -1, self.est_time% str(2) ) 
     371        self.est_time = '*Estimated Computation time : %s' 
     372        self.time_text = wx.StaticText(self, -1, self.est_time% str('2 sec') ) 
    364373        self.orient_combo = self._fill_orient_combo() 
    365374        self.orient_combo.Show(False) 
    … …  
    368377        self.bt_compute.SetToolTipString("Compute 2D Scattering Pattern.") 
    369378        self.button_sizer.AddMany([(self.time_text , 0, wx.LEFT, 20), 
    370                                    (self.orient_combo , 0, wx.LEFT, 20), 
     379                                   (self.orient_combo , 0, wx.LEFT, 40), 
    371380                                   (self.bt_compute, 0, wx.LEFT, 20)]) 
    372381         
    … …  
    378387        n_qbins *= n_qbins 
    379388        n_pixs = float(self.parent.get_npix()) 
     389        if self.is_avg: 
     390            n_pixs *= (n_pixs / 200) 
    380391        x_in = n_qbins * n_pixs / 100000 
    381392        # magic equation: not very accurate 
    … …  
    387398        Set text for est. computation time 
    388399        """ 
     400        unit = 'sec' 
    389401        if self.time_text != None: 
     402            self.time_text.SetForegroundColour('black') 
    390403            etime = self.estimate_ctime() 
    391             self.time_text.SetLabel(self.est_time% str(etime)) 
     404            if etime > 60: 
     405                etime /= 60 
     406                unit = 'min' 
     407                self.time_text.SetForegroundColour('red') 
     408            time_str = str(etime) + ' ' + unit 
     409            self.time_text.SetLabel(self.est_time% time_str) 
    392410         
    393411    def _do_layout(self): 
    … …  
    541559            self.file_name = filename.split('.')[0] 
    542560            self.orient_combo.SetSelection(0) 
     561            self.is_avg = False 
    543562            if self.ext in self.omfreader.ext: 
    544563                gen = sans_gen.OMF2SLD() 
    … …  
    577596        Set sld data helper 
    578597        """ 
    579         is_avg = self.orient_combo.GetCurrentSelection() == 1 
    580         self.model.set_is_avg(is_avg) 
     598        #is_avg = self.orient_combo.GetCurrentSelection() == 1 
     599        self.model.set_is_avg(self.is_avg) 
    581600        self.model.set_sld_data(self.sld_data) 
    582601         
    … …  
    741760            graph_title += " - Magnetic Vector as Arrow -"  
    742761            panel = self.plot_frame.plotpanel 
    743             def _draw_arrow(input=None, elapsed=0.1): 
     762            def _draw_arrow(input=None, update=None): 
    744763                """ 
    745764                draw magnetic vectors w/arrow 
    … …  
    828847            self.model.set_sld_data(self.sld_data) 
    829848            self.set_input_params() 
    830             self._create_default_2d_data() 
    831             i_out = numpy.zeros(len(self.data.data)) 
     849            if self.is_avg or self.is_avg == None: 
     850                self._create_default_1d_data() 
     851                i_out = numpy.zeros(len(self.data.y)) 
     852                inputs = [self.data.x, [], i_out] 
     853            else: 
     854                self._create_default_2d_data() 
     855                i_out = numpy.zeros(len(self.data.data)) 
     856                inputs=[self.data.qx_data, self.data.qy_data, i_out] 
     857                 
    832858            msg = "Computation is in progress..." 
    833859            status_type = 'progress' 
    834860            self._status_info(msg, status_type) 
    835              
    836             cal_out = CalcGen(input=[self.data.qx_data,  
    837                                      self.data.qy_data, i_out],  
     861            cal_out = CalcGen(input=inputs,  
    838862                              completefn=self.complete,  
    839863                              updatefn=self._update) 
    840             cal_out.queue() 
     864            cal_out.queue()               
    841865             
    842866        except: 
    … …  
    900924    def _update(self, time=None): 
    901925        """ 
    902         Update the output of plotting model 
    903         """ 
    904         if self.parent.parent != None: 
    905             msg = "Computation/drawing is in progress..." 
    906             wx.PostEvent(self.parent.parent,  
    907                          StatusEvent(status=msg, type="update")) 
     926        Update the progress bar 
     927        """ 
     928        if self.parent.parent == None: 
     929            return 
     930        type = "progress" 
     931        msg = "Please wait. Computing... (Note: Window may look frozen.)" 
     932        wx.PostEvent(self.parent.parent, StatusEvent(status=msg, 
     933                                                  type=type)) 
    908934                                 
    909     def complete(self, input, elapsed=None):    
     935    def complete(self, input, update=None):    
    910936        """ 
    911937        Gen compute complete function 
    912938        :Param input: input list [qx_data, qy_data, i_out] 
    913939        """ 
    914         out = self.model.runXY(input) 
    915         self._draw2D(out) 
     940        out = numpy.empty(0) 
     941        #s = time.time() 
     942        for ind in range(len(input[0])): 
     943            if self.is_avg: 
     944                if ind % 1 == 0 and update != None: 
     945                    update() 
     946                    time.sleep(0.1) 
     947                inputi = [input[0][ind:ind+1], [], input[2][ind:ind+1]] 
     948                outi = self.model.run(inputi) 
     949                out = numpy.append(out, outi) 
     950            else: 
     951                if ind % 50 == 0  and update != None: 
     952                    update() 
     953                    time.sleep(0.001) 
     954                inputi = [input[0][ind:ind+1], input[1][ind:ind+1],  
     955                          input[2][ind:ind+1]] 
     956                outi = self.model.runXY(inputi) 
     957                out = numpy.append(out, outi) 
     958        #print time.time() - s 
     959        if self.is_avg or self.is_avg == None: 
     960            self._draw1D(out) 
     961        else: 
     962            #out = self.model.runXY(input) 
     963            self._draw2D(out) 
     964             
    916965        msg = "Gen computation completed.\n" 
    917966        status_type = 'stop' 
    918967        self._status_info(msg, status_type) 
    919          
     968                
    920969    def _create_default_2d_data(self): 
    921970        """ 
    … …  
    9771026        self.data.ymin = ymin 
    9781027        self.data.ymax = ymax 
    979          
     1028 
     1029    def _create_default_1d_data(self): 
     1030        """ 
     1031        Create 2D data by default 
     1032        Only when the page is on theory mode. 
     1033        :warning: This data is never plotted. 
     1034                    residuals.x = data_copy.x[index] 
     1035            residuals.dy = numpy.ones(len(residuals.y)) 
     1036            residuals.dx = None 
     1037            residuals.dxl = None 
     1038            residuals.dxw = None 
     1039        """ 
     1040        self.qmax_x = float(self.qmax_ctl.GetValue()) 
     1041        self.npts_x = int(float(self.npt_ctl.GetValue())) 
     1042        qmax = self.qmax_x #/ numpy.sqrt(2) 
     1043        ## Default values 
     1044        xmax = qmax 
     1045        xmin = qmax * _Q1D_MIN 
     1046        qstep = self.npts_x 
     1047        x = numpy.linspace(start=xmin, stop=xmax, num=qstep, endpoint=True) 
     1048        # store x and y bin centers in q space 
     1049        #self.data.source = Source() 
     1050        y = numpy.ones(len(x)) 
     1051        dy = numpy.zeros(len(x)) 
     1052        dx = numpy.zeros(len(x)) 
     1053        self.data = Data1D(x=x, y=y) 
     1054        self.data.dx = dx 
     1055        self.data.dy = dy 
     1056 
     1057    def _draw1D(self, y_out): 
     1058        """ 
     1059        Complete get the result of modelthread and create model 2D 
     1060        that can be plot. 
     1061        """ 
     1062        page_id = self.id 
     1063        data = self.data 
     1064         
     1065        model = self.model 
     1066        state = None 
     1067         
     1068        new_plot = Data1D(x=data.x, y=y_out) 
     1069        new_plot.dx = data.dx 
     1070        new_plot.dy = data.dy 
     1071        new_plot.xaxis('\\rm{Q_{x}}', '\AA^{-1}') 
     1072        new_plot.yaxis('\\rm{Intensity}', 'cm^{-1}') 
     1073        new_plot.is_data = False 
     1074        new_plot.id = str(self.uid) + " GenData1D" 
     1075        new_plot.group_id = str(self.uid) + " Model1D" 
     1076        new_plot.name = model.name + '1d' 
     1077        new_plot.title = "Generic model1D " 
     1078        new_plot.id = str(page_id) + ': ' + self.file_name \ 
     1079                        + ' #%s'% str(self.graph_num) + "_1D" 
     1080        new_plot.group_id = str(page_id) + " Model1D"  +\ 
     1081                             ' #%s'% str(self.graph_num) + "_1D" 
     1082        new_plot.is_data = False 
     1083 
     1084        title = new_plot.title 
     1085        _yaxis, _yunit = new_plot.get_yaxis() 
     1086        _xaxis, _xunit = new_plot.get_xaxis() 
     1087        new_plot.xaxis(str(_xaxis), str(_xunit)) 
     1088        new_plot.yaxis(str(_yaxis), str(_yunit)) 
     1089         
     1090        if new_plot.is_data: 
     1091            data_name = str(new_plot.name) 
     1092        else: 
     1093            data_name = str(model.__class__.__name__) + '1d' 
     1094 
     1095        if len(title) > 1: 
     1096            new_plot.title = "Gen Theory for %s "% model.name + data_name 
     1097        new_plot.name = new_plot.id 
     1098        new_plot.label = new_plot.id 
     1099        #theory_data = deepcopy(new_plot) 
     1100        if self.parent.parent != None: 
     1101            self.parent.parent.update_theory(data_id=new_plot.id, 
     1102                                           theory=new_plot, 
     1103                                           state=state) 
     1104        title = new_plot.title 
     1105        num_graph = str(self.graph_num) 
     1106        wx.CallAfter(self.parent.draw_graph, new_plot,  
     1107                     title="Graph %s: "% num_graph + new_plot.id ) 
     1108        self.graph_num += 1 
     1109                 
    9801110    def _draw2D(self, image): 
    9811111        """ 
    … …  
    9951125        new_plot.title = "Generic model 2D " 
    9961126        new_plot.id = str(page_id) + ': ' + self.file_name \ 
    997                         + ' #%s'% str(self.graph_num) 
    998         new_plot.group_id = str(page_id) + " Model2D" 
     1127                        + ' #%s'% str(self.graph_num) + "_2D" 
     1128        new_plot.group_id = str(page_id) + " Model2D" \ 
     1129                        + ' #%s'% str(self.graph_num) + "_2D" 
    9991130        new_plot.detector = data.detector 
    10001131        new_plot.source = data.source 

  • sanscalculator/src/sans/calculator/sans_gen.py

    r11363ee r9624cda  
    3030    sld_m = factor * mag 
    3131    return sld_m 
     32 
     33def transform_center(pos_x, pos_y, pos_z): 
     34    """ 
     35    re-center  
     36     
     37    :return: posx, posy, posz   [arrays] 
     38    """  
     39    posx = pos_x - (min(pos_x) + max(pos_x)) / 2.0 
     40    posy = pos_y - (min(pos_y) + max(pos_y)) / 2.0 
     41    posz = pos_z - (min(pos_z) + max(pos_z)) / 2.0   
     42    return posx, posy, posz    
    3243 
    3344class GenSAS(BaseComponent): 
    … …  
    101112        :return: function value 
    102113        """ 
    103         len_x = len(self.data_x) 
    104         if self.is_avg: 
     114        pos_x = self.data_x 
     115        pos_y = self.data_y 
     116        pos_z = self.data_z 
     117        len_x = len(pos_x) 
     118        if self.is_avg == None: 
    105119            len_x *= -1 
     120            pos_x, pos_y, pos_z = transform_center(pos_x, pos_y, pos_z) 
    106121        len_q = len(x) 
    107122        sldn = copy.deepcopy(self.data_sldn) 
    108123        sldn -= self.params['solvent_SLD'] 
    109         model = mod.new_GenI(len_x, self.data_x, self.data_y, self.data_z,  
     124        model = mod.new_GenI(len_x, pos_x, pos_y, pos_z,  
    110125                             sldn, self.data_mx, self.data_my,  
    111126                             self.data_mz, self.data_vol, 
    … …  
    113128                             self.params['Up_frac_out'],  
    114129                             self.params['Up_theta']) 
    115  
    116         mod.genicom(model, len_q, x, y, i) 
     130        if y == []: 
     131            mod.genicom(model, len_q, x, i) 
     132        else: 
     133            mod.genicomXY(model, len_q, x, y, i) 
    117134        vol_correction = self.data_total_volume / self.params['total_volume'] 
    118135        return  self.params['scale'] * vol_correction * i + \ 
    … …  
    151168        """ 
    152169        if x.__class__.__name__ == 'list': 
     170            if len(x[1]) > 0: 
     171                msg = "Not a 1D." 
     172                raise ValueError, msg 
    153173            i_out = numpy.zeros_like(x[0]) 
    154             y_in = numpy.zeros_like(x[0]) 
     174            #import time 
     175            #s = time.time() 
    155176            # 1D I is found at y =0 in the 2D pattern 
    156             return self._gen(x[0], y_in, i_out ) 
     177            out = self._gen(x[0], [], i_out ) 
     178            #print "i_out_time", time.time() - s 
     179            return out 
    157180        else: 
    158181            msg = "Q must be given as list of qx's and qy's" 
    … …  
    168191        if x.__class__.__name__ == 'list': 
    169192            i_out = numpy.zeros_like(x[0]) 
    170             import time 
    171             s = time.time() 
     193            #import time 
     194            #s = time.time() 
    172195            out = self._gen(x[0], x[1], i_out) 
    173             print "i_out_time", time.time() - s 
     196            #print "i_out_time", time.time() - s 
    174197            return out 
    175198        else: 
    … …  
    189212        """ 
    190213        if qdist.__class__.__name__ == 'list': 
    191             out = self.runXY(qdist) 
     214            if len(qdist[1]) < 1: 
     215                out = self.run(qdist) 
     216            else: 
     217                out = self.runXY(qdist) 
    192218            return out 
    193219        else: 

  • sansmodels/src/c_gen/sld2i.cpp

    r5917637 r9624cda  
    5050 
    5151/** 
    52  * Compute 
    53  */ 
    54 void GenI :: genicom(int npoints, double *qx, double *qy, double *I_out){ 
     52 * Compute 2D anisotropic 
     53 */ 
     54void GenI :: genicomXY(int npoints, double *qx, double *qy, double *I_out){ 
    5555        //npoints is given negative for angular averaging  
    5656        // Assumes that q doesn't have qz component and sld_n is all real 
    57         double q = 0.0; 
    58         int is_avg = 0; 
     57        //double q = 0.0; 
    5958        //double Pi = 4.0*atan(1.0); 
    6059        polar_sld b_sld; 
    … …  
    6362        complex ephase = cassign(0.0, 0.0); 
    6463        complex comp_sld = cassign(0.0, 0.0); 
    65         complex sumj; 
     64 
    6665        complex sumj_uu; 
    6766        complex sumj_ud; 
    … …  
    7271        double count = 0.0; 
    7372        //check if this computation is for averaging 
    74         if (n_pix < 0 ){ 
    75                 is_avg = 1; 
    76                 n_pix = fabs(n_pix); 
    77         } 
     73 
    7874        //Assume that pixel volumes are given in vol_pix in A^3 unit 
    7975        //int x_size = 0; //in Ang 
    8076        //int y_size = 0; //in Ang 
    8177        //int z_size = 0; //in Ang 
     78         
    8279        // Loop over q-values and multiply apply matrix 
    83  
     80         
    8481        for(int i=0; i<npoints; i++){ 
    8582                //I_out[i] = 0.0; 
    86                 sumj = cassign(0.0, 0.0); 
    8783                sumj_uu = cassign(0.0, 0.0); 
    8884                sumj_ud = cassign(0.0, 0.0); 
    … …  
    9086                sumj_dd = cassign(0.0, 0.0);             
    9187                //printf ("%d ", i); 
    92                 q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]); 
     88                //q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]); 
    9389 
    9490                for(int j=0; j<n_pix; j++){ 
    95                          
    9691                        if (sldn_val[j]!=0.0||mx_val[j]!=0.0||my_val[j]!=0.0||mz_val[j]!=0.0) 
    9792                        {        
     93                                //anisotropic 
    9894                                temp_fi = cassign(0.0, 0.0); 
    9995                                b_sld = cal_msld(0, qx[i], qy[i], sldn_val[j], 
    100                                                                  mx_val[j], my_val[j], mz_val[j], 
    101                                                                  inspin, outspin, stheta); 
    102                                 if (is_avg < 1){ 
    103                                         //anisotropic 
    104                                         qr = (qx[i]*x_val[j] + qy[i]*y_val[j]); 
    105                                         iqr = cassign(0.0, qr); 
    106                                         ephase = cplx_exp(iqr); 
    107                                         } 
    108                                 else{ 
    109                                         //isotropic 
    110                                         qr = sqrt(x_val[j]*x_val[j]+y_val[j]*y_val[j]+z_val[j]*z_val[j]); 
    111                                         qr *= q; 
    112                                         if (qr == 0.0){ 
    113                                                 ephase = cassign(0.0, 1.0); 
    114                                                 } 
    115                                         else{ 
    116                                                 qr = sin(qr) / qr; 
    117                                                 ephase = cassign(qr, 0.0); 
    118                                                 } 
    119                                         } 
     96                                                         mx_val[j], my_val[j], mz_val[j], 
     97                                                         inspin, outspin, stheta); 
     98                                qr = (qx[i]*x_val[j] + qy[i]*y_val[j]); 
     99                                iqr = cassign(0.0, qr); 
     100                                ephase = cplx_exp(iqr); 
     101                                 
    120102                                //Let's multiply pixel(atomic) volume here 
    121103                                ephase = rcmult(vol_pix[j], ephase); 
    … …  
    125107                                        temp_fi = cplx_mult(comp_sld, ephase); 
    126108                                        sumj_uu = cplx_add(sumj_uu, temp_fi); 
    127                                         } 
     109                                } 
    128110                                //down_down 
    129111                                if (inspin < 1.0 && outspin < 1.0){ 
    … …  
    131113                                        temp_fi = cplx_mult(comp_sld, ephase); 
    132114                                        sumj_dd = cplx_add(sumj_dd, temp_fi); 
    133                                         } 
     115                                } 
    134116                                //up_down 
    135117                                if (inspin > 0.0 && outspin < 1.0){ 
    … …  
    137119                                        temp_fi = cplx_mult(comp_sld, ephase); 
    138120                                        sumj_ud = cplx_add(sumj_ud, temp_fi); 
    139                                         } 
     121                                } 
    140122                                //down_up 
    141123                                if (inspin < 1.0 && outspin > 0.0){ 
    … …  
    143125                                        temp_fi = cplx_mult(comp_sld, ephase); 
    144126                                        sumj_du = cplx_add(sumj_du, temp_fi); 
    145                                         } 
     127                                } 
     128 
     129 
    146130                                if (i == 0){ 
    147131                                        count += vol_pix[j]; 
    … …  
    150134                } 
    151135                //printf("aa%d=%g %g %d\n", i, (sumj_uu.re*sumj_uu.re + sumj_uu.im*sumj_uu.im), (sumj_dd.re*sumj_dd.re + sumj_dd.im*sumj_dd.im), count); 
     136 
    152137                I_out[i] = (sumj_uu.re*sumj_uu.re + sumj_uu.im*sumj_uu.im); 
    153138                I_out[i] += (sumj_ud.re*sumj_ud.re + sumj_ud.im*sumj_ud.im); 
    154139                I_out[i] += (sumj_du.re*sumj_du.re + sumj_du.im*sumj_du.im); 
    155140                I_out[i] += (sumj_dd.re*sumj_dd.re + sumj_dd.im*sumj_dd.im); 
     141 
    156142                I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix 
    157143        } 
    158144        //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]); 
    159145} 
     146/** 
     147 * Compute 1D isotropic 
     148 * Isotropic: Assumes all slds are real (no magnetic) 
     149 * Also assumes there is no polarization: No dependency on spin 
     150 */ 
     151void GenI :: genicom(int npoints, double *q, double *I_out){ 
     152        //npoints is given negative for angular averaging  
     153        // Assumes that q doesn't have qz component and sld_n is all real 
     154        //double Pi = 4.0*atan(1.0); 
     155        int is_sym = 0; 
     156        double qr = 0.0; 
     157        double sumj; 
     158        double sld_j = 0.0; 
     159        double count = 0.0; 
     160        if (n_pix < 0 ){ 
     161                is_sym = 1; 
     162                n_pix = n_pix * -1; 
     163        } 
     164        //Assume that pixel volumes are given in vol_pix in A^3 unit 
     165        // Loop over q-values and multiply apply matrix 
     166        for(int i=0; i<npoints; i++){ 
     167                sumj =0.0;               
     168                for(int j=0; j<n_pix; j++){ 
     169                        //Isotropic: Assumes all slds are real (no magnetic) 
     170                        //Also assumes there is no polarization: No dependency on spin 
     171                        if (is_sym == 1){ 
     172                                // approximation for a spherical symmetric particle 
     173                                qr = sqrt(x_val[j]*x_val[j]+y_val[j]*y_val[j]+z_val[j]*z_val[j])*q[i]; 
     174                                if (qr > 0.0){ 
     175                                        qr = sin(qr) / qr; 
     176                                        sumj += sldn_val[j] * vol_pix[j] * qr; 
     177                                } 
     178                                else{ 
     179                                        sumj += sldn_val[j] * vol_pix[j]; 
     180                                } 
     181                        } 
     182                        else{ 
     183                                //full calculation 
     184                                #pragma omp parallel for 
     185                                for(int k=0; k<n_pix; k++){ 
     186                                        sld_j =  sldn_val[j] * sldn_val[k] * vol_pix[j] * vol_pix[k]; 
     187                                        qr = (x_val[j]-x_val[k])*(x_val[j]-x_val[k])+ 
     188                                                      (y_val[j]-y_val[k])*(y_val[j]-y_val[k])+  
     189                                                      (z_val[j]-z_val[k])*(z_val[j]-z_val[k]); 
     190                                        qr = sqrt(qr) * q[i]; 
     191                                        if (qr > 0.0){ 
     192                                                sumj += sld_j*sin(qr)/qr; 
     193                                        } 
     194                                        else{ 
     195                                                sumj += sld_j; 
     196                                        } 
     197                                } 
     198                        } 
     199                        if (i == 0){ 
     200                                count += vol_pix[j]; 
     201                        } 
     202                } 
     203                I_out[i] = sumj; 
     204                if (is_sym == 1){ 
     205                        I_out[i] *= sumj; 
     206                } 
     207                I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix 
     208        } 
     209        //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]); 
     210} 

  • sansmodels/src/c_gen/sld2i.hh

    rb1174ec r9624cda  
    3939                        double s_theta); 
    4040        // compute function 
    41         virtual void genicom(int npoints, double* qx, double* qy, double *I_out); 
     41        virtual void genicomXY(int npoints, double* qx, double* qy, double *I_out); 
     42        virtual void genicom(int npoints, double* q, double *I_out); 
    4243}; 
    4344 

  • sansmodels/src/c_gen/sld2i_module.cpp

    rb1174ec r9624cda  
    7575 
    7676/** 
    77  * GenI the given input according to a given object 
     77 * GenI the given input (2D) according to a given object 
    7878 */ 
    79 PyObject * genicom_input(PyObject *, PyObject *args) { 
     79PyObject * genicom_inputXY(PyObject *, PyObject *args) { 
    8080        int npoints; 
    8181        PyObject *qx_obj; 
    … …  
    100100        GenI* s = static_cast<GenI *>(temp); 
    101101 
    102         s->genicom(npoints, qx, qy, I_out); 
     102        s->genicomXY(npoints, qx, qy, I_out); 
    103103        //return PyCObject_FromVoidPtr(s, del_genicom); 
    104104        return Py_BuildValue("i",1); 
    105105} 
    106106 
     107/** 
     108 * GenI the given 1D input according to a given object 
     109 */ 
     110PyObject * genicom_input(PyObject *, PyObject *args) { 
     111        int npoints; 
     112        PyObject *q_obj; 
     113        double *q; 
     114        PyObject *I_out_obj; 
     115        Py_ssize_t n_out; 
     116        double *I_out; 
     117        PyObject *gen_obj; 
     118 
     119        if (!PyArg_ParseTuple(args, "OiOO",  &gen_obj, &npoints, &q_obj, &I_out_obj)) return NULL; 
     120        OUTVECTOR(q_obj, q, n_out); 
     121        OUTVECTOR(I_out_obj, I_out, n_out); 
     122 
     123        // Sanity check 
     124        //if(n_in!=n_out) return Py_BuildValue("i",-1); 
     125 
     126        // Set the array pointers 
     127        void *temp = PyCObject_AsVoidPtr(gen_obj); 
     128        GenI* s = static_cast<GenI *>(temp); 
     129 
     130        s->genicom(npoints, q, I_out); 
     131        //return PyCObject_FromVoidPtr(s, del_genicom); 
     132        return Py_BuildValue("i",1); 
     133} 
    107134 
    108135/** 
    … …  
    113140                  "Create a new GenI object"}, 
    114141        {"genicom",(PyCFunction)genicom_input, METH_VARARGS, 
    115                   "genicom the given input arrays"}, 
     142                  "genicom the given 1d input arrays"}, 
     143        {"genicomXY",(PyCFunction)genicom_inputXY, METH_VARARGS, 
     144                  "genicomXY the given 2d input arrays"}, 
    116145    {NULL} 
    117146}; 
Note: See TracChangeset for help on using the changeset viewer.