gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

PYTHON: Correct calculation of electric field.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2022-10-23 15:55:08 -03:00
parent a70a74d231
commit 41bfe45764
1 changed files with 55 additions and 60 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

115
python/amunpy/src/amunpy/amun.py
View File

@ -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':