source: sasview/sanscalculator/src/sans/calculator/sans_gen.py @ 318b5bbb

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 318b5bbb was 318b5bbb, checked in by Jae Cho <jhjcho@…>, 11 years ago

Added polarization and magnetic stuffs

  • Property mode set to 100644
File size: 30.9 KB
Line 
1"""
2SANS generic computation and sld file readers
3"""
4from sans.models.BaseComponent import BaseComponent
5import sans.models.sans_extension.sld2i as mod
6import numpy
7import os
8
9MFactor_AM = 2.853E-12
10MFactor_mT = 2.3164E-9
11METER2ANG = 1.0E+10
12
13def mag2sld(mag, v_unit=None):
14    """
15    Convert magnetization to magnatic SLD
16    sldm = Dm * mag where Dm = gamma * classical elec. radius/(2*Bohr magneton)
17    Dm ~ 2.853E-12 [A^(-2)] ==> Shouldn't be 2.90636E-12 [A^(-2)]???
18    """ 
19    if v_unit == "A/m":
20        factor = MFactor_AM
21    elif v_unit == "mT":
22        factor = MFactor_mT
23    else:
24        raise ValueError, "Invalid valueunit"
25    sld_m = factor * mag
26    return sld_m
27
28class GenSAS(BaseComponent):
29    """
30    Generic SAS computation Model based on sld (n & m) arrays
31    """
32    def __init__(self):
33        """
34        Init
35       
36        :Params sld_data: MagSLD object
37        """
38        # Initialize BaseComponent
39        BaseComponent.__init__(self)
40        self.sld_data = None
41        self.data_pos_unit = None
42        self.data_x = None
43        self.data_y = None
44        self.data_z = None
45        self.data_sldn = None
46        self.data_mx = None
47        self.data_my = None
48        self.data_mz = None
49        self.volume = 1000 #[A^3]
50        ## Name of the model
51        self.name = "GenSAS"
52        ## Define parameters
53        self.params = {}
54        self.params['scale']       = 1.0
55        self.params['background']  = 0.0
56        self.params['Up_frac_i']     = 1.0
57        self.params['Up_frac_f']    = 1.0
58        self.params['Up_theta']  = 0.0
59        self.description='GenSAS'
60        ## Parameter details [units, min, max]
61        self.details = {}
62        self.details['scale']      = ['', None, None]
63        self.details['background'] = ['[1/cm]', None, None]
64        self.details['Up_frac_i']    = ['[u/(u+d)]', None, None]
65        self.details['Up_frac_f']   = ['[u/(u+d)]', None, None]
66        self.details['Up_theta'] = ['[deg]', None, None]
67        # fixed parameters
68        self.fixed=[]
69       
70    def set_pixel_volume(self, volume):     
71        """
72        Set the volume of a pixel in (A^3) unit
73        """
74        self.volume = volume
75       
76    def _gen(self, x, y, i):
77        """
78        Evaluate the function
79        :Param x: x-value
80        :return: function value
81        """
82        len_x = len(self.data_x)
83        len_q = len(x)
84        model = mod.new_GenI(len_x, self.data_x, self.data_y, self.data_z, 
85                             self.data_sldn, self.data_mx, self.data_my, 
86                             self.data_mz, self.params['Up_frac_i'], 
87                             self.params['Up_frac_f'], self.params['Up_theta'])
88        mod.genicom(model, len_q, x, y, i)
89        return  self.params['scale'] * i * self.volume \
90                + self.params['background']
91       
92    def set_sld_data(self, sld_data=None):   
93        """
94        Sets sld_data
95        """
96        self.sld_data = sld_data
97        self.data_pos_unit = sld_data.pos_unit
98        self.data_x = sld_data.pos_x
99        self.data_y = sld_data.pos_y
100        self.data_z = sld_data.pos_z
101        self.data_sldn = sld_data.sld_n
102        self.data_mx = sld_data.sld_mx
103        self.data_my = sld_data.sld_my
104        self.data_mz = sld_data.sld_mz
105       
106    def getProfile(self):
107        """
108        Get SLD profile
109       
110        : return: sld_data
111        """
112        return self.sld_data
113     
114    def run(self, x = 0.0):
115        """
116        Evaluate the model
117        :param x: simple value
118        :return: (I value)
119        """
120        if x.__class__.__name__ == 'list':
121            i_out = numpy.zeros_like(x[0])
122            y_in = numpy.zero_like(x[0])
123            # 1D I is found at y =0 in the 2D pattern
124            return self._gen(x[0],y_in, i_out )
125        else:
126            msg = "Q must be given as list of qx's and qy's"
127            raise ValueError, msg
128   
129    def runXY(self, x = 0.0):
130        """
131        Evaluate the model
132        :param x: simple value
133        :return: I value
134        :Use this runXY() for the computation
135        """
136        if x.__class__.__name__ == 'list':
137            i_out = numpy.zeros_like(x[0])
138            import time
139            s=time.time()
140            out = self._gen(x[0], x[1], i_out)
141            print "i_out_time",time.time()-s
142            return out
143        else:
144            msg = "Q must be given as list of qx's and qy's"
145            raise ValueError, msg
146       
147    def evalDistribution(self, qdist):
148        """
149        Evaluate a distribution of q-values.
150        * For 1D, a numpy array is expected as input:
151            evalDistribution(q)   
152        where q is a numpy array.
153        * For 2D, a list of numpy arrays are expected: [qx_prime,qy_prime],
154          where 1D arrays, 
155        :param qdist: ndarray of scalar q-values or list [qx,qy]
156                    where qx,qy are 1D ndarrays
157        """
158        if qdist.__class__.__name__ == 'list':
159            out = self.runXY(qdist)
160            return out
161        else:
162            mesg = "evalDistribution is expecting an ndarray of "
163            mesg += "a list [qx,qy] where qx,qy are arrays."
164            raise RuntimeError, mesg
165     
166class OMF2SLD:
167    """
168    Convert OMFData to MAgData
169    """
170    def __init__(self): 
171        """
172        Init
173        """
174        self.pos_x = None
175        self.pos_y = None
176        self.pos_z = None
177        self.mx = None
178        self.my = None
179        self.mz = None
180        self.sld_n = None
181        self.output = None
182        self.omfdata = None
183   
184    def set_data(self, omfdata, shape='rectangular'):
185        """
186        Set all data
187        """
188        self.omfdata = omfdata
189        length = int(omfdata.xnodes * omfdata.ynodes * omfdata.znodes)
190        pos_x = numpy.arange(omfdata.xmin, 
191                             omfdata.xnodes*omfdata.xstepsize + omfdata.xmin, 
192                             omfdata.xstepsize)
193        pos_y = numpy.arange(omfdata.ymin, 
194                             omfdata.ynodes*omfdata.ystepsize + omfdata.ymin, 
195                             omfdata.ystepsize)
196        pos_z = numpy.arange(omfdata.zmin, 
197                             omfdata.znodes*omfdata.zstepsize + omfdata.zmin, 
198                             omfdata.zstepsize)
199        self.pos_x = numpy.tile(pos_x, int(omfdata.ynodes * omfdata.znodes))
200        self.pos_y = pos_y.repeat(int(omfdata.xnodes))
201        self.pos_y = numpy.tile(self.pos_y, int(omfdata.znodes))
202        self.pos_z = pos_z.repeat(int(omfdata.xnodes * omfdata.ynodes))
203        self.mx = omfdata.mx
204        self.my = omfdata.my
205        self.mz = omfdata.mz
206        self.sld_n = numpy.zeros(length)
207           
208        if omfdata.mx == None:
209            self.mx = numpy.zeros(length)
210        if omfdata.my == None:
211            self.my = numpy.zeros(length)
212        if omfdata.mz == None:
213            self.mz = numpy.zeros(length)
214       
215        self._check_data_length(length)
216        self.remove_null_points()
217        mask = numpy.ones(len(self.sld_n), dtype=bool)
218        if shape == 'ellipsoid':
219            try:
220                # Pixel (step) size included
221                x_r = (max(self.pos_x) - min(self.pos_x) + \
222                       omfdata.xstepsize) / 2.0
223                y_r = (max(self.pos_y) - min(self.pos_y) + \
224                       omfdata.ystepsize) / 2.0
225                z_r = (max(self.pos_z) - min(self.pos_z) + \
226                       omfdata.zstepsize) / 2.0
227                x_dir2 = (self.pos_x / x_r) * (self.pos_x / x_r)
228                y_dir2 = (self.pos_y / y_r) * (self.pos_y / y_r)
229                z_dir2 = (self.pos_z / z_r) * (self.pos_z / z_r)
230                mask = (x_dir2 + y_dir2 + z_dir2) <= 1
231            except:
232                pass
233        self.output = MagSLD(self.pos_x[mask], self.pos_y[mask], 
234                             self.pos_z[mask], self.sld_n[mask], 
235                             self.mx[mask], self.my[mask], self.mz[mask])
236       
237    def get_omfdata(self):
238        """
239        Return all data
240        """
241        return self.omfdata
242   
243    def get_output(self):
244        """
245        Return output
246        """
247        return self.output
248   
249    def _check_data_length(self, length):
250        """
251        Check if the data lengths are consistent
252        :Params length: data length
253        """
254        msg = "Error: Inconsistent data length."
255        if len(self.pos_x) != length:
256            raise ValueError, msg
257        if len(self.pos_y) != length:
258            raise ValueError, msg
259        if len(self.pos_z) != length:
260            raise ValueError, msg
261        if len(self.mx) != length:
262            raise ValueError, msg
263        if len(self.my) != length:
264            raise ValueError, msg
265        if len(self.mz) != length:
266            raise ValueError, msg
267       
268    def remove_null_points(self, remove=False):
269        """
270        Removes any mx, my, and mz = 0 points
271        """
272        if remove:
273            is_nonzero = (numpy.fabs(self.mx) + numpy.fabs(self.my) + 
274                          numpy.fabs(self.mz)).nonzero() 
275            if len(is_nonzero[0]) > 0: 
276                self.pos_x = self.pos_x[is_nonzero]
277                self.pos_y = self.pos_y[is_nonzero]
278                self.pos_z = self.pos_z[is_nonzero]
279                self.sld_n = self.sld_n[is_nonzero]
280                self.mx = self.mx[is_nonzero]
281                self.my = self.my[is_nonzero]
282                self.mz = self.mz[is_nonzero]
283        self.pos_x -= (min(self.pos_x) + max(self.pos_x)) / 2.0
284        self.pos_y -= (min(self.pos_y) + max(self.pos_y)) / 2.0
285        self.pos_z -= (min(self.pos_z) + max(self.pos_z)) / 2.0
286 
287    def get_magsld(self):
288        """
289        return MagSLD
290        """
291        return self.output
292   
293       
294class OMFReader:
295    """
296    Class to load omf/ascii files (3 columns w/header).
297    """
298    ## File type
299    type_name = "OMF ASCII"
300   
301    ## Wildcards
302    type = ["OMF files (*.OMF, *.omf)|*.omf"]
303    ## List of allowed extensions
304    ext = ['.omf', '.OMF']
305       
306    def read(self, path):
307        """
308        Load data file
309        :param path: file path
310        :return: x, y, z, sld_n, sld_mx, sld_my, sld_mz
311        """
312        data_conv_m = None
313        desc = ""
314        mx = numpy.zeros(0)
315        my = numpy.zeros(0)
316        mz = numpy.zeros(0)
317        try:
318            data_conv_val = None
319            data_conv_mesh = None
320            input_f = open(path, 'rb')
321            buff = input_f.read()
322            lines = buff.split('\n')
323            input_f.close()
324            output = OMFData()
325            valueunit = None
326            for line in lines:
327                toks = line.split()
328                # Read data
329                try:
330                    _mx = float(toks[0])
331                    _my = float(toks[1])
332                    _mz = float(toks[2])
333                    _mx = mag2sld(_mx, valueunit)
334                    _my = mag2sld(_my, valueunit)
335                    _mz = mag2sld(_mz, valueunit)
336                    mx = numpy.append(mx, _mx)
337                    my = numpy.append(my, _my)
338                    mz = numpy.append(mz, _mz)
339                except:
340                    # Skip non-data lines
341                    pass
342                #Reading Header; Segment count ignored
343                s_line = line.split(":", 1)
344                if s_line[0].lower().count("oommf") > 0:
345                    oommf = s_line[1].lstrip()
346                if s_line[0].lower().count("title") > 0:
347                    title = s_line[1].lstrip()
348                if s_line[0].lower().count("desc") > 0:
349                    desc += s_line[1].lstrip()
350                    desc += '\n'
351                if s_line[0].lower().count("meshtype") > 0:
352                    meshtype = s_line[1].lstrip()
353                if s_line[0].lower().count("meshunit") > 0:
354                    meshunit = s_line[1].lstrip()
355                    if meshunit.count("m") < 1:
356                        msg = "Error: \n"
357                        msg += "We accept only m as meshunit"
358                        raise ValueError, msg
359                if s_line[0].lower().count("xbase") > 0:
360                    xbase = s_line[1].lstrip()
361                if s_line[0].lower().count("ybase") > 0:
362                    ybase = s_line[1].lstrip()
363                if s_line[0].lower().count("zbase") > 0:
364                    zbase = s_line[1].lstrip()
365                if s_line[0].lower().count("xstepsize") > 0:
366                    xstepsize = s_line[1].lstrip()
367                if s_line[0].lower().count("ystepsize") > 0:
368                    ystepsize = s_line[1].lstrip()
369                if s_line[0].lower().count("zstepsize") > 0:
370                    zstepsize = s_line[1].lstrip()
371                if s_line[0].lower().count("xnodes") > 0:
372                    xnodes = s_line[1].lstrip()
373                if s_line[0].lower().count("ynodes") > 0:
374                    ynodes = s_line[1].lstrip()
375                if s_line[0].lower().count("znodes") > 0:
376                    znodes = s_line[1].lstrip()
377                if s_line[0].lower().count("xmin") > 0:
378                    xmin = s_line[1].lstrip()
379                if s_line[0].lower().count("ymin") > 0:
380                    ymin = s_line[1].lstrip()
381                if s_line[0].lower().count("zmin") > 0:
382                    zmin = s_line[1].lstrip()
383                if s_line[0].lower().count("xmax") > 0:
384                    xmax = s_line[1].lstrip()
385                if s_line[0].lower().count("ymax") > 0:
386                    ymax = s_line[1].lstrip()
387                if s_line[0].lower().count("zmax") > 0:
388                    zmax = s_line[1].lstrip()
389                if s_line[0].lower().count("valueunit") > 0:
390                    valueunit = s_line[1].lstrip().rstrip()
391                if s_line[0].lower().count("valuemultiplier") > 0:
392                    valuemultiplier = s_line[1].lstrip()
393                if s_line[0].lower().count("valuerangeminmag") > 0:
394                    valuerangeminmag = s_line[1].lstrip()
395                if s_line[0].lower().count("valuerangemaxmag") > 0:
396                    valuerangemaxmag = s_line[1].lstrip()
397                if s_line[0].lower().count("end") > 0:
398                    #output.set_sldms(mx, my, mz)   
399                    output.filename = os.path.basename(path)
400                    output.oommf = oommf
401                    output.title = title
402                    output.desc = desc
403                    output.meshtype = meshtype
404                    output.xbase = float(xbase) * METER2ANG
405                    output.ybase = float(ybase) * METER2ANG
406                    output.zbase = float(zbase) * METER2ANG
407                    output.xstepsize = float(xstepsize) * METER2ANG
408                    output.ystepsize = float(ystepsize) * METER2ANG
409                    output.zstepsize = float(zstepsize) * METER2ANG
410                    output.xnodes = float(xnodes)
411                    output.ynodes = float(ynodes)
412                    output.znodes = float(znodes)
413                    output.xmin = float(xmin) * METER2ANG
414                    output.ymin = float(ymin) * METER2ANG
415                    output.zmin = float(zmin) * METER2ANG
416                    output.xmax = float(xmax) * METER2ANG
417                    output.ymax = float(ymax) * METER2ANG
418                    output.zmax = float(zmax) * METER2ANG
419                    output.valuemultiplier = valuemultiplier
420                    output.valuerangeminmag = mag2sld(float(valuerangeminmag),\
421                                                      valueunit)
422                    output.valuerangemaxmag = mag2sld(float(valuerangemaxmag),\
423                                                      valueunit)
424            output.set_m(mx, my, mz)
425            return output
426        except:
427            msg = "%s is not supported: \n" % path
428            msg += "We accept only Text format OMF file."
429            raise RuntimeError, msg
430
431class SLDReader:
432    """
433    Class to load ascii files (7 columns).
434    """
435    ## File type
436    type_name = "SLD ASCII"
437   
438    ## Wildcards
439    type = ["sld files (*.SLD, *.sld)|*.sld",
440            "txt files (*.TXT, *.txt)|*.txt"]
441    ## List of allowed extensions
442    ext = ['.txt', '.TXT', '.sld', 'SLD', '.*']
443   
444    def read(self, path):
445        """
446        Load data file
447       
448        :param path: file path
449       
450        :return: x, y, z, sld_n, sld_mx, sld_my, sld_mz
451       
452        :raise RuntimeError: when the file can't be opened
453        :raise ValueError: when the length of the data vectors are inconsistent
454        """
455        pos_x = numpy.zeros(0)
456        pos_y = numpy.zeros(0)
457        pos_z = numpy.zeros(0)
458        sld_n = numpy.zeros(0)
459        sld_mx = numpy.zeros(0)
460        sld_my = numpy.zeros(0)
461        sld_mz = numpy.zeros(0)
462        try:
463            input_f = open(path, 'rb')
464            buff = input_f.read()
465            lines = buff.split('\n')
466            input_f.close()
467            for line in lines:
468                toks = line.split()
469                try:
470                    _pos_x = float(toks[0])
471                    _pos_y = float(toks[1])
472                    _pos_z = float(toks[2])
473                    _sld_n = float(toks[3])
474                    _sld_mx = float(toks[4])
475                    _sld_my = float(toks[5])
476                    _sld_mz = float(toks[6])
477
478                    pos_x = numpy.append(pos_x, _pos_x)
479                    pos_y = numpy.append(pos_y, _pos_y)
480                    pos_z = numpy.append(pos_z, _pos_z)
481                    sld_n = numpy.append(sld_n, _sld_n)
482                    sld_mx = numpy.append(sld_mx, _sld_mx)
483                    sld_my = numpy.append(sld_my, _sld_my)
484                    sld_mz = numpy.append(sld_mz, _sld_mz)
485                except:
486                    # Skip non-data lines
487                    pass
488            pos_x -= (min(pos_x) + max(pos_x)) / 2.0
489            pos_y -= (min(pos_y) + max(pos_y)) / 2.0
490            pos_z -= (min(pos_z) + max(pos_z)) / 2.0
491
492            output = MagSLD(pos_x, pos_y, pos_z, sld_n, sld_mx, sld_my, sld_mz)
493            output.filename = os.path.basename(path)
494            return output
495        except:
496            RuntimeError, "%s is not a sld file" % path
497       
498    def write(self, path, data):
499        """
500        Write sld file
501       
502        :Param path: file path
503        :Param data: MagSLD data object
504        """
505        if path == None:
506            raise ValueError, "Missing the file path."
507        if data == None:
508            raise ValueError, "Missing the data to save."
509       
510        x_val = data.pos_x
511        y_val = data.pos_y
512        z_val = data.pos_z
513       
514        length = len(x_val)
515       
516        sld_n = data.sld_n
517        if sld_n == None:
518            sld_n = numpy.zerros(length)
519           
520        sld_mx = data.sld_mx
521        if sld_mx == None:
522            sld_mx = numpy.zerros(length)
523            sld_my = numpy.zerros(length)
524            sld_mz = numpy.zerros(length)
525        else:
526            sld_my = data.sld_my
527            sld_mz = data.sld_mz
528           
529        out = open(path, 'w')
530        # First Line: Column names
531        out.write("X  Y  Z  SLDN  SLDMx  SLDMy  SLDMz")
532        for ind in range(length):
533            out.write("\n%g  %g  %g  %g  %g  %g  %g" % (x_val[ind], y_val[ind], 
534                                                    z_val[ind], sld_n[ind], 
535                                                    sld_mx[ind], sld_my[ind], 
536                                                    sld_mz[ind])) 
537        out.close()       
538           
539           
540class OMFData:
541    """
542    OMF Data.
543    """
544    _meshunit = "A"
545    _valueunit = "A^(-2)"
546    def __init__(self):
547        """
548        Init for mag SLD
549        """
550        self.filename = 'default'
551        self.oommf = ''
552        self.title = ''
553        self.desc = ''
554        self.meshtype = ''
555        self.meshunit = self._meshunit
556        self.valueunit = self._valueunit
557        self.xbase = 0.0
558        self.ybase = 0.0
559        self.zbase = 0.0
560        self.xstepsize = 6.0
561        self.ystepsize = 6.0
562        self.zstepsize = 6.0
563        self.xnodes = 10.0
564        self.ynodes = 10.0
565        self.znodes = 10.0
566        self.xmin = 0.0
567        self.ymin = 0.0
568        self.zmin = 0.0
569        self.xmax = 60.0
570        self.ymax = 60.0
571        self.zmax = 60.0
572        self.mx = None
573        self.my = None
574        self.mz = None
575        self.valuemultiplier = 1.
576        self.valuerangeminmag = 0
577        self.valuerangemaxmag = 0
578       
579    def __str__(self):
580        """
581        doc strings
582        """
583        _str =  "Type:            %s\n" % self.__class__.__name__
584        _str += "File:            %s\n" % self.filename
585        _str += "OOMMF:           %s\n" % self.oommf
586        _str += "Title:           %s\n" % self.title
587        _str += "Desc:            %s\n" % self.desc
588        _str += "meshtype:        %s\n" % self.meshtype
589        _str += "meshunit:        %s\n" % str(self.meshunit)
590        _str += "xbase:           %s [%s]\n" % (str(self.xbase), self.meshunit)
591        _str += "ybase:           %s [%s]\n" % (str(self.ybase), self.meshunit)
592        _str += "zbase:           %s [%s]\n" % (str(self.zbase), self.meshunit)
593        _str += "xstepsize:       %s [%s]\n" % (str(self.xstepsize), 
594                                                self.meshunit)
595        _str += "ystepsize:       %s [%s]\n" % (str(self.ystepsize), 
596                                                self.meshunit)
597        _str += "zstepsize:       %s [%s]\n" % (str(self.zstepsize), 
598                                                self.meshunit)
599        _str += "xnodes:          %s\n" % str(self.xnodes)
600        _str += "ynodes:          %s\n" % str(self.ynodes)
601        _str += "znodes:          %s\n" % str(self.znodes)
602        _str += "xmin:            %s [%s]\n" % (str(self.xmin), self.meshunit)
603        _str += "ymin:            %s [%s]\n" % (str(self.ymin), self.meshunit)
604        _str += "zmin:            %s [%s]\n" % (str(self.zmin), self.meshunit)
605        _str += "xmax:            %s [%s]\n" % (str(self.xmax), self.meshunit)
606        _str += "ymax:            %s [%s]\n" % (str(self.ymax), self.meshunit)
607        _str += "zmax:            %s [%s]\n" % (str(self.zmax), self.meshunit)
608        _str += "valueunit:       %s\n" % self.valueunit
609        _str += "valuemultiplier: %s\n" % str(self.valuemultiplier)
610        _str += "ValueRangeMinMag:%s [%s]\n" % (str(self.valuerangeminmag), 
611                                                 self.valueunit)
612        _str += "ValueRangeMaxMag:%s [%s]\n" % (str(self.valuerangemaxmag), 
613                                                 self.valueunit)
614        return _str
615   
616    def set_m(self, mx, my, mz):
617        """
618        Set the Mx, My, Mz values
619        """
620        self.mx = mx
621        self.my = my
622        self.mz = mz
623               
624class MagSLD:
625    """
626    Magnetic SLD.
627    """
628    pos_x = None
629    pos_y = None
630    pos_z = None
631    sld_n = None
632    sld_mx = None
633    sld_my = None
634    sld_mz = None
635    # Units
636    _pos_unit = 'A'
637    _sld_unit = '1/A^(2)'
638   
639    def __init__(self, pos_x, pos_y, pos_z, sld_n=None, 
640                 sld_mx=None, sld_my=None, sld_mz=None):
641        """
642        Init for mag SLD
643        :params : All should be numpy 1D array
644        """
645        self.is_data = True
646        self.filename = ''
647        self.xstepsize = 6.0
648        self.ystepsize = 6.0
649        self.zstepsize = 6.0
650        self.xnodes = 10.0
651        self.ynodes = 10.0
652        self.znodes = 10.0
653        self.has_stepsize = False
654        self.pos_unit = self._pos_unit
655        self.sld_unit = self._sld_unit
656        self.pos_x = pos_x
657        self.pos_y = pos_y
658        self.pos_z = pos_z
659        self.sld_n = sld_n
660        self.sld_mx = sld_mx
661        self.sld_my = sld_my
662        self.sld_mz = sld_mz
663        self.sld_m = None
664        self.sld_phi = None
665        self.sld_theta = None
666        if sld_mx != None and sld_my != None and sld_mz != None:
667            self.set_sldms(sld_mx, sld_my, sld_mz)
668        self.set_nodes()
669                   
670    def __str__(self):
671        """
672        doc strings
673        """
674        _str =  "Type:       %s\n" % self.__class__.__name__
675        _str += "File:       %s\n" % self.filename
676        _str += "Axis_unit:  %s\n" % self.pos_unit
677        _str += "SLD_unit:   %s\n" % self.sld_unit
678        return _str
679   
680    def set_sldn(self, sld_n):
681        """
682        Sets neutron SLD
683        """
684        if sld_n.__class__.__name__ == 'float':
685            if self.is_data:
686                # For data, put the value to only the pixels w non-zero M
687                self.sld_n = numpy.zeros(len(self.pos_x))
688                is_nonzero = (numpy.fabs(self.sld_mx) + 
689                              numpy.fabs(self.sld_my) + 
690                              numpy.fabs(self.sld_mz)).nonzero() 
691                self.sld_n[is_nonzero] = sld_n
692            else:
693                # For non-data, put the value to all the pixels
694                self.sld_n = numpy.ones(len(self.pos_x)) * sld_n
695        else:
696            self.sld_n = sld_n
697           
698    def set_sldms(self, sld_mx, sld_my, sld_mz):
699        """
700        Sets (|m|, m_theta, m_phi)
701        """
702        if sld_mx.__class__.__name__ == 'float':
703            self.sld_mx = numpy.ones(len(self.pos_x)) * sld_mx
704        else:
705            self.sld_mx = sld_mx
706        if sld_my.__class__.__name__ == 'float':
707            self.sld_my = numpy.ones(len(self.pos_x)) * sld_my
708        else:
709            self.sld_my = sld_my
710        if sld_mz.__class__.__name__ == 'float':
711            self.sld_mz = numpy.ones(len(self.pos_x)) * sld_mz
712        else:
713            self.sld_mz = sld_mz
714
715        sld_m = numpy.sqrt(sld_mx * sld_mx + sld_my * sld_my + \
716                                sld_mz * sld_mz) 
717        self.sld_m = sld_m
718       
719    def get_sldn(self):
720        """
721        Returns nuclear sld
722        """
723        return self.sld_n
724   
725    def get_sld_mxyz(self):
726        """
727        Returns (sld_mx, sldmy, sld_mz)
728        """
729        return (self.sld_mx, self.sldmy, self.sld_mz)
730   
731    def get_sld_m(self):
732        """
733        Returns (sldm, sld_theta, sld_phi)
734        """
735        return (self.sldm, self.sld_theta, self.sld_phi)
736   
737    def set_nodes(self):
738        """
739        Set xnodes, ynodes, and znodes
740        """
741        if not self.has_stepsize:
742            self.set_stepsize()
743        try:
744            xdist = (max(self.pos_x) - min(self.pos_x)) / self.xstepsize
745            ydist = (max(self.pos_y) - min(self.pos_y)) / self.ystepsize
746            zdist = (max(self.pos_z) - min(self.pos_z)) / self.zstepsize
747            self.xnodes = int(xdist) + 1
748            self.ynodes = int(ydist) + 1
749            self.znodes = int(zdist) + 1
750        except:
751            self.xnodes = None
752            self.ynodes = None
753            self.znodes = None
754       
755    def set_stepsize(self):
756        """
757        Set xtepsize, ystepsize, and zstepsize
758        """
759        try:
760            xpos_pre = self.pos_x[0]
761            ypos_pre = self.pos_y[0]
762            zpos_pre = self.pos_z[0]
763            for x_pos in self.pos_x:
764                if xpos_pre != x_pos:
765                    self.xstepsize = numpy.fabs(x_pos - xpos_pre)
766                    break
767            for y_pos in self.pos_y:
768                if ypos_pre != y_pos:
769                    self.ystepsize = numpy.fabs(y_pos - ypos_pre)
770                    break
771            for z_pos in self.pos_z:
772                if zpos_pre != z_pos:
773                    self.zstepsize = numpy.fabs(z_pos - zpos_pre)
774                    break
775            self.has_stepsize = True
776        except:
777            self.xstepsize = None
778            self.ystepsize = None
779            self.zstepsize = None
780            self.has_stepsize = False
781        return self.xstepsize, self.ystepsize, self.zstepsize
782
783   
784def test_load():
785    from danse.common.plottools.arrow3d import Arrow3D
786    import os
787    dir = os.path.abspath(os.path.curdir)
788    print dir
789    for i in range(6):
790        dir, _ = os.path.split(dir)
791        tfile = os.path.join(dir, "test", "CoreXY_ShellZ.txt")
792        ofile = os.path.join(dir, "test", "A_Raw_Example-1.omf")
793        if os.path.isfile(tfile):
794            tfpath = tfile
795            ofpath = ofile
796            break
797    reader = SLDReader()
798    oreader = OMFReader()
799    output = reader.read(tfpath)
800    ooutput = oreader.read(ofpath)
801    foutput = OMF2SLD()
802    foutput.set_data(ooutput)
803
804    import matplotlib.pyplot as plt
805    from mpl_toolkits.mplot3d import Axes3D
806    fig = plt.figure()
807    ax = fig.gca(projection='3d')
808    ax.plot(output.pos_x, output.pos_y, output.pos_z, '.', c="g", 
809            alpha = 0.7, markeredgecolor='gray',rasterized=True)
810    gap = 7
811    max_mx = max(output.sld_mx)
812    max_my = max(output.sld_my)
813    max_mz = max(output.sld_mz)
814    max_m = max(max_mx,max_my,max_mz)
815    x2 = output.pos_x+output.sld_mx/max_m * gap
816    y2 = output.pos_y+output.sld_my/max_m * gap
817    z2 = output.pos_z+output.sld_mz/max_m * gap
818    x_arrow = numpy.column_stack((output.pos_x,x2))
819    y_arrow = numpy.column_stack((output.pos_y,y2))
820    z_arrow = numpy.column_stack((output.pos_z,z2))
821    unit_x2 = output.sld_mx / max_m
822    unit_y2 = output.sld_my / max_m
823    unit_z2 = output.sld_mz / max_m
824    color_x = numpy.fabs(unit_x2 * 0.8)
825    color_y = numpy.fabs(unit_y2 * 0.8)
826    color_z = numpy.fabs(unit_z2 * 0.8)
827    colors =  numpy.column_stack((color_x, color_y, color_z))
828    a = Arrow3D(None, x_arrow,y_arrow,z_arrow, colors, mutation_scale=10, lw=1, 
829                arrowstyle="->", color="y", alpha = 0.5)
830
831    ax.add_artist(a)
832    plt.show()
833   
834def test():
835    import os
836    dir = os.path.abspath(os.path.curdir)
837    for i in range(3):
838        dir, _ = os.path.split(dir)
839        #tfile = os.path.join(dir, "test", "C_Example_Converted.txt")
840        ofile = os.path.join(dir, "test", "A_Raw_Example-1.omf")
841        if os.path.isfile(ofile):
842            #tfpath = tfile
843            ofpath = ofile
844            break
845    #reader = SLDReader()
846    oreader = OMFReader()
847    #output = reader.read(tfpath)
848    ooutput = oreader.read(ofpath)
849    foutput = OMF2SLD() 
850    foutput.set_data(ooutput) 
851    writer =  SLDReader() 
852    writer.write(os.path.join(os.path.dirname(ofpath), "out.txt"), 
853                 foutput.output)
854    model = GenSAS()
855    model.set_sld_data(foutput.output) 
856    x = numpy.arange(1000)/10000. + 1e-5
857    y = numpy.arange(1000)/10000. + 1e-5
858    z = numpy.arange(1000)/10000. + 1e-5
859    i = numpy.zeros(1000)
860    out = model.runXY([x,y,i])
861       
862if __name__ == "__main__": 
863    #test()
864    test_load()
Note: See TracBrowser for help on using the repository browser.