[b6627d9] | 1 | # pylint: disable=invalid-name |
---|
[51f14603] | 2 | """ |
---|
[fd5ac0d] | 3 | SAS generic computation and sld file readers |
---|
[51f14603] | 4 | """ |
---|
[a1b8fee] | 5 | from __future__ import print_function |
---|
| 6 | |
---|
[51f14603] | 7 | import os |
---|
[b6627d9] | 8 | import sys |
---|
[1d014cb] | 9 | import copy |
---|
[b6627d9] | 10 | import logging |
---|
[51f14603] | 11 | |
---|
[1d014cb] | 12 | from periodictable import formula |
---|
| 13 | from periodictable import nsf |
---|
| 14 | import numpy as np |
---|
| 15 | |
---|
| 16 | from .core import sld2i as mod |
---|
| 17 | from .BaseComponent import BaseComponent |
---|
| 18 | |
---|
[463e7ffc] | 19 | logger = logging.getLogger(__name__) |
---|
[c155a16] | 20 | |
---|
[1d014cb] | 21 | if sys.version_info[0] < 3: |
---|
| 22 | def decode(s): |
---|
| 23 | return s |
---|
| 24 | else: |
---|
| 25 | def decode(s): |
---|
| 26 | return s.decode() if isinstance(s, bytes) else s |
---|
| 27 | |
---|
[e3f77d8b] | 28 | MFACTOR_AM = 2.853E-12 |
---|
| 29 | MFACTOR_MT = 2.3164E-9 |
---|
[51f14603] | 30 | METER2ANG = 1.0E+10 |
---|
| 31 | #Avogadro constant [1/mol] |
---|
| 32 | NA = 6.02214129e+23 |
---|
| 33 | |
---|
| 34 | def mag2sld(mag, v_unit=None): |
---|
| 35 | """ |
---|
| 36 | Convert magnetization to magnatic SLD |
---|
| 37 | sldm = Dm * mag where Dm = gamma * classical elec. radius/(2*Bohr magneton) |
---|
| 38 | Dm ~ 2.853E-12 [A^(-2)] ==> Shouldn't be 2.90636E-12 [A^(-2)]??? |
---|
[e3f77d8b] | 39 | """ |
---|
[51f14603] | 40 | if v_unit == "A/m": |
---|
[e3f77d8b] | 41 | factor = MFACTOR_AM |
---|
[51f14603] | 42 | elif v_unit == "mT": |
---|
[e3f77d8b] | 43 | factor = MFACTOR_MT |
---|
[51f14603] | 44 | else: |
---|
[574adc7] | 45 | raise ValueError("Invalid valueunit") |
---|
[51f14603] | 46 | sld_m = factor * mag |
---|
| 47 | return sld_m |
---|
| 48 | |
---|
| 49 | def transform_center(pos_x, pos_y, pos_z): |
---|
| 50 | """ |
---|
[e3f77d8b] | 51 | re-center |
---|
[51f14603] | 52 | :return: posx, posy, posz [arrays] |
---|
[e3f77d8b] | 53 | """ |
---|
[51f14603] | 54 | posx = pos_x - (min(pos_x) + max(pos_x)) / 2.0 |
---|
| 55 | posy = pos_y - (min(pos_y) + max(pos_y)) / 2.0 |
---|
[e3f77d8b] | 56 | posz = pos_z - (min(pos_z) + max(pos_z)) / 2.0 |
---|
| 57 | return posx, posy, posz |
---|
[51f14603] | 58 | |
---|
| 59 | class GenSAS(BaseComponent): |
---|
| 60 | """ |
---|
| 61 | Generic SAS computation Model based on sld (n & m) arrays |
---|
| 62 | """ |
---|
| 63 | def __init__(self): |
---|
| 64 | """ |
---|
| 65 | Init |
---|
| 66 | :Params sld_data: MagSLD object |
---|
| 67 | """ |
---|
| 68 | # Initialize BaseComponent |
---|
| 69 | BaseComponent.__init__(self) |
---|
| 70 | self.sld_data = None |
---|
| 71 | self.data_pos_unit = None |
---|
| 72 | self.data_x = None |
---|
| 73 | self.data_y = None |
---|
| 74 | self.data_z = None |
---|
| 75 | self.data_sldn = None |
---|
| 76 | self.data_mx = None |
---|
| 77 | self.data_my = None |
---|
| 78 | self.data_mz = None |
---|
| 79 | self.data_vol = None #[A^3] |
---|
| 80 | self.is_avg = False |
---|
| 81 | ## Name of the model |
---|
| 82 | self.name = "GenSAS" |
---|
| 83 | ## Define parameters |
---|
| 84 | self.params = {} |
---|
[e3f77d8b] | 85 | self.params['scale'] = 1.0 |
---|
| 86 | self.params['background'] = 0.0 |
---|
| 87 | self.params['solvent_SLD'] = 0.0 |
---|
[51f14603] | 88 | self.params['total_volume'] = 1.0 |
---|
[e3f77d8b] | 89 | self.params['Up_frac_in'] = 1.0 |
---|
| 90 | self.params['Up_frac_out'] = 1.0 |
---|
| 91 | self.params['Up_theta'] = 0.0 |
---|
[51f14603] | 92 | self.description = 'GenSAS' |
---|
| 93 | ## Parameter details [units, min, max] |
---|
| 94 | self.details = {} |
---|
[9a5097c] | 95 | self.details['scale'] = ['', 0.0, np.inf] |
---|
| 96 | self.details['background'] = ['[1/cm]', 0.0, np.inf] |
---|
| 97 | self.details['solvent_SLD'] = ['1/A^(2)', -np.inf, np.inf] |
---|
| 98 | self.details['total_volume'] = ['A^(3)', 0.0, np.inf] |
---|
[cb4ef58] | 99 | self.details['Up_frac_in'] = ['[u/(u+d)]', 0.0, 1.0] |
---|
| 100 | self.details['Up_frac_out'] = ['[u/(u+d)]', 0.0, 1.0] |
---|
[9a5097c] | 101 | self.details['Up_theta'] = ['[deg]', -np.inf, np.inf] |
---|
[51f14603] | 102 | # fixed parameters |
---|
| 103 | self.fixed = [] |
---|
[e3f77d8b] | 104 | |
---|
| 105 | def set_pixel_volumes(self, volume): |
---|
[51f14603] | 106 | """ |
---|
| 107 | Set the volume of a pixel in (A^3) unit |
---|
| 108 | :Param volume: pixel volume [float] |
---|
| 109 | """ |
---|
[235f514] | 110 | if self.data_vol is None: |
---|
[574adc7] | 111 | raise TypeError("data_vol is missing") |
---|
[51f14603] | 112 | self.data_vol = volume |
---|
[e3f77d8b] | 113 | |
---|
| 114 | def set_is_avg(self, is_avg=False): |
---|
[51f14603] | 115 | """ |
---|
| 116 | Sets is_avg: [bool] |
---|
| 117 | """ |
---|
| 118 | self.is_avg = is_avg |
---|
[e3f77d8b] | 119 | |
---|
[54b0650] | 120 | def _gen(self, qx, qy): |
---|
[51f14603] | 121 | """ |
---|
| 122 | Evaluate the function |
---|
| 123 | :Param x: array of x-values |
---|
| 124 | :Param y: array of y-values |
---|
[e3f77d8b] | 125 | :Param i: array of initial i-value |
---|
[51f14603] | 126 | :return: function value |
---|
| 127 | """ |
---|
| 128 | pos_x = self.data_x |
---|
| 129 | pos_y = self.data_y |
---|
| 130 | pos_z = self.data_z |
---|
[235f514] | 131 | if self.is_avg is None: |
---|
[51f14603] | 132 | pos_x, pos_y, pos_z = transform_center(pos_x, pos_y, pos_z) |
---|
| 133 | sldn = copy.deepcopy(self.data_sldn) |
---|
| 134 | sldn -= self.params['solvent_SLD'] |
---|
[144e032a] | 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) |
---|
[54b0650] | 148 | if len(qy): |
---|
| 149 | qx, qy = _vec(qx), _vec(qy) |
---|
| 150 | I_out = np.empty_like(qx) |
---|
[a1daf86] | 151 | #print("npoints", qx.shape, "npixels", pos_x.shape) |
---|
[54b0650] | 152 | mod.genicomXY(model, qx, qy, I_out) |
---|
[a1daf86] | 153 | #print("I_out after", I_out) |
---|
[51f14603] | 154 | else: |
---|
[54b0650] | 155 | qx = _vec(qx) |
---|
| 156 | I_out = np.empty_like(qx) |
---|
| 157 | mod.genicom(model, qx, I_out) |
---|
[51f14603] | 158 | vol_correction = self.data_total_volume / self.params['total_volume'] |
---|
[54b0650] | 159 | result = (self.params['scale'] * vol_correction * I_out |
---|
| 160 | + self.params['background']) |
---|
| 161 | return result |
---|
[e3f77d8b] | 162 | |
---|
| 163 | def set_sld_data(self, sld_data=None): |
---|
[51f14603] | 164 | """ |
---|
| 165 | Sets sld_data |
---|
| 166 | """ |
---|
| 167 | self.sld_data = sld_data |
---|
| 168 | self.data_pos_unit = sld_data.pos_unit |
---|
[54b0650] | 169 | self.data_x = _vec(sld_data.pos_x) |
---|
| 170 | self.data_y = _vec(sld_data.pos_y) |
---|
| 171 | self.data_z = _vec(sld_data.pos_z) |
---|
| 172 | self.data_sldn = _vec(sld_data.sld_n) |
---|
| 173 | self.data_mx = _vec(sld_data.sld_mx) |
---|
| 174 | self.data_my = _vec(sld_data.sld_my) |
---|
| 175 | self.data_mz = _vec(sld_data.sld_mz) |
---|
| 176 | self.data_vol = _vec(sld_data.vol_pix) |
---|
[51f14603] | 177 | self.data_total_volume = sum(sld_data.vol_pix) |
---|
| 178 | self.params['total_volume'] = sum(sld_data.vol_pix) |
---|
[e3f77d8b] | 179 | |
---|
[51f14603] | 180 | def getProfile(self): |
---|
| 181 | """ |
---|
[e3f77d8b] | 182 | Get SLD profile |
---|
[51f14603] | 183 | : return: sld_data |
---|
| 184 | """ |
---|
| 185 | return self.sld_data |
---|
[e3f77d8b] | 186 | |
---|
| 187 | def run(self, x=0.0): |
---|
| 188 | """ |
---|
[51f14603] | 189 | Evaluate the model |
---|
| 190 | :param x: simple value |
---|
| 191 | :return: (I value) |
---|
| 192 | """ |
---|
[54b0650] | 193 | if isinstance(x, list): |
---|
[51f14603] | 194 | if len(x[1]) > 0: |
---|
| 195 | msg = "Not a 1D." |
---|
[574adc7] | 196 | raise ValueError(msg) |
---|
[51f14603] | 197 | # 1D I is found at y =0 in the 2D pattern |
---|
[54b0650] | 198 | out = self._gen(x[0], []) |
---|
[51f14603] | 199 | return out |
---|
| 200 | else: |
---|
| 201 | msg = "Q must be given as list of qx's and qy's" |
---|
[574adc7] | 202 | raise ValueError(msg) |
---|
[e3f77d8b] | 203 | |
---|
| 204 | def runXY(self, x=0.0): |
---|
| 205 | """ |
---|
[51f14603] | 206 | Evaluate the model |
---|
| 207 | :param x: simple value |
---|
| 208 | :return: I value |
---|
| 209 | :Use this runXY() for the computation |
---|
| 210 | """ |
---|
[54b0650] | 211 | if isinstance(x, list): |
---|
| 212 | return self._gen(x[0], x[1]) |
---|
[51f14603] | 213 | else: |
---|
| 214 | msg = "Q must be given as list of qx's and qy's" |
---|
[574adc7] | 215 | raise ValueError(msg) |
---|
[e3f77d8b] | 216 | |
---|
[51f14603] | 217 | def evalDistribution(self, qdist): |
---|
| 218 | """ |
---|
| 219 | Evaluate a distribution of q-values. |
---|
[ac7be54] | 220 | |
---|
| 221 | :param qdist: ndarray of scalar q-values (for 1D) or list [qx,qy] |
---|
| 222 | where qx,qy are 1D ndarrays (for 2D). |
---|
[51f14603] | 223 | """ |
---|
[54b0650] | 224 | if isinstance(qdist, list): |
---|
| 225 | return self.run(qdist) if len(qdist[1]) < 1 else self.runXY(qdist) |
---|
[51f14603] | 226 | else: |
---|
| 227 | mesg = "evalDistribution is expecting an ndarray of " |
---|
| 228 | mesg += "a list [qx,qy] where qx,qy are arrays." |
---|
[574adc7] | 229 | raise RuntimeError(mesg) |
---|
[e3f77d8b] | 230 | |
---|
[54b0650] | 231 | def _vec(v): |
---|
| 232 | return np.ascontiguousarray(v, 'd') |
---|
| 233 | |
---|
[e3f77d8b] | 234 | class OMF2SLD(object): |
---|
[51f14603] | 235 | """ |
---|
| 236 | Convert OMFData to MAgData |
---|
| 237 | """ |
---|
[e3f77d8b] | 238 | def __init__(self): |
---|
[51f14603] | 239 | """ |
---|
| 240 | Init |
---|
| 241 | """ |
---|
| 242 | self.pos_x = None |
---|
| 243 | self.pos_y = None |
---|
| 244 | self.pos_z = None |
---|
| 245 | self.mx = None |
---|
| 246 | self.my = None |
---|
| 247 | self.mz = None |
---|
| 248 | self.sld_n = None |
---|
| 249 | self.vol_pix = None |
---|
| 250 | self.output = None |
---|
| 251 | self.omfdata = None |
---|
[e3f77d8b] | 252 | |
---|
[51f14603] | 253 | def set_data(self, omfdata, shape='rectangular'): |
---|
| 254 | """ |
---|
| 255 | Set all data |
---|
| 256 | """ |
---|
| 257 | self.omfdata = omfdata |
---|
| 258 | length = int(omfdata.xnodes * omfdata.ynodes * omfdata.znodes) |
---|
[9a5097c] | 259 | pos_x = np.arange(omfdata.xmin, |
---|
[e3f77d8b] | 260 | omfdata.xnodes*omfdata.xstepsize + omfdata.xmin, |
---|
[51f14603] | 261 | omfdata.xstepsize) |
---|
[9a5097c] | 262 | pos_y = np.arange(omfdata.ymin, |
---|
[e3f77d8b] | 263 | omfdata.ynodes*omfdata.ystepsize + omfdata.ymin, |
---|
[51f14603] | 264 | omfdata.ystepsize) |
---|
[9a5097c] | 265 | pos_z = np.arange(omfdata.zmin, |
---|
[e3f77d8b] | 266 | omfdata.znodes*omfdata.zstepsize + omfdata.zmin, |
---|
[51f14603] | 267 | omfdata.zstepsize) |
---|
[9a5097c] | 268 | self.pos_x = np.tile(pos_x, int(omfdata.ynodes * omfdata.znodes)) |
---|
[51f14603] | 269 | self.pos_y = pos_y.repeat(int(omfdata.xnodes)) |
---|
[9a5097c] | 270 | self.pos_y = np.tile(self.pos_y, int(omfdata.znodes)) |
---|
[51f14603] | 271 | self.pos_z = pos_z.repeat(int(omfdata.xnodes * omfdata.ynodes)) |
---|
| 272 | self.mx = omfdata.mx |
---|
| 273 | self.my = omfdata.my |
---|
| 274 | self.mz = omfdata.mz |
---|
[9a5097c] | 275 | self.sld_n = np.zeros(length) |
---|
[e3f77d8b] | 276 | |
---|
[235f514] | 277 | if omfdata.mx is None: |
---|
[9a5097c] | 278 | self.mx = np.zeros(length) |
---|
[235f514] | 279 | if omfdata.my is None: |
---|
[9a5097c] | 280 | self.my = np.zeros(length) |
---|
[235f514] | 281 | if omfdata.mz is None: |
---|
[9a5097c] | 282 | self.mz = np.zeros(length) |
---|
[e3f77d8b] | 283 | |
---|
[51f14603] | 284 | self._check_data_length(length) |
---|
| 285 | self.remove_null_points(False, False) |
---|
[9a5097c] | 286 | mask = np.ones(len(self.sld_n), dtype=bool) |
---|
[51f14603] | 287 | if shape.lower() == 'ellipsoid': |
---|
| 288 | try: |
---|
| 289 | # Pixel (step) size included |
---|
| 290 | x_c = max(self.pos_x) + min(self.pos_x) |
---|
| 291 | y_c = max(self.pos_y) + min(self.pos_y) |
---|
| 292 | z_c = max(self.pos_z) + min(self.pos_z) |
---|
| 293 | x_d = max(self.pos_x) - min(self.pos_x) |
---|
| 294 | y_d = max(self.pos_y) - min(self.pos_y) |
---|
| 295 | z_d = max(self.pos_z) - min(self.pos_z) |
---|
| 296 | x_r = (x_d + omfdata.xstepsize) / 2.0 |
---|
| 297 | y_r = (y_d + omfdata.ystepsize) / 2.0 |
---|
| 298 | z_r = (z_d + omfdata.zstepsize) / 2.0 |
---|
| 299 | x_dir2 = ((self.pos_x - x_c / 2.0) / x_r) |
---|
| 300 | x_dir2 *= x_dir2 |
---|
| 301 | y_dir2 = ((self.pos_y - y_c / 2.0) / y_r) |
---|
| 302 | y_dir2 *= y_dir2 |
---|
| 303 | z_dir2 = ((self.pos_z - z_c / 2.0) / z_r) |
---|
| 304 | z_dir2 *= z_dir2 |
---|
| 305 | mask = (x_dir2 + y_dir2 + z_dir2) <= 1.0 |
---|
[1d014cb] | 306 | except Exception: |
---|
[c155a16] | 307 | logger.error(sys.exc_value) |
---|
[e3f77d8b] | 308 | self.output = MagSLD(self.pos_x[mask], self.pos_y[mask], |
---|
| 309 | self.pos_z[mask], self.sld_n[mask], |
---|
[51f14603] | 310 | self.mx[mask], self.my[mask], self.mz[mask]) |
---|
| 311 | self.output.set_pix_type('pixel') |
---|
| 312 | self.output.set_pixel_symbols('pixel') |
---|
[e3f77d8b] | 313 | |
---|
[51f14603] | 314 | def get_omfdata(self): |
---|
| 315 | """ |
---|
| 316 | Return all data |
---|
| 317 | """ |
---|
| 318 | return self.omfdata |
---|
[e3f77d8b] | 319 | |
---|
[51f14603] | 320 | def get_output(self): |
---|
| 321 | """ |
---|
| 322 | Return output |
---|
| 323 | """ |
---|
| 324 | return self.output |
---|
[e3f77d8b] | 325 | |
---|
[51f14603] | 326 | def _check_data_length(self, length): |
---|
| 327 | """ |
---|
| 328 | Check if the data lengths are consistent |
---|
| 329 | :Params length: data length |
---|
| 330 | """ |
---|
[574adc7] | 331 | parts = (self.pos_x, self.pos_y, self.pos_z, self.mx, self.my, self.mz) |
---|
| 332 | if any(len(v) != length for v in parts): |
---|
| 333 | raise ValueError("Error: Inconsistent data length.") |
---|
[e3f77d8b] | 334 | |
---|
[51f14603] | 335 | def remove_null_points(self, remove=False, recenter=False): |
---|
| 336 | """ |
---|
| 337 | Removes any mx, my, and mz = 0 points |
---|
| 338 | """ |
---|
| 339 | if remove: |
---|
[9a5097c] | 340 | is_nonzero = (np.fabs(self.mx) + np.fabs(self.my) + |
---|
| 341 | np.fabs(self.mz)).nonzero() |
---|
[e3f77d8b] | 342 | if len(is_nonzero[0]) > 0: |
---|
[51f14603] | 343 | self.pos_x = self.pos_x[is_nonzero] |
---|
| 344 | self.pos_y = self.pos_y[is_nonzero] |
---|
| 345 | self.pos_z = self.pos_z[is_nonzero] |
---|
| 346 | self.sld_n = self.sld_n[is_nonzero] |
---|
| 347 | self.mx = self.mx[is_nonzero] |
---|
| 348 | self.my = self.my[is_nonzero] |
---|
| 349 | self.mz = self.mz[is_nonzero] |
---|
| 350 | if recenter: |
---|
| 351 | self.pos_x -= (min(self.pos_x) + max(self.pos_x)) / 2.0 |
---|
| 352 | self.pos_y -= (min(self.pos_y) + max(self.pos_y)) / 2.0 |
---|
| 353 | self.pos_z -= (min(self.pos_z) + max(self.pos_z)) / 2.0 |
---|
[e3f77d8b] | 354 | |
---|
[51f14603] | 355 | def get_magsld(self): |
---|
| 356 | """ |
---|
| 357 | return MagSLD |
---|
| 358 | """ |
---|
| 359 | return self.output |
---|
[e3f77d8b] | 360 | |
---|
| 361 | |
---|
| 362 | class OMFReader(object): |
---|
[51f14603] | 363 | """ |
---|
| 364 | Class to load omf/ascii files (3 columns w/header). |
---|
| 365 | """ |
---|
| 366 | ## File type |
---|
| 367 | type_name = "OMF ASCII" |
---|
[e3f77d8b] | 368 | |
---|
[51f14603] | 369 | ## Wildcards |
---|
| 370 | type = ["OMF files (*.OMF, *.omf)|*.omf"] |
---|
| 371 | ## List of allowed extensions |
---|
| 372 | ext = ['.omf', '.OMF'] |
---|
[e3f77d8b] | 373 | |
---|
[51f14603] | 374 | def read(self, path): |
---|
| 375 | """ |
---|
| 376 | Load data file |
---|
| 377 | :param path: file path |
---|
| 378 | :return: x, y, z, sld_n, sld_mx, sld_my, sld_mz |
---|
| 379 | """ |
---|
| 380 | desc = "" |
---|
[9a5097c] | 381 | mx = np.zeros(0) |
---|
| 382 | my = np.zeros(0) |
---|
| 383 | mz = np.zeros(0) |
---|
[51f14603] | 384 | try: |
---|
| 385 | input_f = open(path, 'rb') |
---|
[1d014cb] | 386 | buff = decode(input_f.read()) |
---|
[51f14603] | 387 | lines = buff.split('\n') |
---|
| 388 | input_f.close() |
---|
| 389 | output = OMFData() |
---|
| 390 | valueunit = None |
---|
| 391 | for line in lines: |
---|
[1d014cb] | 392 | line = line.strip() |
---|
[51f14603] | 393 | # Read data |
---|
[1d014cb] | 394 | if line and not line.startswith('#'): |
---|
| 395 | try: |
---|
| 396 | toks = line.split() |
---|
| 397 | _mx = float(toks[0]) |
---|
| 398 | _my = float(toks[1]) |
---|
| 399 | _mz = float(toks[2]) |
---|
| 400 | _mx = mag2sld(_mx, valueunit) |
---|
| 401 | _my = mag2sld(_my, valueunit) |
---|
| 402 | _mz = mag2sld(_mz, valueunit) |
---|
| 403 | mx = np.append(mx, _mx) |
---|
| 404 | my = np.append(my, _my) |
---|
| 405 | mz = np.append(mz, _mz) |
---|
| 406 | except Exception as exc: |
---|
| 407 | # Skip non-data lines |
---|
| 408 | logger.error(str(exc)+" when processing %r"%line) |
---|
[51f14603] | 409 | #Reading Header; Segment count ignored |
---|
| 410 | s_line = line.split(":", 1) |
---|
| 411 | if s_line[0].lower().count("oommf") > 0: |
---|
| 412 | oommf = s_line[1].lstrip() |
---|
| 413 | if s_line[0].lower().count("title") > 0: |
---|
| 414 | title = s_line[1].lstrip() |
---|
| 415 | if s_line[0].lower().count("desc") > 0: |
---|
| 416 | desc += s_line[1].lstrip() |
---|
| 417 | desc += '\n' |
---|
| 418 | if s_line[0].lower().count("meshtype") > 0: |
---|
| 419 | meshtype = s_line[1].lstrip() |
---|
| 420 | if s_line[0].lower().count("meshunit") > 0: |
---|
| 421 | meshunit = s_line[1].lstrip() |
---|
| 422 | if meshunit.count("m") < 1: |
---|
| 423 | msg = "Error: \n" |
---|
| 424 | msg += "We accept only m as meshunit" |
---|
[574adc7] | 425 | raise ValueError(msg) |
---|
[51f14603] | 426 | if s_line[0].lower().count("xbase") > 0: |
---|
| 427 | xbase = s_line[1].lstrip() |
---|
| 428 | if s_line[0].lower().count("ybase") > 0: |
---|
| 429 | ybase = s_line[1].lstrip() |
---|
| 430 | if s_line[0].lower().count("zbase") > 0: |
---|
| 431 | zbase = s_line[1].lstrip() |
---|
| 432 | if s_line[0].lower().count("xstepsize") > 0: |
---|
| 433 | xstepsize = s_line[1].lstrip() |
---|
| 434 | if s_line[0].lower().count("ystepsize") > 0: |
---|
| 435 | ystepsize = s_line[1].lstrip() |
---|
| 436 | if s_line[0].lower().count("zstepsize") > 0: |
---|
| 437 | zstepsize = s_line[1].lstrip() |
---|
| 438 | if s_line[0].lower().count("xnodes") > 0: |
---|
| 439 | xnodes = s_line[1].lstrip() |
---|
| 440 | if s_line[0].lower().count("ynodes") > 0: |
---|
| 441 | ynodes = s_line[1].lstrip() |
---|
| 442 | if s_line[0].lower().count("znodes") > 0: |
---|
| 443 | znodes = s_line[1].lstrip() |
---|
| 444 | if s_line[0].lower().count("xmin") > 0: |
---|
| 445 | xmin = s_line[1].lstrip() |
---|
| 446 | if s_line[0].lower().count("ymin") > 0: |
---|
| 447 | ymin = s_line[1].lstrip() |
---|
| 448 | if s_line[0].lower().count("zmin") > 0: |
---|
| 449 | zmin = s_line[1].lstrip() |
---|
| 450 | if s_line[0].lower().count("xmax") > 0: |
---|
| 451 | xmax = s_line[1].lstrip() |
---|
| 452 | if s_line[0].lower().count("ymax") > 0: |
---|
| 453 | ymax = s_line[1].lstrip() |
---|
| 454 | if s_line[0].lower().count("zmax") > 0: |
---|
| 455 | zmax = s_line[1].lstrip() |
---|
| 456 | if s_line[0].lower().count("valueunit") > 0: |
---|
| 457 | valueunit = s_line[1].lstrip().rstrip() |
---|
| 458 | if s_line[0].lower().count("valuemultiplier") > 0: |
---|
| 459 | valuemultiplier = s_line[1].lstrip() |
---|
| 460 | if s_line[0].lower().count("valuerangeminmag") > 0: |
---|
| 461 | valuerangeminmag = s_line[1].lstrip() |
---|
| 462 | if s_line[0].lower().count("valuerangemaxmag") > 0: |
---|
| 463 | valuerangemaxmag = s_line[1].lstrip() |
---|
| 464 | if s_line[0].lower().count("end") > 0: |
---|
| 465 | output.filename = os.path.basename(path) |
---|
| 466 | output.oommf = oommf |
---|
| 467 | output.title = title |
---|
| 468 | output.desc = desc |
---|
| 469 | output.meshtype = meshtype |
---|
| 470 | output.xbase = float(xbase) * METER2ANG |
---|
| 471 | output.ybase = float(ybase) * METER2ANG |
---|
| 472 | output.zbase = float(zbase) * METER2ANG |
---|
| 473 | output.xstepsize = float(xstepsize) * METER2ANG |
---|
| 474 | output.ystepsize = float(ystepsize) * METER2ANG |
---|
| 475 | output.zstepsize = float(zstepsize) * METER2ANG |
---|
| 476 | output.xnodes = float(xnodes) |
---|
| 477 | output.ynodes = float(ynodes) |
---|
| 478 | output.znodes = float(znodes) |
---|
| 479 | output.xmin = float(xmin) * METER2ANG |
---|
| 480 | output.ymin = float(ymin) * METER2ANG |
---|
| 481 | output.zmin = float(zmin) * METER2ANG |
---|
| 482 | output.xmax = float(xmax) * METER2ANG |
---|
| 483 | output.ymax = float(ymax) * METER2ANG |
---|
| 484 | output.zmax = float(zmax) * METER2ANG |
---|
| 485 | output.valuemultiplier = valuemultiplier |
---|
| 486 | output.valuerangeminmag = mag2sld(float(valuerangeminmag), \ |
---|
| 487 | valueunit) |
---|
| 488 | output.valuerangemaxmag = mag2sld(float(valuerangemaxmag), \ |
---|
| 489 | valueunit) |
---|
| 490 | output.set_m(mx, my, mz) |
---|
| 491 | return output |
---|
[1d014cb] | 492 | except Exception: |
---|
[51f14603] | 493 | msg = "%s is not supported: \n" % path |
---|
| 494 | msg += "We accept only Text format OMF file." |
---|
[574adc7] | 495 | raise RuntimeError(msg) |
---|
[51f14603] | 496 | |
---|
[e3f77d8b] | 497 | class PDBReader(object): |
---|
[51f14603] | 498 | """ |
---|
| 499 | PDB reader class: limited for reading the lines starting with 'ATOM' |
---|
| 500 | """ |
---|
| 501 | type_name = "PDB" |
---|
| 502 | ## Wildcards |
---|
| 503 | type = ["pdb files (*.PDB, *.pdb)|*.pdb"] |
---|
| 504 | ## List of allowed extensions |
---|
[e3f77d8b] | 505 | ext = ['.pdb', '.PDB'] |
---|
| 506 | |
---|
[51f14603] | 507 | def read(self, path): |
---|
| 508 | """ |
---|
| 509 | Load data file |
---|
[e3f77d8b] | 510 | |
---|
[51f14603] | 511 | :param path: file path |
---|
| 512 | :return: MagSLD |
---|
| 513 | :raise RuntimeError: when the file can't be opened |
---|
| 514 | """ |
---|
[9a5097c] | 515 | pos_x = np.zeros(0) |
---|
| 516 | pos_y = np.zeros(0) |
---|
| 517 | pos_z = np.zeros(0) |
---|
| 518 | sld_n = np.zeros(0) |
---|
| 519 | sld_mx = np.zeros(0) |
---|
| 520 | sld_my = np.zeros(0) |
---|
| 521 | sld_mz = np.zeros(0) |
---|
| 522 | vol_pix = np.zeros(0) |
---|
| 523 | pix_symbol = np.zeros(0) |
---|
[51f14603] | 524 | x_line = [] |
---|
| 525 | y_line = [] |
---|
| 526 | z_line = [] |
---|
| 527 | x_lines = [] |
---|
| 528 | y_lines = [] |
---|
[e3f77d8b] | 529 | z_lines = [] |
---|
[51f14603] | 530 | try: |
---|
| 531 | input_f = open(path, 'rb') |
---|
[1d014cb] | 532 | buff = decode(input_f.read()) |
---|
[51f14603] | 533 | lines = buff.split('\n') |
---|
| 534 | input_f.close() |
---|
| 535 | num = 0 |
---|
| 536 | for line in lines: |
---|
| 537 | try: |
---|
| 538 | # check if line starts with "ATOM" |
---|
| 539 | if line[0:6].strip().count('ATM') > 0 or \ |
---|
| 540 | line[0:6].strip() == 'ATOM': |
---|
| 541 | # define fields of interest |
---|
| 542 | atom_name = line[12:16].strip() |
---|
| 543 | try: |
---|
| 544 | float(line[12]) |
---|
| 545 | atom_name = atom_name[1].upper() |
---|
[1d014cb] | 546 | except Exception: |
---|
[51f14603] | 547 | if len(atom_name) == 4: |
---|
| 548 | atom_name = atom_name[0].upper() |
---|
[e3f77d8b] | 549 | elif line[12] != ' ': |
---|
[51f14603] | 550 | atom_name = atom_name[0].upper() + \ |
---|
[e3f77d8b] | 551 | atom_name[1].lower() |
---|
[51f14603] | 552 | else: |
---|
| 553 | atom_name = atom_name[0].upper() |
---|
| 554 | _pos_x = float(line[30:38].strip()) |
---|
| 555 | _pos_y = float(line[38:46].strip()) |
---|
| 556 | _pos_z = float(line[46:54].strip()) |
---|
[9a5097c] | 557 | pos_x = np.append(pos_x, _pos_x) |
---|
| 558 | pos_y = np.append(pos_y, _pos_y) |
---|
| 559 | pos_z = np.append(pos_z, _pos_z) |
---|
[51f14603] | 560 | try: |
---|
| 561 | val = nsf.neutron_sld(atom_name)[0] |
---|
| 562 | # sld in Ang^-2 unit |
---|
| 563 | val *= 1.0e-6 |
---|
[9a5097c] | 564 | sld_n = np.append(sld_n, val) |
---|
[51f14603] | 565 | atom = formula(atom_name) |
---|
| 566 | # cm to A units |
---|
| 567 | vol = 1.0e+24 * atom.mass / atom.density / NA |
---|
[9a5097c] | 568 | vol_pix = np.append(vol_pix, vol) |
---|
[1d014cb] | 569 | except Exception: |
---|
[b58265c3] | 570 | logger.error("Error: set the sld of %s to zero"% atom_name) |
---|
[9a5097c] | 571 | sld_n = np.append(sld_n, 0.0) |
---|
| 572 | sld_mx = np.append(sld_mx, 0) |
---|
| 573 | sld_my = np.append(sld_my, 0) |
---|
| 574 | sld_mz = np.append(sld_mz, 0) |
---|
| 575 | pix_symbol = np.append(pix_symbol, atom_name) |
---|
[e3f77d8b] | 576 | elif line[0:6].strip().count('CONECT') > 0: |
---|
[51f14603] | 577 | toks = line.split() |
---|
[e3f77d8b] | 578 | num = int(toks[1]) - 1 |
---|
| 579 | val_list = [] |
---|
[51f14603] | 580 | for val in toks[2:]: |
---|
| 581 | try: |
---|
| 582 | int_val = int(val) |
---|
[1d014cb] | 583 | except Exception: |
---|
[51f14603] | 584 | break |
---|
| 585 | if int_val == 0: |
---|
| 586 | break |
---|
| 587 | val_list.append(int_val) |
---|
[e3f77d8b] | 588 | #need val_list ordered |
---|
[51f14603] | 589 | for val in val_list: |
---|
| 590 | index = val - 1 |
---|
| 591 | if (pos_x[index], pos_x[num]) in x_line and \ |
---|
| 592 | (pos_y[index], pos_y[num]) in y_line and \ |
---|
| 593 | (pos_z[index], pos_z[num]) in z_line: |
---|
| 594 | continue |
---|
| 595 | x_line.append((pos_x[num], pos_x[index])) |
---|
| 596 | y_line.append((pos_y[num], pos_y[index])) |
---|
| 597 | z_line.append((pos_z[num], pos_z[index])) |
---|
| 598 | if len(x_line) > 0: |
---|
[e3f77d8b] | 599 | x_lines.append(x_line) |
---|
[51f14603] | 600 | y_lines.append(y_line) |
---|
[e3f77d8b] | 601 | z_lines.append(z_line) |
---|
[1d014cb] | 602 | except Exception: |
---|
[c155a16] | 603 | logger.error(sys.exc_value) |
---|
[e3f77d8b] | 604 | |
---|
[51f14603] | 605 | output = MagSLD(pos_x, pos_y, pos_z, sld_n, sld_mx, sld_my, sld_mz) |
---|
| 606 | output.set_conect_lines(x_line, y_line, z_line) |
---|
| 607 | output.filename = os.path.basename(path) |
---|
| 608 | output.set_pix_type('atom') |
---|
| 609 | output.set_pixel_symbols(pix_symbol) |
---|
| 610 | output.set_nodes() |
---|
| 611 | output.set_pixel_volumes(vol_pix) |
---|
| 612 | output.sld_unit = '1/A^(2)' |
---|
| 613 | return output |
---|
[1d014cb] | 614 | except Exception: |
---|
[574adc7] | 615 | raise RuntimeError("%s is not a sld file" % path) |
---|
[51f14603] | 616 | |
---|
| 617 | def write(self, path, data): |
---|
| 618 | """ |
---|
| 619 | Write |
---|
| 620 | """ |
---|
[9c3d784] | 621 | print("Not implemented... ") |
---|
[51f14603] | 622 | |
---|
[e3f77d8b] | 623 | class SLDReader(object): |
---|
[51f14603] | 624 | """ |
---|
| 625 | Class to load ascii files (7 columns). |
---|
| 626 | """ |
---|
| 627 | ## File type |
---|
| 628 | type_name = "SLD ASCII" |
---|
| 629 | ## Wildcards |
---|
| 630 | type = ["sld files (*.SLD, *.sld)|*.sld", |
---|
| 631 | "txt files (*.TXT, *.txt)|*.txt", |
---|
| 632 | "all files (*.*)|*.*"] |
---|
| 633 | ## List of allowed extensions |
---|
| 634 | ext = ['.sld', '.SLD', '.txt', '.TXT', '.*'] |
---|
| 635 | def read(self, path): |
---|
| 636 | """ |
---|
| 637 | Load data file |
---|
| 638 | :param path: file path |
---|
| 639 | :return MagSLD: x, y, z, sld_n, sld_mx, sld_my, sld_mz |
---|
| 640 | :raise RuntimeError: when the file can't be opened |
---|
| 641 | :raise ValueError: when the length of the data vectors are inconsistent |
---|
| 642 | """ |
---|
| 643 | try: |
---|
[9a5097c] | 644 | pos_x = np.zeros(0) |
---|
| 645 | pos_y = np.zeros(0) |
---|
| 646 | pos_z = np.zeros(0) |
---|
| 647 | sld_n = np.zeros(0) |
---|
| 648 | sld_mx = np.zeros(0) |
---|
| 649 | sld_my = np.zeros(0) |
---|
| 650 | sld_mz = np.zeros(0) |
---|
[51f14603] | 651 | try: |
---|
| 652 | # Use numpy to speed up loading |
---|
[9a5097c] | 653 | input_f = np.loadtxt(path, dtype='float', skiprows=1, |
---|
[51f14603] | 654 | ndmin=1, unpack=True) |
---|
[9a5097c] | 655 | pos_x = np.array(input_f[0]) |
---|
| 656 | pos_y = np.array(input_f[1]) |
---|
| 657 | pos_z = np.array(input_f[2]) |
---|
| 658 | sld_n = np.array(input_f[3]) |
---|
| 659 | sld_mx = np.array(input_f[4]) |
---|
| 660 | sld_my = np.array(input_f[5]) |
---|
| 661 | sld_mz = np.array(input_f[6]) |
---|
[51f14603] | 662 | ncols = len(input_f) |
---|
| 663 | if ncols == 8: |
---|
[9a5097c] | 664 | vol_pix = np.array(input_f[7]) |
---|
[51f14603] | 665 | elif ncols == 7: |
---|
[e3f77d8b] | 666 | vol_pix = None |
---|
[1d014cb] | 667 | except Exception: |
---|
[51f14603] | 668 | # For older version of numpy |
---|
| 669 | input_f = open(path, 'rb') |
---|
[1d014cb] | 670 | buff = decode(input_f.read()) |
---|
[51f14603] | 671 | lines = buff.split('\n') |
---|
| 672 | input_f.close() |
---|
| 673 | for line in lines: |
---|
| 674 | toks = line.split() |
---|
| 675 | try: |
---|
| 676 | _pos_x = float(toks[0]) |
---|
| 677 | _pos_y = float(toks[1]) |
---|
| 678 | _pos_z = float(toks[2]) |
---|
| 679 | _sld_n = float(toks[3]) |
---|
| 680 | _sld_mx = float(toks[4]) |
---|
| 681 | _sld_my = float(toks[5]) |
---|
| 682 | _sld_mz = float(toks[6]) |
---|
[9a5097c] | 683 | pos_x = np.append(pos_x, _pos_x) |
---|
| 684 | pos_y = np.append(pos_y, _pos_y) |
---|
| 685 | pos_z = np.append(pos_z, _pos_z) |
---|
| 686 | sld_n = np.append(sld_n, _sld_n) |
---|
| 687 | sld_mx = np.append(sld_mx, _sld_mx) |
---|
| 688 | sld_my = np.append(sld_my, _sld_my) |
---|
| 689 | sld_mz = np.append(sld_mz, _sld_mz) |
---|
[51f14603] | 690 | try: |
---|
[e3f77d8b] | 691 | _vol_pix = float(toks[7]) |
---|
[9a5097c] | 692 | vol_pix = np.append(vol_pix, _vol_pix) |
---|
[1d014cb] | 693 | except Exception: |
---|
[51f14603] | 694 | vol_pix = None |
---|
[1d014cb] | 695 | except Exception: |
---|
[51f14603] | 696 | # Skip non-data lines |
---|
[c155a16] | 697 | logger.error(sys.exc_value) |
---|
[e3f77d8b] | 698 | output = MagSLD(pos_x, pos_y, pos_z, sld_n, |
---|
[51f14603] | 699 | sld_mx, sld_my, sld_mz) |
---|
| 700 | output.filename = os.path.basename(path) |
---|
| 701 | output.set_pix_type('pixel') |
---|
| 702 | output.set_pixel_symbols('pixel') |
---|
[7432acb] | 703 | if vol_pix is not None: |
---|
[51f14603] | 704 | output.set_pixel_volumes(vol_pix) |
---|
| 705 | return output |
---|
[1d014cb] | 706 | except Exception: |
---|
[574adc7] | 707 | raise RuntimeError("%s is not a sld file" % path) |
---|
[e3f77d8b] | 708 | |
---|
[51f14603] | 709 | def write(self, path, data): |
---|
| 710 | """ |
---|
| 711 | Write sld file |
---|
[e3f77d8b] | 712 | :Param path: file path |
---|
[51f14603] | 713 | :Param data: MagSLD data object |
---|
| 714 | """ |
---|
[235f514] | 715 | if path is None: |
---|
[574adc7] | 716 | raise ValueError("Missing the file path.") |
---|
[235f514] | 717 | if data is None: |
---|
[574adc7] | 718 | raise ValueError("Missing the data to save.") |
---|
[51f14603] | 719 | x_val = data.pos_x |
---|
| 720 | y_val = data.pos_y |
---|
| 721 | z_val = data.pos_z |
---|
| 722 | vol_pix = data.vol_pix |
---|
| 723 | length = len(x_val) |
---|
| 724 | sld_n = data.sld_n |
---|
[235f514] | 725 | if sld_n is None: |
---|
[9a5097c] | 726 | sld_n = np.zeros(length) |
---|
[51f14603] | 727 | sld_mx = data.sld_mx |
---|
[235f514] | 728 | if sld_mx is None: |
---|
[9a5097c] | 729 | sld_mx = np.zeros(length) |
---|
| 730 | sld_my = np.zeros(length) |
---|
| 731 | sld_mz = np.zeros(length) |
---|
[51f14603] | 732 | else: |
---|
| 733 | sld_my = data.sld_my |
---|
| 734 | sld_mz = data.sld_mz |
---|
| 735 | out = open(path, 'w') |
---|
| 736 | # First Line: Column names |
---|
| 737 | out.write("X Y Z SLDN SLDMx SLDMy SLDMz VOLUMEpix") |
---|
| 738 | for ind in range(length): |
---|
[e3f77d8b] | 739 | out.write("\n%g %g %g %g %g %g %g %g" % \ |
---|
| 740 | (x_val[ind], y_val[ind], z_val[ind], sld_n[ind], |
---|
| 741 | sld_mx[ind], sld_my[ind], sld_mz[ind], vol_pix[ind])) |
---|
| 742 | out.close() |
---|
| 743 | |
---|
| 744 | |
---|
| 745 | class OMFData(object): |
---|
[51f14603] | 746 | """ |
---|
| 747 | OMF Data. |
---|
| 748 | """ |
---|
| 749 | _meshunit = "A" |
---|
| 750 | _valueunit = "A^(-2)" |
---|
| 751 | def __init__(self): |
---|
| 752 | """ |
---|
| 753 | Init for mag SLD |
---|
| 754 | """ |
---|
| 755 | self.filename = 'default' |
---|
| 756 | self.oommf = '' |
---|
| 757 | self.title = '' |
---|
| 758 | self.desc = '' |
---|
| 759 | self.meshtype = '' |
---|
| 760 | self.meshunit = self._meshunit |
---|
| 761 | self.valueunit = self._valueunit |
---|
| 762 | self.xbase = 0.0 |
---|
| 763 | self.ybase = 0.0 |
---|
| 764 | self.zbase = 0.0 |
---|
| 765 | self.xstepsize = 6.0 |
---|
| 766 | self.ystepsize = 6.0 |
---|
| 767 | self.zstepsize = 6.0 |
---|
| 768 | self.xnodes = 10.0 |
---|
| 769 | self.ynodes = 10.0 |
---|
| 770 | self.znodes = 10.0 |
---|
| 771 | self.xmin = 0.0 |
---|
| 772 | self.ymin = 0.0 |
---|
| 773 | self.zmin = 0.0 |
---|
| 774 | self.xmax = 60.0 |
---|
| 775 | self.ymax = 60.0 |
---|
| 776 | self.zmax = 60.0 |
---|
| 777 | self.mx = None |
---|
| 778 | self.my = None |
---|
| 779 | self.mz = None |
---|
| 780 | self.valuemultiplier = 1. |
---|
| 781 | self.valuerangeminmag = 0 |
---|
| 782 | self.valuerangemaxmag = 0 |
---|
[e3f77d8b] | 783 | |
---|
[51f14603] | 784 | def __str__(self): |
---|
| 785 | """ |
---|
| 786 | doc strings |
---|
| 787 | """ |
---|
[e3f77d8b] | 788 | _str = "Type: %s\n" % self.__class__.__name__ |
---|
[51f14603] | 789 | _str += "File: %s\n" % self.filename |
---|
| 790 | _str += "OOMMF: %s\n" % self.oommf |
---|
| 791 | _str += "Title: %s\n" % self.title |
---|
| 792 | _str += "Desc: %s\n" % self.desc |
---|
| 793 | _str += "meshtype: %s\n" % self.meshtype |
---|
| 794 | _str += "meshunit: %s\n" % str(self.meshunit) |
---|
| 795 | _str += "xbase: %s [%s]\n" % (str(self.xbase), self.meshunit) |
---|
| 796 | _str += "ybase: %s [%s]\n" % (str(self.ybase), self.meshunit) |
---|
| 797 | _str += "zbase: %s [%s]\n" % (str(self.zbase), self.meshunit) |
---|
[e3f77d8b] | 798 | _str += "xstepsize: %s [%s]\n" % (str(self.xstepsize), |
---|
[51f14603] | 799 | self.meshunit) |
---|
[e3f77d8b] | 800 | _str += "ystepsize: %s [%s]\n" % (str(self.ystepsize), |
---|
[51f14603] | 801 | self.meshunit) |
---|
[e3f77d8b] | 802 | _str += "zstepsize: %s [%s]\n" % (str(self.zstepsize), |
---|
[51f14603] | 803 | self.meshunit) |
---|
| 804 | _str += "xnodes: %s\n" % str(self.xnodes) |
---|
| 805 | _str += "ynodes: %s\n" % str(self.ynodes) |
---|
| 806 | _str += "znodes: %s\n" % str(self.znodes) |
---|
| 807 | _str += "xmin: %s [%s]\n" % (str(self.xmin), self.meshunit) |
---|
| 808 | _str += "ymin: %s [%s]\n" % (str(self.ymin), self.meshunit) |
---|
| 809 | _str += "zmin: %s [%s]\n" % (str(self.zmin), self.meshunit) |
---|
| 810 | _str += "xmax: %s [%s]\n" % (str(self.xmax), self.meshunit) |
---|
| 811 | _str += "ymax: %s [%s]\n" % (str(self.ymax), self.meshunit) |
---|
| 812 | _str += "zmax: %s [%s]\n" % (str(self.zmax), self.meshunit) |
---|
| 813 | _str += "valueunit: %s\n" % self.valueunit |
---|
| 814 | _str += "valuemultiplier: %s\n" % str(self.valuemultiplier) |
---|
[e3f77d8b] | 815 | _str += "ValueRangeMinMag:%s [%s]\n" % (str(self.valuerangeminmag), |
---|
| 816 | self.valueunit) |
---|
| 817 | _str += "ValueRangeMaxMag:%s [%s]\n" % (str(self.valuerangemaxmag), |
---|
| 818 | self.valueunit) |
---|
[51f14603] | 819 | return _str |
---|
[e3f77d8b] | 820 | |
---|
[51f14603] | 821 | def set_m(self, mx, my, mz): |
---|
| 822 | """ |
---|
| 823 | Set the Mx, My, Mz values |
---|
| 824 | """ |
---|
| 825 | self.mx = mx |
---|
| 826 | self.my = my |
---|
| 827 | self.mz = mz |
---|
[e3f77d8b] | 828 | |
---|
| 829 | class MagSLD(object): |
---|
[51f14603] | 830 | """ |
---|
| 831 | Magnetic SLD. |
---|
| 832 | """ |
---|
| 833 | pos_x = None |
---|
| 834 | pos_y = None |
---|
| 835 | pos_z = None |
---|
| 836 | sld_n = None |
---|
| 837 | sld_mx = None |
---|
| 838 | sld_my = None |
---|
| 839 | sld_mz = None |
---|
| 840 | # Units |
---|
| 841 | _pos_unit = 'A' |
---|
| 842 | _sld_unit = '1/A^(2)' |
---|
| 843 | _pix_type = 'pixel' |
---|
[e3f77d8b] | 844 | |
---|
| 845 | def __init__(self, pos_x, pos_y, pos_z, sld_n=None, |
---|
[51f14603] | 846 | sld_mx=None, sld_my=None, sld_mz=None, vol_pix=None): |
---|
| 847 | """ |
---|
| 848 | Init for mag SLD |
---|
| 849 | :params : All should be numpy 1D array |
---|
| 850 | """ |
---|
| 851 | self.is_data = True |
---|
| 852 | self.filename = '' |
---|
| 853 | self.xstepsize = 6.0 |
---|
| 854 | self.ystepsize = 6.0 |
---|
| 855 | self.zstepsize = 6.0 |
---|
| 856 | self.xnodes = 10.0 |
---|
| 857 | self.ynodes = 10.0 |
---|
| 858 | self.znodes = 10.0 |
---|
| 859 | self.has_stepsize = False |
---|
| 860 | self.has_conect = False |
---|
| 861 | self.pos_unit = self._pos_unit |
---|
| 862 | self.sld_unit = self._sld_unit |
---|
| 863 | self.pix_type = 'pixel' |
---|
| 864 | self.pos_x = pos_x |
---|
| 865 | self.pos_y = pos_y |
---|
| 866 | self.pos_z = pos_z |
---|
| 867 | self.sld_n = sld_n |
---|
| 868 | self.line_x = None |
---|
| 869 | self.line_y = None |
---|
| 870 | self.line_z = None |
---|
| 871 | self.sld_mx = sld_mx |
---|
| 872 | self.sld_my = sld_my |
---|
| 873 | self.sld_mz = sld_mz |
---|
| 874 | self.vol_pix = vol_pix |
---|
| 875 | self.sld_m = None |
---|
| 876 | self.sld_phi = None |
---|
| 877 | self.sld_theta = None |
---|
| 878 | self.pix_symbol = None |
---|
[7432acb] | 879 | if sld_mx is not None and sld_my is not None and sld_mz is not None: |
---|
[51f14603] | 880 | self.set_sldms(sld_mx, sld_my, sld_mz) |
---|
| 881 | self.set_nodes() |
---|
[e3f77d8b] | 882 | |
---|
[51f14603] | 883 | def __str__(self): |
---|
| 884 | """ |
---|
| 885 | doc strings |
---|
| 886 | """ |
---|
[e3f77d8b] | 887 | _str = "Type: %s\n" % self.__class__.__name__ |
---|
[51f14603] | 888 | _str += "File: %s\n" % self.filename |
---|
| 889 | _str += "Axis_unit: %s\n" % self.pos_unit |
---|
| 890 | _str += "SLD_unit: %s\n" % self.sld_unit |
---|
| 891 | return _str |
---|
[e3f77d8b] | 892 | |
---|
[51f14603] | 893 | def set_pix_type(self, pix_type): |
---|
| 894 | """ |
---|
| 895 | Set pixel type |
---|
| 896 | :Param pix_type: string, 'pixel' or 'atom' |
---|
| 897 | """ |
---|
| 898 | self.pix_type = pix_type |
---|
[e3f77d8b] | 899 | |
---|
[51f14603] | 900 | def set_sldn(self, sld_n): |
---|
| 901 | """ |
---|
| 902 | Sets neutron SLD |
---|
| 903 | """ |
---|
| 904 | if sld_n.__class__.__name__ == 'float': |
---|
| 905 | if self.is_data: |
---|
[e3f77d8b] | 906 | # For data, put the value to only the pixels w non-zero M |
---|
[9a5097c] | 907 | is_nonzero = (np.fabs(self.sld_mx) + |
---|
| 908 | np.fabs(self.sld_my) + |
---|
| 909 | np.fabs(self.sld_mz)).nonzero() |
---|
| 910 | self.sld_n = np.zeros(len(self.pos_x)) |
---|
[51f14603] | 911 | if len(self.sld_n[is_nonzero]) > 0: |
---|
| 912 | self.sld_n[is_nonzero] = sld_n |
---|
| 913 | else: |
---|
| 914 | self.sld_n.fill(sld_n) |
---|
| 915 | else: |
---|
| 916 | # For non-data, put the value to all the pixels |
---|
[9a5097c] | 917 | self.sld_n = np.ones(len(self.pos_x)) * sld_n |
---|
[51f14603] | 918 | else: |
---|
| 919 | self.sld_n = sld_n |
---|
[e3f77d8b] | 920 | |
---|
[51f14603] | 921 | def set_sldms(self, sld_mx, sld_my, sld_mz): |
---|
[ac7be54] | 922 | r""" |
---|
[5ed76f8] | 923 | Sets mx, my, mz and abs(m). |
---|
| 924 | """ # Note: escaping |
---|
[51f14603] | 925 | if sld_mx.__class__.__name__ == 'float': |
---|
[9a5097c] | 926 | self.sld_mx = np.ones(len(self.pos_x)) * sld_mx |
---|
[51f14603] | 927 | else: |
---|
| 928 | self.sld_mx = sld_mx |
---|
| 929 | if sld_my.__class__.__name__ == 'float': |
---|
[9a5097c] | 930 | self.sld_my = np.ones(len(self.pos_x)) * sld_my |
---|
[51f14603] | 931 | else: |
---|
| 932 | self.sld_my = sld_my |
---|
| 933 | if sld_mz.__class__.__name__ == 'float': |
---|
[9a5097c] | 934 | self.sld_mz = np.ones(len(self.pos_x)) * sld_mz |
---|
[51f14603] | 935 | else: |
---|
| 936 | self.sld_mz = sld_mz |
---|
| 937 | |
---|
[9a5097c] | 938 | sld_m = np.sqrt(sld_mx * sld_mx + sld_my * sld_my + \ |
---|
[e3f77d8b] | 939 | sld_mz * sld_mz) |
---|
[51f14603] | 940 | self.sld_m = sld_m |
---|
[e3f77d8b] | 941 | |
---|
| 942 | def set_pixel_symbols(self, symbol='pixel'): |
---|
[51f14603] | 943 | """ |
---|
| 944 | Set pixel |
---|
| 945 | :Params pixel: str; pixel or atomic symbol, or array of strings |
---|
| 946 | """ |
---|
[235f514] | 947 | if self.sld_n is None: |
---|
[51f14603] | 948 | return |
---|
| 949 | if symbol.__class__.__name__ == 'str': |
---|
[9a5097c] | 950 | self.pix_symbol = np.repeat(symbol, len(self.sld_n)) |
---|
[51f14603] | 951 | else: |
---|
| 952 | self.pix_symbol = symbol |
---|
[e3f77d8b] | 953 | |
---|
| 954 | def set_pixel_volumes(self, vol): |
---|
[51f14603] | 955 | """ |
---|
| 956 | Set pixel volumes |
---|
| 957 | :Params pixel: str; pixel or atomic symbol, or array of strings |
---|
| 958 | """ |
---|
[235f514] | 959 | if self.sld_n is None: |
---|
[51f14603] | 960 | return |
---|
| 961 | if vol.__class__.__name__ == 'ndarray': |
---|
| 962 | self.vol_pix = vol |
---|
| 963 | elif vol.__class__.__name__.count('float') > 0: |
---|
[9a5097c] | 964 | self.vol_pix = np.repeat(vol, len(self.sld_n)) |
---|
[51f14603] | 965 | else: |
---|
| 966 | self.vol_pix = None |
---|
| 967 | |
---|
| 968 | def get_sldn(self): |
---|
| 969 | """ |
---|
| 970 | Returns nuclear sld |
---|
| 971 | """ |
---|
| 972 | return self.sld_n |
---|
[e3f77d8b] | 973 | |
---|
[51f14603] | 974 | def set_nodes(self): |
---|
| 975 | """ |
---|
| 976 | Set xnodes, ynodes, and znodes |
---|
| 977 | """ |
---|
| 978 | self.set_stepsize() |
---|
| 979 | if self.pix_type == 'pixel': |
---|
| 980 | try: |
---|
| 981 | xdist = (max(self.pos_x) - min(self.pos_x)) / self.xstepsize |
---|
| 982 | ydist = (max(self.pos_y) - min(self.pos_y)) / self.ystepsize |
---|
| 983 | zdist = (max(self.pos_z) - min(self.pos_z)) / self.zstepsize |
---|
| 984 | self.xnodes = int(xdist) + 1 |
---|
| 985 | self.ynodes = int(ydist) + 1 |
---|
| 986 | self.znodes = int(zdist) + 1 |
---|
[1d014cb] | 987 | except Exception: |
---|
[51f14603] | 988 | self.xnodes = None |
---|
| 989 | self.ynodes = None |
---|
| 990 | self.znodes = None |
---|
| 991 | else: |
---|
| 992 | self.xnodes = None |
---|
| 993 | self.ynodes = None |
---|
| 994 | self.znodes = None |
---|
[e3f77d8b] | 995 | |
---|
[51f14603] | 996 | def set_stepsize(self): |
---|
| 997 | """ |
---|
| 998 | Set xtepsize, ystepsize, and zstepsize |
---|
| 999 | """ |
---|
| 1000 | if self.pix_type == 'pixel': |
---|
| 1001 | try: |
---|
| 1002 | xpos_pre = self.pos_x[0] |
---|
| 1003 | ypos_pre = self.pos_y[0] |
---|
| 1004 | zpos_pre = self.pos_z[0] |
---|
| 1005 | for x_pos in self.pos_x: |
---|
| 1006 | if xpos_pre != x_pos: |
---|
[9a5097c] | 1007 | self.xstepsize = np.fabs(x_pos - xpos_pre) |
---|
[51f14603] | 1008 | break |
---|
| 1009 | for y_pos in self.pos_y: |
---|
| 1010 | if ypos_pre != y_pos: |
---|
[9a5097c] | 1011 | self.ystepsize = np.fabs(y_pos - ypos_pre) |
---|
[51f14603] | 1012 | break |
---|
| 1013 | for z_pos in self.pos_z: |
---|
| 1014 | if zpos_pre != z_pos: |
---|
[9a5097c] | 1015 | self.zstepsize = np.fabs(z_pos - zpos_pre) |
---|
[51f14603] | 1016 | break |
---|
| 1017 | #default pix volume |
---|
[9a5097c] | 1018 | self.vol_pix = np.ones(len(self.pos_x)) |
---|
[51f14603] | 1019 | vol = self.xstepsize * self.ystepsize * self.zstepsize |
---|
| 1020 | self.set_pixel_volumes(vol) |
---|
| 1021 | self.has_stepsize = True |
---|
[1d014cb] | 1022 | except Exception: |
---|
[51f14603] | 1023 | self.xstepsize = None |
---|
| 1024 | self.ystepsize = None |
---|
| 1025 | self.zstepsize = None |
---|
| 1026 | self.vol_pix = None |
---|
| 1027 | self.has_stepsize = False |
---|
| 1028 | else: |
---|
| 1029 | self.xstepsize = None |
---|
| 1030 | self.ystepsize = None |
---|
| 1031 | self.zstepsize = None |
---|
| 1032 | self.has_stepsize = True |
---|
| 1033 | return self.xstepsize, self.ystepsize, self.zstepsize |
---|
| 1034 | |
---|
| 1035 | def set_conect_lines(self, line_x, line_y, line_z): |
---|
| 1036 | """ |
---|
| 1037 | Set bonding line data if taken from pdb |
---|
| 1038 | """ |
---|
| 1039 | if line_x.__class__.__name__ != 'list' or len(line_x) < 1: |
---|
| 1040 | return |
---|
| 1041 | if line_y.__class__.__name__ != 'list' or len(line_y) < 1: |
---|
| 1042 | return |
---|
| 1043 | if line_z.__class__.__name__ != 'list' or len(line_z) < 1: |
---|
| 1044 | return |
---|
| 1045 | self.has_conect = True |
---|
[e3f77d8b] | 1046 | self.line_x = line_x |
---|
| 1047 | self.line_y = line_y |
---|
[51f14603] | 1048 | self.line_z = line_z |
---|
| 1049 | |
---|
[54b0650] | 1050 | def _get_data_path(*path_parts): |
---|
| 1051 | from os.path import realpath, join as joinpath, dirname, abspath |
---|
| 1052 | # in sas/sascalc/calculator; want sas/sasview/test |
---|
| 1053 | return joinpath(dirname(realpath(__file__)), |
---|
| 1054 | '..', '..', 'sasview', 'test', *path_parts) |
---|
| 1055 | |
---|
[51f14603] | 1056 | def test_load(): |
---|
[b6627d9] | 1057 | """ |
---|
| 1058 | Test code |
---|
| 1059 | """ |
---|
[06aaa75d] | 1060 | from mpl_toolkits.mplot3d import Axes3D |
---|
[54b0650] | 1061 | tfpath = _get_data_path("1d_data", "CoreXY_ShellZ.txt") |
---|
| 1062 | ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf") |
---|
| 1063 | if not os.path.isfile(tfpath) or not os.path.isfile(ofpath): |
---|
| 1064 | raise ValueError("file(s) not found: %r, %r"%(tfpath, ofpath)) |
---|
[51f14603] | 1065 | reader = SLDReader() |
---|
| 1066 | oreader = OMFReader() |
---|
[f926abb] | 1067 | output = reader.read(tfpath) |
---|
| 1068 | ooutput = oreader.read(ofpath) |
---|
[51f14603] | 1069 | foutput = OMF2SLD() |
---|
| 1070 | foutput.set_data(ooutput) |
---|
| 1071 | |
---|
| 1072 | import matplotlib.pyplot as plt |
---|
| 1073 | fig = plt.figure() |
---|
[06aaa75d] | 1074 | ax = Axes3D(fig) |
---|
[e3f77d8b] | 1075 | ax.plot(output.pos_x, output.pos_y, output.pos_z, '.', c="g", |
---|
| 1076 | alpha=0.7, markeredgecolor='gray', rasterized=True) |
---|
[51f14603] | 1077 | gap = 7 |
---|
| 1078 | max_mx = max(output.sld_mx) |
---|
| 1079 | max_my = max(output.sld_my) |
---|
| 1080 | max_mz = max(output.sld_mz) |
---|
[e3f77d8b] | 1081 | max_m = max(max_mx, max_my, max_mz) |
---|
[51f14603] | 1082 | x2 = output.pos_x+output.sld_mx/max_m * gap |
---|
| 1083 | y2 = output.pos_y+output.sld_my/max_m * gap |
---|
| 1084 | z2 = output.pos_z+output.sld_mz/max_m * gap |
---|
[9a5097c] | 1085 | x_arrow = np.column_stack((output.pos_x, x2)) |
---|
| 1086 | y_arrow = np.column_stack((output.pos_y, y2)) |
---|
| 1087 | z_arrow = np.column_stack((output.pos_z, z2)) |
---|
[51f14603] | 1088 | unit_x2 = output.sld_mx / max_m |
---|
| 1089 | unit_y2 = output.sld_my / max_m |
---|
| 1090 | unit_z2 = output.sld_mz / max_m |
---|
[9a5097c] | 1091 | color_x = np.fabs(unit_x2 * 0.8) |
---|
| 1092 | color_y = np.fabs(unit_y2 * 0.8) |
---|
| 1093 | color_z = np.fabs(unit_z2 * 0.8) |
---|
| 1094 | colors = np.column_stack((color_x, color_y, color_z)) |
---|
[51f14603] | 1095 | plt.show() |
---|
[b6627d9] | 1096 | |
---|
[f926abb] | 1097 | def test_save(): |
---|
| 1098 | ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf") |
---|
| 1099 | if not os.path.isfile(ofpath): |
---|
| 1100 | raise ValueError("file(s) not found: %r"%(ofpath,)) |
---|
| 1101 | oreader = OMFReader() |
---|
| 1102 | omfdata = oreader.read(ofpath) |
---|
| 1103 | omf2sld = OMF2SLD() |
---|
| 1104 | omf2sld.set_data(omfdata) |
---|
| 1105 | writer = SLDReader() |
---|
| 1106 | writer.write("out.txt", omf2sld.output) |
---|
| 1107 | |
---|
[51f14603] | 1108 | def test(): |
---|
[b6627d9] | 1109 | """ |
---|
| 1110 | Test code |
---|
| 1111 | """ |
---|
[54b0650] | 1112 | ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf") |
---|
| 1113 | if not os.path.isfile(ofpath): |
---|
| 1114 | raise ValueError("file(s) not found: %r"%(ofpath,)) |
---|
[51f14603] | 1115 | oreader = OMFReader() |
---|
[f926abb] | 1116 | omfdata = oreader.read(ofpath) |
---|
| 1117 | omf2sld = OMF2SLD() |
---|
| 1118 | omf2sld.set_data(omfdata) |
---|
[51f14603] | 1119 | model = GenSAS() |
---|
[f926abb] | 1120 | model.set_sld_data(omf2sld.output) |
---|
| 1121 | x = np.linspace(0, 0.1, 11)[1:] |
---|
[a1daf86] | 1122 | return model.runXY([x, x]) |
---|
[e3f77d8b] | 1123 | |
---|
| 1124 | if __name__ == "__main__": |
---|
[54b0650] | 1125 | #test_load() |
---|
[f926abb] | 1126 | #test_save() |
---|
[a1daf86] | 1127 | #print(test()) |
---|
[51f14603] | 1128 | test() |
---|