In [1]:
# !/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on 20191125
@author: zhangji
test the linear relationship
U_t =?= U_sh + U_wm
U_t is the total velocity
U_sh is the velocity induced by shear flow
U_wm is the active velocity.
"""
# %pylab inline
# pylab.rcParams['figure.figsize'] = (25, 11)
# fontsize = 40
# import numpy as np
# import scipy as sp
# from scipy.optimize import leastsq, curve_fit
# from scipy import interpolate
# from scipy.interpolate import interp1d
# from scipy.io import loadmat, savemat
# # import scipy.misc
# import matplotlib
# from matplotlib import pyplot as plt
# from matplotlib import animation, rc
# import matplotlib.ticker as mtick
# from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes
# from mpl_toolkits.mplot3d import Axes3D, axes3d
# from sympy import symbols, simplify, series, exp
# from sympy.matrices import Matrix
# from sympy.solvers import solve
# from IPython.display import display, HTML
# from tqdm import tqdm_notebook as tqdm
# import pandas as pd
# import re
# from scanf import scanf
# import os
# import glob
# from codeStore import support_fun as spf
# from src.support_class import *
# from src import stokes_flow as sf
# rc('animation', html='html5')
# PWD = os.getcwd()
# font = {'size': 20}
# matplotlib.rc('font', **font)
# np.set_printoptions(linewidth=90, precision=5)
%load_ext autoreload
%autoreload 2
import os
import glob
import re
import pandas as pd
from scanf import scanf
import natsort
import numpy as np
import scipy as sp
from scipy.optimize import leastsq, curve_fit
from scipy import interpolate, spatial, sparse
# from scipy.interpolate import interp1d
from scipy.io import loadmat, savemat
# import scipy.misc
import importlib
from IPython.display import display, HTML
import pickle
import matplotlib
from matplotlib import pyplot as plt
from matplotlib import colors as mcolors
from matplotlib import animation, rc
import matplotlib.ticker as mtick
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes
from mpl_toolkits.mplot3d import Axes3D, axes3d
from matplotlib import ticker, cm
from mpl_toolkits.axes_grid1 import make_axes_locatable
from time import time
from src import support_class as spc
from src import jeffery_model as jm
from codeStore import support_fun as spf
from codeStore import support_fun_table as spf_tb
# %matplotlib notebook
%matplotlib inline
rc('animation', html='html5')
rc('text', usetex=True)
params = {'text.latex.preamble': [r'\usepackage{bm}', r'\usepackage{amsmath}']}
plt.rcParams.update(params)
fontsize = 40
figsize = (30, 16)
PWD = os.getcwd()
In [2]:
fig = plt.figure(figsize=(2, 2))
fig.patch.set_facecolor('white')
ax0 = fig.add_subplot(1, 1, 1)
In [3]:
fft_pickle = 'ecoC01B05_tao1_wm0'
In [4]:
# do 2D IFFT myself and compare with numpy version
use_idx_max = 6
# ((psi, ((theta, phi, ui))))
with open('%s.pickle' % fft_pickle, 'rb') as handle:
pickle_data = pickle.load(handle)
ttheta, tphi = pickle_data[0][1][0][:2]
tpsi = np.array([ti[0] for ti in pickle_data])
tw = pickle_data[0][1][3][2].values
tw_fft = np.fft.fft2(tw)
th_freq = np.fft.fftfreq(ttheta.size, 1 / ttheta.size)
ph_freq = np.fft.fftshift(np.fft.fftfreq(tphi.size, 1 / tphi.size))
idx_max = np.dstack(np.unravel_index(np.argsort(np.abs(tw_fft).ravel()), tw_fft.shape))[0][::-1]
# numpy version IFFT
idx = np.zeros_like(tw_fft, dtype=bool)
for ti in idx_max[:use_idx_max]:
idx[ti[0], ti[1]] = 1
tw_fft2 = tw_fft * idx
tw2 = np.fft.ifft2(tw_fft2)
# do IFFT myself
tw3 = np.zeros_like(tw2)
tM, tN = tw_fft.shape
tm, tn = np.meshgrid(np.arange(tM), np.arange(tN), indexing='ij')
for tk, tl in idx_max[:use_idx_max]:
tAkl = tw_fft[tk, tl]
tw3 = tw3 + tAkl * np.exp(2 * np.pi * 1j * (tm * tk / tM + tn * tl / tN)) / (tM * tN)
print(tk, tl, tAkl)
print((np.abs(tw3.imag) / np.abs(tw3)).max())
print((tw3 - tw).max(), (tw3 - tw2).max())
In [5]:
# print values of primary frequence
print_idx_max = 20
# ((psi, ((theta, phi, ui))))
with open('%s.pickle' % fft_pickle, 'rb') as handle:
pickle_data = pickle.load(handle)
ttheta, tphi = pickle_data[0][1][0][:2]
tpsi = np.array([ti[0] for ti in pickle_data])
tw = pickle_data[0][1][3][2].values
tw_fft = np.fft.fft2(tw)
th_freq = np.fft.fftfreq(ttheta.size, 1 / ttheta.size)
ph_freq = np.fft.fftshift(np.fft.fftfreq(tphi.size, 1 / tphi.size))
idx_max = np.dstack(np.unravel_index(np.argsort(np.abs(tw_fft).ravel()), tw_fft.shape))[0][::-1]
for i0, (tk, tl) in enumerate(idx_max[:print_idx_max]):
if not np.mod(i0, 2):
print('*', tk, tl, tw_fft[tk, tl])
else:
print(' ', tk, tl, tw_fft[tk, tl])
In [6]:
# 2D fft of Wx, Wy, Wz
importlib.reload(spf_tb)
use_uidx = 3 # 3: wx, 4: wy, 5: wz
tktl_list = [[0, 62], [1, 62], [31, 62], [1, 0]]
# use_idx_list = [[0], [1, 2], [3, 4], [5, 6]]
# ((psi, ((theta, phi, ui))))
with open('%s.pickle' % fft_pickle, 'rb') as handle:
pickle_data = pickle.load(handle)
ttheta, tphi = pickle_data[0][1][0][:2]
tpsi = np.array([ti[0] for ti in pickle_data])
tw = pickle_data[0][1][use_uidx][2].values
for tktl in tktl_list:
if np.sum(tktl) < 1:
spf_tb.show_fft_major(tw, (tktl, ), ttheta, tphi)
else:
spf_tb.show_fft_fit(tw, tktl, ttheta, tphi)
spf_tb.show_fft_major(tw, tktl_list, ttheta, tphi)
Out[6]:
In [7]:
# 2D fft of Wx, Wy, Wz, investigate Amp and \alpha as function of \psi
importlib.reload(spf_tb)
use_uidx = 3 # 3: wx, 4: wy, 5: wz
tktl = [0, 62]
# ((psi, ((theta, phi, ui))))
with open('%s.pickle' % fft_pickle, 'rb') as handle:
pickle_data = pickle.load(handle)
ttheta, tphi = pickle_data[0][1][0][:2]
tpsi = np.array([ti[0] for ti in pickle_data])
for tpsi, table_psi_data in pickle_data[:]:
tw = table_psi_data[use_uidx][2].values
spf_tb.show_fft_fit(tw, tktl, ttheta, tphi)
In [ ]:
In [384]:
fft_pickle = 'ecoC01B05_tau1_wm0_b'
In [385]:
# do 3D IFFT myself and compare with numpy version
importlib.reload(spf_tb)
use_idx_max = 6
use_uidx = 5 # 3: wx, 4: wy, 5: wz
# ((psi, ((theta, phi, ui), )))
with open('%s.pickle' % fft_pickle, 'rb') as handle:
pickle_data = pickle.load(handle)
ttheta, tphi = pickle_data[0][1][0][:2]
tpsi = np.array([ti[0] for ti in pickle_data])
U_all = [[] for i in range(6)]
for tpsi, table_psi_data in pickle_data:
for (ty, tx, tU), Ui in zip(table_psi_data, U_all):
Ui.append((tpsi, ty, tx, tU))
tw = np.dstack([tU.values for _, _, _, tU in U_all[use_uidx]])
tw_fft = np.fft.fftn(tw)
idx_max = np.dstack(np.unravel_index(np.argsort(np.abs(tw_fft).ravel()), tw_fft.shape))[0][::-1]
# numpy version IFFT
idx = np.zeros_like(tw_fft, dtype=bool)
for ti in idx_max[:use_idx_max]:
idx[ti[0], ti[1], ti[2]] = 1
tw_fft2 = tw_fft * idx
tw2 = np.fft.ifftn(tw_fft2)
print()
# do IFFT myself
tw3 = np.zeros_like(tw2)
tM, tN, tO = tw_fft.shape
tm, tn, to = np.meshgrid(np.arange(tM), np.arange(tN), np.arange(tO), indexing='ij')
for tk, tl, tj in idx_max[:use_idx_max]:
tAkl = tw_fft[tk, tl, tj]
tw3 = tw3 + tAkl * np.exp(2 * np.pi * 1j * tm * tk / tM) / tM * \
np.exp(2 * np.pi * 1j * tn * tl / tN) / tN * \
np.exp(2 * np.pi * 1j * to * tj / tO) / tO
print(tk, tl, tj, tAkl)
print((np.abs(tw3.imag) / np.abs(tw3)).max())
print(np.abs(tw3 - tw2).max())
In [386]:
# print values of primary frequence
importlib.reload(spf_tb)
print_idx_max = 100
use_uidx = 5 # 3: wx, 4: wy, 5: wz
# ((psi, ((theta, phi, ui), )))
with open('%s.pickle' % fft_pickle, 'rb') as handle:
pickle_data = pickle.load(handle)
ttheta, tphi = pickle_data[0][1][0][:2]
tpsi = np.array([ti[0] for ti in pickle_data])
U_all = [[] for i in range(6)]
for tpsi, table_psi_data in pickle_data:
for (ty, tx, tU), Ui in zip(table_psi_data, U_all):
Ui.append((tpsi, ty, tx, tU))
tw = np.dstack([tU.values for _, _, _, tU in U_all[use_uidx]])
tw_fft = np.fft.fftn(tw)
idx_max = np.dstack(np.unravel_index(np.argsort(np.abs(tw_fft).ravel()), tw_fft.shape))[0][::-1]
t1 = []
for i0, (tk, tl, tj) in enumerate(idx_max[:print_idx_max]):
t1.append(np.abs(tw_fft[tk, tl, tj]))
if not np.mod(i0, 2):
print('%5d * [%d, %d, %d] ' % (i0, tk, tl, tj), tw_fft[tk, tl, tj], np.abs(tw_fft[tk, tl, tj]))
else:
print(' [%d, %d, %d] ' % (tk, tl, tj), tw_fft[tk, tl, tj], np.abs(tw_fft[tk, tl, tj]))
In [388]:
tktltj_list3 = ([0, 0, 0], [1, 1, 0], [31, 1, 0], [1, 63, 2], [31, 1, 2], [31, 63, 2], [1, 1, 2], [0, 63, 14],
[0, 63, 2], [0, 63, 0], [30, 63, 0], [30, 1, 0], [30, 63, 14])
t1 = [np.abs(tw_fft[tk, tl, tj]) for tk, tl, tj in idx_max]
plt.plot(t1, '.')
plt.xlim(-1, 50)
# plt.ylim(10, 100)
Out[388]:
In [390]:
importlib.reload(spf_tb)
use_uidx = 4 # 3: wx, 4: wy, 5: wz
tktltj_list1 = ([0, 0, 0], [0, 62, 0], [1, 62, 0], [31, 62, 0], [0, 62, 14], [0, 2, 14], [1, 0, 0], [31, 0, 2],
[31, 0, 14], [31, 62, 2], [31, 2, 14], [1, 2, 2], [31, 2, 2])
tktltj_list2 = ([0, 0, 0], [31, 0, 0], [0, 62, 0], [0, 0, 2], [1, 62, 0], [31, 62, 0], [0, 62, 14], [1, 0, 2],
[31, 0, 2], [0, 2, 14], [31, 62, 2], [31, 2, 14], [1, 2, 2])
tktltj_list3 = ([0, 0, 0], [1, 1, 0], [31, 1, 0], [1, 63, 2], [31, 1, 2], [31, 63, 2], [1, 1, 2], [0, 63, 14],
[0, 63, 2], [0, 63, 0], [30, 63, 0], [30, 1, 0], [30, 63, 14])
tktltj_list = [tktltj_list1, tktltj_list2, tktltj_list3][use_uidx - 3]
# tktltj_list = ([0, 62, 0], [1, 62, 0] )
# ((psi, ((theta, phi, ui), )))
with open('%s.pickle' % fft_pickle, 'rb') as handle:
pickle_data = pickle.load(handle)
ttheta_all, tphi_all = pickle_data[0][1][0][:2]
tpsi_all = np.array([ti[0] for ti in pickle_data])
U_all = [[] for i in range(6)]
for _, table_psi_data in pickle_data:
for (_, _, tU), Ui in zip(table_psi_data, U_all):
Ui.append(tU)
tw = np.dstack(U_all[use_uidx])
# tw_fft, tw2, tw_fft2 = spf_tb.do_3dfft_major_conj(tw, tktltj_list)
tw_fft, tw2, tw_fft2 = spf_tb.do_3dfft_major(tw, tktltj_list)
# spf_tb.show_ui_psi(tw, ttheta_all, tphi_all, tpsi_all)
spf_tb.show_ui_psi(tw2, ttheta_all, tphi_all, tpsi_all)
Out[390]:
In [393]:
tw2[0, 0, 0]
Out[393]:
In [313]:
# export frequence dates to mathematica
importlib.reload(spf_tb)
# ((psi, ((theta, phi, ui), )))
with open('%s.pickle' % fft_pickle, 'rb') as handle:
pickle_data = pickle.load(handle)
ttheta_all, tphi_all = pickle_data[0][1][0][:2]
tpsi_all = np.array([ti[0] for ti in pickle_data])
U_all = [[] for i in range(6)]
for _, table_psi_data in pickle_data:
for (_, _, tU), Ui in zip(table_psi_data, U_all):
Ui.append(tU)
for i0 in (3, 4, 5):
tw = np.dstack(U_all[i0])
tw_fft = np.fft.fftn(tw)
tname = '%s_w%d_fft.mat' % (fft_pickle, i0 - 2)
savemat(tname, {'tw_fft': np.moveaxis(tw_fft[:, :, :], 0, -1)})
print()
print('save fft data of %dth component to %s. ' % (i0, tname))
# print(tw_fft[0, 0, 0])
tw_fft[0, 0, 0] = 0
idx_max = np.dstack(np.unravel_index(np.argsort(np.abs(tw_fft).ravel()), tw_fft.shape))[0][::-1]
use_idx_list1 = idx_max[0::2]
use_idx_list2 = idx_max[1::2]
t1 = np.hstack([tw_fft[tk, tl, tj] for tk, tl, tj in use_idx_list1])
t2 = np.hstack([tw_fft[tk, tl, tj] for tk, tl, tj in use_idx_list2])
tidx = np.argwhere(np.logical_not(np.isclose(t1, np.conj(t2))))
tname = '%s_w%d_fft_idx.mat' % (fft_pickle, i0 - 2)
savemat(tname, {'use_idx_list': use_idx_list1})
print('save major index of fft data of %dth component to %s, unmatched case: ' % (i0, tname))
for i0 in tidx[:1]:
print('%5d: %3d, %3d, %3d, %15.5e, %15.5e' % (i0, use_idx_list1[i0, 0], use_idx_list1[i0, 1], use_idx_list1[i0, 2],
np.real(t1[i0]), np.imag(t1[i0])))
print(' %3d, %3d, %3d, %15.5e, %15.5e' % (use_idx_list2[i0, 0], use_idx_list2[i0, 1], use_idx_list2[i0, 2],
np.real(t2[i0]), np.imag(t2[i0])))
In [ ]:
# chech the conjunction pairs of FFT of rotation velocity.
idx_max = np.dstack(np.unravel_index(np.argsort(np.abs(tw_fft).ravel()), tw_fft.shape))[0][::-1]
use_idx_list1 = idx_max[0::2]
use_idx_list2 = idx_max[1::2]
t1 = np.hstack([tw_fft[tk, tl, tj] for tk, tl, tj in use_idx_list1])
t2 = np.hstack([tw_fft[tk, tl, tj] for tk, tl, tj in use_idx_list2])
tidx = np.argwhere(np.logical_not(np.isclose(t1, np.conj(t2))))
for i0 in range(10):
tk1, tl1, tj1 = use_idx_list1[i0]
tk2 = tM - tk1 if tk1 > 0 else tk1
tl2 = tN - tl1 if tl1 > 0 else tl1
tj2 = tO - tj1 if tj1 > 0 else tj1
print('%5d: %3d, %3d, %3d, %15.5e, %15.5e' % (i0, tk1, tl1, tj1,
np.real(t1[i0]), np.imag(t1[i0])))
print(' %3d, %3d, %3d, %15.5e, %15.5e' % (tk2, tl2, tj2,
np.real(t1[i0]), np.imag(t1[i0])))
print(' %3d, %3d, %3d, %15.5e, %15.5e' % (use_idx_list2[i0, 0], use_idx_list2[i0, 1], use_idx_list2[i0, 2],
np.real(t2[i0]), np.imag(t2[i0])))
In [405]:
from sympy.parsing import mathematica
import sympy
from sympy.printing.latex import LatexPrinter, print_latex
from sympy.utilities.lambdify import lambdify, lambdastr
import inspect
def mm2py(mm_string, sym_list, modules):
t12 = mathematica.mathematica(mm_string.replace('*^', '*10^'))
t13 = sympy.lambdify(sym_list, t12, modules=modules)
return t13
sym_list = (sympy.Symbol('theta'), sympy.Symbol('phi'), sympy.Symbol('psi'))
modules = ['numpy', ]
tw1_apx = '-0.0001866096570410253 + (-1.7029153048081684*^-6*Cos[62*phi] - 1960.88611819068*Sin[62*phi])/16384 + (279.32173923353105*Cos[2*phi + 2*psi] + 577.5209169432512*Sin[2*phi + 2*psi])/16384 + (-276.61411203034004*Cos[62*phi + 2*psi] - 571.7247600604117*Sin[62*phi + 2*psi])/16384 + (-6.114817144105495*Cos[2*theta] + 382.9878449974278*Sin[2*theta])/16384 + (-342.6737383563244*Cos[2*psi + 2*theta] + 166.50007489335445*Sin[2*psi + 2*theta])/16384 + (-326.49299534579876*Cos[2*phi + 2*psi + 2*theta] - 136.54707368952933*Sin[2*phi + 2*psi + 2*theta])/16384 + (-4.588689623494809*^-6*Cos[2*phi + 62*theta] - 1213.5905914763723*Sin[2*phi + 62*theta])/16384 + (0.000012727231708422586*Cos[62*phi + 62*theta] + 747.848109951239*Sin[62*phi + 62*theta])/16384 + (342.67373891716625*Cos[2*psi + 62*theta] - 166.5000738452281*Sin[2*psi + 62*theta])/16384 + (94.47535863228539*Cos[2*phi + 2*psi + 62*theta] - 339.70890817594295*Sin[2*phi + 2*psi + 62*theta])/16384 + (329.20060461853905*Cos[62*phi + 2*psi + 62*theta] + 142.34323812113675*Sin[62*phi + 2*psi + 62*theta])/16384 + (-91.76774684966742*Cos[2*phi + 14*psi + 62*theta] - 345.5050691031397*Sin[2*phi + 14*psi + 62*theta])/16384'
tw2_apx = '0.6197200974539849 + (-1962.0948612095603*Cos[62*phi] + 1.707528116724706*^-6*Sin[62*phi])/16384 + (1149.590231471158*Cos[14*psi] + 556.0920670954531*Sin[14*psi])/16384 + (-577.8673478225172*Cos[2*phi + 2*psi] + 279.47840108665457*Sin[2*phi + 2*psi])/16384 + (-572.0709983620775*Cos[62*phi + 2*psi] + 276.77041538443757*Sin[62*phi + 2*psi])/16384 + (5887.066706999647*Cos[62*theta] - 0.000026515346007727203*Sin[62*theta])/16384 + (136.62192184194265*Cos[2*phi + 2*psi + 2*theta] - 326.75582207414834*Sin[2*phi + 2*psi + 2*theta])/16384 + (-574.8788716391416*Cos[14*psi + 2*theta] - 278.08471606050534*Sin[14*psi + 2*theta])/16384 + (1213.9851268953007*Cos[2*phi + 62*theta] - 4.571595740233921*^-6*Sin[2*phi + 62*theta])/16384 + (748.0700251351213*Cos[62*phi + 62*theta] - 0.000012726440764286632*Sin[62*phi + 62*theta])/16384 + (142.4182788688823*Cos[62*phi + 2*psi + 62*theta] - 329.463789838122*Sin[62*phi + 2*psi + 62*theta])/16384 + (-574.8788879671743*Cos[14*psi + 62*theta] - 278.0847288505968*Sin[14*psi + 62*theta])/16384 + (345.6246033372179*Cos[2*phi + 14*psi + 62*theta] - 91.95740219001547*Sin[2*phi + 14*psi + 62*theta])/16384'
tw3_apx = '1.9074525687068844*^-16 + (518.7366380862838*Cos[63*phi] - 4.1392973060100844*^-6*Sin[63*phi])/16384 + (235.2128422040162*Cos[63*phi + 2*psi] + 486.31273548670134*Sin[63*phi + 2*psi])/16384 + (235.2128439364814*Cos[63*phi + 14*psi] - 486.3127504338104*Sin[63*phi + 14*psi])/16384 + (2099.122594069163*Cos[phi + 2*theta] + 5.5981210027771*^-6*Sin[phi + 2*theta])/16384 + (-496.40643828187825*Cos[phi + 2*psi + 2*theta] + 321.27367065113106*Sin[phi + 2*psi + 2*theta])/16384 + (560.2854005256127*Cos[63*phi + 2*psi + 2*theta] - 191.41794598717235*Sin[63*phi + 2*psi + 2*theta])/16384 + (-213.32024088220518*Cos[phi + 60*theta] - 5.69643947721086*^-6*Sin[phi + 60*theta])/16384 + (-213.45250720383183*Cos[63*phi + 60*theta] - 5.6899100784152866*^-6*Sin[63*phi + 60*theta])/16384 + (-79.60775091498382*Cos[63*phi + 14*psi + 60*theta] + 164.03928821590327*Sin[63*phi + 14*psi + 60*theta])/16384 + (-1827.2186834657691*Cos[phi + 62*theta] + 3.5884381293276077*^-6*Sin[phi + 62*theta])/16384 + (560.2853940215208*Cos[phi + 2*psi + 62*theta] - 191.41792073134854*Sin[phi + 2*psi + 62*theta])/16384 + (-496.4064319219047*Cos[63*phi + 2*psi + 62*theta] + 321.2736806311276*Sin[63*phi + 2*psi + 62*theta])/16384'
fun_w1_apx = lambda theta, phi, psi: mm2py(tw1_apx, sym_list, modules)(theta, phi, psi)
fun_w2_apx = lambda theta, phi, psi: mm2py(tw2_apx, sym_list, modules)(theta, phi, psi)
fun_w3_apx = lambda theta, phi, psi: mm2py(tw3_apx, sym_list, modules)(theta, phi, psi)
In [406]:
importlib.reload(spf_tb)
ttheta_all = np.linspace(0, np.pi, 33)
tphi_all = np.linspace(0, 2 * np.pi, 65)
tpsi_all = np.linspace(0, 2 * np.pi, 16, endpoint=False)
tw1 = np.dstack([fun_w1_apx(*np.meshgrid(ttheta_all, tphi_all, indexing='ij'), tpsi) for tpsi in tpsi_all])
tw2 = np.dstack([fun_w2_apx(*np.meshgrid(ttheta_all, tphi_all, indexing='ij'), tpsi) for tpsi in tpsi_all])
tw3 = np.dstack([fun_w3_apx(*np.meshgrid(ttheta_all, tphi_all, indexing='ij'), tpsi) for tpsi in tpsi_all])
In [402]:
spf_tb.show_ui_psi(tw1, ttheta_all, tphi_all, tpsi_all)
Out[402]:
In [403]:
spf_tb.show_ui_psi(tw2, ttheta_all, tphi_all, tpsi_all)
Out[403]:
In [404]:
spf_tb.show_ui_psi(tw3, ttheta_all, tphi_all, tpsi_all)
Out[404]:
In [419]:
sym_list = (sympy.Symbol('z'), sympy.Symbol('th0'), sympy.Symbol('ph0'), sympy.Symbol('ps0'),
sympy.Symbol('c01'), sympy.Symbol('c02'), sympy.Symbol('c03'))
modules = ['numpy', ]
theta_apx = 'th0 - 0.6197200974539849*c01*z*Cos[ph0] + 0.000036887909511729156*c01*z*Cos[61*ph0] + 0.11971987852173585*c01*z*Cos[63*ph0] - 0.1796590181579482*c01*z*Cos[ph0 - 62*th0] - 0.21669486384509337*c01*z*Cos[ph0 + 62*th0] + 0.03703584568714515*c01*z*Cos[3*ph0 + 62*th0] - 0.0001866096570410253*c01*z*Sin[ph0] - 1.4077184803886813*^-13*c01*z*Sin[61*ph0] - 1.0407847355752181*^-10*c01*z*Sin[63*ph0] - 8.091841433022218*^-10*c01*z*Sin[ph0 - 62*th0] + 9.492198373786016*^-10*c01*z*Sin[ph0 + 62*th0] - 1.4003569407637966*^-10*c01*z*Sin[3*ph0 + 62*th0]'
fun_theta_apx = lambda theta, phi, psi: mm2py(theta_apx, sym_list, modules)(theta, phi, psi)
In [423]:
# mm2py(theta_apx, sym_list, modules)
mm2py(theta_apx, sym_list, modules)(1, 1, 1, 1, -1, -1, -1)
Out[423]:
In [ ]:
In [635]:
importlib.reload(spf_tb)
use_uidx = 4 # 3: wx, 4: wy, 5: wz
with open('%s.pickle' % fft_pickle, 'rb') as handle:
pickle_data = pickle.load(handle)
ttheta_all, tphi_all = pickle_data[0][1][0][:2]
if tphi_all[-1] < (2 * np.pi):
tphi_all = np.hstack((tphi_all, 2 * np.pi))
if ttheta_all[-1] < (np.pi):
ttheta_all = np.hstack((ttheta_all, np.pi))
tpsi_all = np.array([ti[0] for ti in pickle_data])
U_all = [[] for i in range(6)]
for _, table_psi_data in pickle_data:
for (ttheta, tphi, tU), Ui in zip(table_psi_data, U_all):
if tphi[-1] < (2 * np.pi):
tU[2 * np.pi] = tU[0]
if ttheta[-1] < (np.pi):
tU = tU.append(tU.loc[0].rename(np.pi))
Ui.append(tU)
tw = np.dstack(U_all[use_uidx])
# spf_tb.show_ui_psi(tw, ttheta_all, tphi_all, tpsi_all)
In [725]:
th2, ph2 = np.meshgrid(ttheta_all, tphi_all, indexing='ij')
p2 = np.dstack((np.sin(th2) * np.cos(ph2), np.sin(th2) * np.sin(ph2), np.cos(th2)))
p21 = np.dstack((p2[..., 1], p2[..., 2], p2[..., 0]))
th21 = np.arccos(p21[..., 2])
ph21 = np.arctan2(p21[..., 1], p21[..., 0])
ph21[ph21 < 0] = ph21[ph21 < 0] + 2 * np.pi
%matplotlib inline
figsize = np.array((16, 9)) * 0.3
dpi = 100
fig, axs = plt.subplots(1, 2, figsize=figsize, dpi=dpi)
fig.patch.set_facecolor('white')
axi = axs[0]
axi.contourf(th2, ph2, th21, 50, cmap = plt.get_cmap('jet'))
axi = axs[1]
axi.contourf(th2, ph2, ph21, 50, cmap = plt.get_cmap('jet'))
fig.tight_layout()
# fig, axs = plt.subplots(1, 2, figsize=figsize, dpi=dpi)
# fig.patch.set_facecolor('white')
# axi = axs[0]
# axi.contourf(th2, ph2, tw[..., 0], 50, cmap = plt.get_cmap('jet'))
# axi = axs[1]
# axi.contourf(th21, ph21, tw[..., 0], 50, cmap = plt.get_cmap('jet'))
# fig.tight_layout()
In [707]:
np.hstack((p_map, np.array(p_map) + 3)) % 3
Out[707]:
In [724]:
# p_map = [0, 2, 1]
# p_sign = [1, -1, 1]
# rot_psi = np.pi / 2
# p_map = [0, 1, 2]
# p_sign = [1, 1, 1]
# rot_psi = 0
# p_map = [1, 2, 0]
# p_sign = [1, 1, 1]
# rot_psi = np.pi
# p_map = [1, 0, 2]
# p_sign = [-1, 1, 1]
# rot_psi = np.pi / 2
p_map = [2, 1, 0]
p_sign = [1, -1, 1]
rot_psi = np.pi
# p_map = [2, 0, 1]
# p_sign = [1, 1, 1]
# rot_psi = np.pi / 2
figsize = np.array((16, 9)) * 0.9
dpi = 100
th2, ph2 = np.meshgrid(ttheta_all, tphi_all, indexing='ij')
p2 = np.dstack((np.sin(th2) * np.cos(ph2), np.sin(th2) * np.sin(ph2), np.cos(th2)))
p21 = np.dstack((p_sign[0] * p2[..., p_map[0]], p_sign[1] * p2[..., p_map[1]], p_sign[2] * p2[..., p_map[2]]))
th21 = np.arccos(p21[..., 2])
ph21 = np.arctan2(p21[..., 1], p21[..., 0])
ph21[ph21 < 0] = ph21[ph21 < 0] + 2 * np.pi
# fig, axs = plt.subplots(1, 2, figsize=figsize, dpi=dpi)
# fig.patch.set_facecolor('white')
# axi = axs[0]
# axi.contourf(th2, ph2, th21, 50, cmap = plt.get_cmap('jet'))
# axi = axs[1]
# axi.contourf(th2, ph2, ph21, 50, cmap = plt.get_cmap('jet'))
# fig.tight_layout()
# tpsi_all2 = (tpsi_all + rot_psi) % (2 * np.pi)
# tidx = np.argsort(tpsi_all2).astype(int)
# tpsi_all2 = tpsi_all
tidx = np.hstack((p_map, np.array(p_map) + 3))
for use_uidx in tidx: # 3: wx, 4: wy, 5: wz
tw = p_sign[0] * p_sign[use_uidx % 3] * np.dstack(U_all[use_uidx])
tw2 = []
for i0 in range(tw.shape[-1]):
intp_fun = interpolate.RectBivariateSpline(ttheta_all, tphi_all, tw[..., i0], kx=3, ky=3)
tw2i = np.hstack([intp_fun(t1, t2) for t1, t2 in zip(th21.ravel(), ph21.ravel())]).reshape(th21.shape)
tw2.append(tw2i)
tw2 = np.dstack(tw2)
spf_tb.show_ui_psi(tw2, ttheta_all, tphi_all, tpsi_all)
# spf_tb.show_ui_psi(tw2, ttheta_all, tphi_all, tpsi_all, polar=True)