diff --git a/python/amunpy/src/amunpy/amun.py b/python/amunpy/src/amunpy/amun.py index 962dcb4..9e612ad 100644 --- a/python/amunpy/src/amunpy/amun.py +++ b/python/amunpy/src/amunpy/amun.py @@ -316,144 +316,139 @@ class Amun: bz = self.__read_binary_data__('magz', chunk_number) vy = self.__read_binary_data__('vely', chunk_number) vz = self.__read_binary_data__('velz', chunk_number) - dset = vy * bz - vz * by + dset = vz * by - vy * bz if self.attributes['resistivity'] > 0: + iy = self.attributes['ndims'] - 1 + tmp = (numpy.roll(bz, -1, axis=iy) - numpy.roll(bz, 1, axis=iy)) if self.attributes['ndims'] == 3: - tmp = (numpy.roll(by, 1, axis=1) - numpy.roll(by, -1, axis=1)) - tmp += (numpy.roll(bz, -1, axis=2) - numpy.roll(bz, 1, axis=2)) - else: - tmp = (numpy.roll(bz, -1, axis=1) - numpy.roll(bz, 1, axis=1)) + iz = self.attributes['ndims'] - 2 + tmp += (numpy.roll(by, 1, axis=iz) - numpy.roll(by, -1, axis=iz)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - dset -= 0.5 * self.attributes['resistivity'] * tmp + dset += (self.attributes['resistivity'] / 2) * tmp elif dataset == 'eley': bx = self.__read_binary_data__('magx', chunk_number) bz = self.__read_binary_data__('magz', chunk_number) vx = self.__read_binary_data__('velx', chunk_number) vz = self.__read_binary_data__('velz', chunk_number) - dset = vz * bx - vx * bz + dset = vx * bz - vz * bx if self.attributes['resistivity'] > 0: + ix = self.attributes['ndims'] + tmp = (numpy.roll(bz, 1, axis=ix) - numpy.roll(bz, -1, axis=ix)) if self.attributes['ndims'] == 3: - tmp = (numpy.roll(bx, -1, axis=1) - numpy.roll(bx, 1, axis=1)) - tmp += (numpy.roll(bz, 1, axis=3) - numpy.roll(bz, -1, axis=3)) - else: - tmp = (numpy.roll(bz, 1, axis=2) - numpy.roll(bz, -1, axis=2)) + iz = self.attributes['ndims'] - 2 + tmp += (numpy.roll(bx, -1, axis=iz) - numpy.roll(bx, 1, axis=iz)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - dset -= 0.5 * self.attributes['resistivity'] * tmp + dset += (self.attributes['resistivity'] / 2) * tmp elif dataset == 'elez': bx = self.__read_binary_data__('magx', chunk_number) by = self.__read_binary_data__('magy', chunk_number) vx = self.__read_binary_data__('velx', chunk_number) vy = self.__read_binary_data__('vely', chunk_number) - dset = vx * by - vy * bx + dset = vy * bx - vx * by if self.attributes['resistivity'] > 0: - if self.attributes['ndims'] == 3: - tmp = (numpy.roll(bx, 1, axis=2) - numpy.roll(bx, -1, axis=2)) - tmp += (numpy.roll(by, -1, axis=3) - numpy.roll(by, 1, axis=3)) - else: - tmp = (numpy.roll(bx, 1, axis=1) - numpy.roll(bx, -1, axis=1)) - tmp += (numpy.roll(by, -1, axis=2) - numpy.roll(by, 1, axis=2)) + ix = self.attributes['ndims'] + iy = self.attributes['ndims'] - 1 + tmp = (numpy.roll(by, -1, axis=ix) - numpy.roll(by, 1, axis=ix)) + tmp += (numpy.roll(bx, 1, axis=iy) - numpy.roll(bx, -1, axis=iy)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - dset -= 0.5 * self.attributes['resistivity'] * tmp + dset += (self.attributes['resistivity'] / 2) * tmp elif dataset == 'elec': b1 = self.__read_binary_data__('magy', chunk_number) b2 = self.__read_binary_data__('magz', chunk_number) v1 = self.__read_binary_data__('vely', chunk_number) v2 = self.__read_binary_data__('velz', chunk_number) - dtmp = v1 * b2 - v2 * b1 + dtmp = v2 * b1 - v1 * b2 if self.attributes['resistivity'] > 0: + iy = self.attributes['ndims'] - 1 + tmp = (numpy.roll(b2, -1, axis=iy) - numpy.roll(b2, 1, axis=iy)) if self.attributes['ndims'] == 3: - tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1)) - tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2)) - else: - tmp = (numpy.roll(b2, -1, axis=1) - numpy.roll(b2, 1, axis=1)) + iz = self.attributes['ndims'] - 2 + tmp += (numpy.roll(b1, 1, axis=iz) - numpy.roll(b1, -1, axis=iz)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - dtmp -= 0.5 * self.attributes['resistivity'] * tmp + dtmp += (self.attributes['resistivity'] / 2) * tmp dset = dtmp**2 b1 = self.__read_binary_data__('magx', chunk_number) v1 = self.__read_binary_data__('velx', chunk_number) - dtmp = v2 * b1 - v1 * b2 + dtmp = v1 * b2 - v2 * b1 if self.attributes['resistivity'] > 0: + ix = self.attributes['ndims'] + tmp = (numpy.roll(b2, 1, axis=ix) - numpy.roll(b2, -1, axis=ix)) if self.attributes['ndims'] == 3: - tmp = (numpy.roll(b1, -1, axis=1) - numpy.roll(b1, 1, axis=1)) - tmp += (numpy.roll(b2, 1, axis=3) - numpy.roll(b2, -1, axis=3)) - else: - tmp = (numpy.roll(b2, 1, axis=2) - numpy.roll(b2, -1, axis=2)) + iz = self.attributes['ndims'] - 2 + tmp += (numpy.roll(b1, -1, axis=iz) - numpy.roll(b1, 1, axis=iz)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - dtmp -= 0.5 * self.attributes['resistivity'] * tmp + dtmp += (self.attributes['resistivity'] / 2) * tmp dset += dtmp**2 b2 = self.__read_binary_data__('magy', chunk_number) v2 = self.__read_binary_data__('vely', chunk_number) - dtmp = v1 * b2 - v2 * b1 + dtmp = v2 * b1 - v1 * b2 if self.attributes['resistivity'] > 0: - if self.attributes['ndims'] == 3: - tmp = (numpy.roll(b1, 1, axis=2) - numpy.roll(b1, -1, axis=2)) - tmp += (numpy.roll(b2, -1, axis=3) - numpy.roll(b2, 1, axis=3)) - else: - tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1)) - tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2)) + ix = self.attributes['ndims'] + iy = self.attributes['ndims'] - 1 + tmp = (numpy.roll(b2, -1, axis=ix) - numpy.roll(b2, 1, axis=ix)) + tmp += (numpy.roll(b1, 1, axis=iy) - numpy.roll(b1, -1, axis=iy)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - dtmp -= 0.5 * self.attributes['resistivity'] * tmp + dtmp += (self.attributes['resistivity'] / 2) * tmp + dset += dtmp**2 - dset = numpy.sqrt(dset) + dset = numpy.sqrt(dset) elif dataset == 'evec': b1 = self.__read_binary_data__('magy', chunk_number) b2 = self.__read_binary_data__('magz', chunk_number) v1 = self.__read_binary_data__('vely', chunk_number) v2 = self.__read_binary_data__('velz', chunk_number) - wx = v1 * b2 - v2 * b1 + wx = v2 * b1 - v1 * b2 if self.attributes['resistivity'] > 0: + iy = self.attributes['ndims'] - 1 + tmp = (numpy.roll(b2, -1, axis=iy) - numpy.roll(b2, 1, axis=iy)) if self.attributes['ndims'] == 3: - tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1)) - tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2)) - else: - tmp = (numpy.roll(b2, -1, axis=1) - numpy.roll(b2, 1, axis=1)) + iz = self.attributes['ndims'] - 2 + tmp += (numpy.roll(b1, 1, axis=iz) - numpy.roll(b1, -1, axis=iz)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - wx -= 0.5 * self.attributes['resistivity'] * tmp + wx += (self.attributes['resistivity'] / 2) * tmp b1 = self.__read_binary_data__('magx', chunk_number) v1 = self.__read_binary_data__('velx', chunk_number) - wy = v2 * b1 - v1 * b2 + wy = v1 * b2 - v2 * b1 if self.attributes['resistivity'] > 0: + ix = self.attributes['ndims'] + tmp = (numpy.roll(b2, 1, axis=ix) - numpy.roll(b2, -1, axis=ix)) if self.attributes['ndims'] == 3: - tmp = (numpy.roll(b1, -1, axis=1) - numpy.roll(b1, 1, axis=1)) - tmp += (numpy.roll(b2, 1, axis=3) - numpy.roll(b2, -1, axis=3)) - else: - tmp = (numpy.roll(b2, 1, axis=2) - numpy.roll(b2, -1, axis=2)) + iz = self.attributes['ndims'] - 2 + tmp += (numpy.roll(b1, -1, axis=iz) - numpy.roll(b1, 1, axis=iz)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - wy -= 0.5 * self.attributes['resistivity'] * tmp + wy += (self.attributes['resistivity'] / 2) * tmp b2 = self.__read_binary_data__('magy', chunk_number) v2 = self.__read_binary_data__('vely', chunk_number) wz = v1 * b2 - v2 * b1 if self.attributes['resistivity'] > 0: - if self.attributes['ndims'] == 3: - tmp = (numpy.roll(b1, 1, axis=2) - numpy.roll(b1, -1, axis=2)) - tmp += (numpy.roll(b2, -1, axis=3) - numpy.roll(b2, 1, axis=3)) - else: - tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1)) - tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2)) + ix = self.attributes['ndims'] + iy = self.attributes['ndims'] - 1 + tmp = (numpy.roll(b2, -1, axis=ix) - numpy.roll(b2, 1, axis=ix)) + tmp += (numpy.roll(b1, 1, axis=iy) - numpy.roll(b1, -1, axis=iy)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - wz -= 0.5 * self.attributes['resistivity'] * tmp + wz += (self.attributes['resistivity'] / 2) * tmp dset = [wx, wy, wz] elif dataset == 'eint':