Skip to content
Snippets Groups Projects
Commit 951769f8 authored by Guillaume Demesy's avatar Guillaume Demesy
Browse files

check - python and matplotlib problem

parent 2553061f
No related branches found
No related tags found
No related merge requests found
Pipeline #4710 passed
......@@ -531,8 +531,8 @@ PostProcessing {
{ Name H_scat ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
{ Name Mnm_source~{pe} ; Value { Local { [ Mnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
For po In {1:p_max}
{ Name a~{pe}~{po} ; Value { Integral { [ {u}*Conj[Mnm_out~{pe}[]]]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
{ Name norma~{pe}~{po} ; Value { Integral { [Mnm_out~{pe}[]*Conj[Mnm_out~{pe}[]]]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
{ Name a~{pe}~{po} ; Value { Integral { [ (Cart2Sph[XYZ[]]*{u})*(Cart2Sph[XYZ[]]*Conj[Mnm_out~{pe}[]])]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
{ Name norma~{pe}~{po} ; Value { Integral { [(Cart2Sph[XYZ[]]*Mnm_out~{pe}[])*(Cart2Sph[XYZ[]]*Conj[Mnm_out~{pe}[]])]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
EndFor
}
}
......@@ -654,6 +654,7 @@ PostOperation {
Name StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]]];
For po In {1:p_max}
Print[a~{pe}~{po}[SurfInt], OnGlobal , Format Table, File StrCat[myDir,StrCat[StrCat["a_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
Print[norma~{pe}~{po}[SurfInt], OnGlobal , Format Table, File StrCat[myDir,StrCat[StrCat["norma_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
EndFor
}
}
......
......@@ -11,8 +11,8 @@ import sys,os
import numpy as np
from scipy.special import jv, yv, hankel1, hankel2
sys.path.append(os.getcwd())
from matplotlib import cm
import pylab as pl
# from matplotlib import cm
# import pylab as pl
from scattering_tmp import *
np.set_printoptions(precision=2)
pi = np.pi
......@@ -218,43 +218,48 @@ def field_VSH_expansion(post_filename):
orth3 = np.trapz(np.trapz((np.sin(theta_sph)*ZnmconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
aM_nm[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*EdotconjMnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) / normalize_aM_nm2
bN_nm[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*EdotconjNnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) / normalize_bN_nm2
print(np.loadtxt('run_results/a_pe%gpo%g.dat'%(1,int(po))))
print(aM_nm)
print('anm not normalized' , np.trapz(np.trapz((np.sin(theta_sph)*EdotconjMnm).transpose(),theta_sph[:,0]),phi_sph[0,:]))
print('normalization ' , normalize_aM_nm2)
# print(np.loadtxt('run_results/a_pe%gpo%g.dat'%(1,int(po))))
# print(aM_nm)
return [fhnm_X,fenm_Y,fenm_Z,aM_nm,bN_nm,FF_Xnm_t,FF_Xnm_p,FF_erCrossXnm_t,FF_erCrossXnm_p]
###########################
def plot_farfield(far_field_sph,filename):
vmax_sph = np.max(far_field_sph)
vmin_sph = np.min(far_field_sph)
x_sph = far_field_sph/vmax_sph*np.sin(theta_sph) * np.cos(phi_sph)
y_sph = far_field_sph/vmax_sph*np.sin(theta_sph) * np.sin(phi_sph)
z_sph = far_field_sph/vmax_sph*np.cos(theta_sph)
fig = pl.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_aspect('equal')
surf=ax.plot_surface(x_sph,y_sph,z_sph,\
facecolors=cm.viridis(far_field_sph/vmax_sph),\
rstride=2,\
cstride=2,\
linewidth=0,\
vmin = vmin_sph,\
vmax = vmax_sph,\
shade=True,\
alpha=0.5,\
antialiased=False)
cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='z', offset=-1)
cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='x', offset=-1)
cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='y', offset=1 )
surf.set_edgecolor('k')
max_range = 0.5*np.max(np.array([x_sph.max()-x_sph.min(), y_sph.max()-y_sph.min(), z_sph.max()-z_sph.min()]))
mid_x = (x_sph.max()+x_sph.min())*0.5
mid_y = (y_sph.max()+y_sph.min())*0.5
mid_z = (z_sph.max()+z_sph.min())*0.5
ax.set_xlim(mid_x-max_range, mid_x+max_range)
ax.set_ylim(mid_y-max_range, mid_y+max_range)
ax.set_zlim(mid_z-max_range, mid_z+max_range)
ax.xaxis.set_ticklabels([]);ax.xaxis.set_label('x')
ax.yaxis.set_ticklabels([]);ax.yaxis.set_label('y')
ax.zaxis.set_ticklabels([]);ax.zaxis.set_label('z')
pl.savefig(filename,bbox_inches='tight')
pl.close('all')
# # ###########################
# # def plot_farfield(far_field_sph,filename):
# # vmax_sph = np.max(far_field_sph)
# # vmin_sph = np.min(far_field_sph)
# # x_sph = far_field_sph/vmax_sph*np.sin(theta_sph) * np.cos(phi_sph)
# # y_sph = far_field_sph/vmax_sph*np.sin(theta_sph) * np.sin(phi_sph)
# # z_sph = far_field_sph/vmax_sph*np.cos(theta_sph)
# # fig = pl.figure()
# # ax = fig.add_subplot(111, projection='3d')
# # ax.set_aspect('equal')
# # surf=ax.plot_surface(x_sph,y_sph,z_sph,\
# # facecolors=cm.viridis(far_field_sph/vmax_sph),\
# # rstride=2,\
# # cstride=2,\
# # linewidth=0,\
# # vmin = vmin_sph,\
# # vmax = vmax_sph,\
# # shade=True,\
# # alpha=0.5,\
# # antialiased=False)
# # cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='z', offset=-1)
# # cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='x', offset=-1)
# # cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='y', offset=1 )
# # surf.set_edgecolor('k')
# # max_range = 0.5*np.max(np.array([x_sph.max()-x_sph.min(), y_sph.max()-y_sph.min(), z_sph.max()-z_sph.min()]))
# # mid_x = (x_sph.max()+x_sph.min())*0.5
# # mid_y = (y_sph.max()+y_sph.min())*0.5
# # mid_z = (z_sph.max()+z_sph.min())*0.5
# # ax.set_xlim(mid_x-max_range, mid_x+max_range)
# # ax.set_ylim(mid_y-max_range, mid_y+max_range)
# # ax.set_zlim(mid_z-max_range, mid_z+max_range)
# # ax.xaxis.set_ticklabels([]);ax.xaxis.set_label('x')
# # ax.yaxis.set_ticklabels([]);ax.yaxis.set_label('y')
# # ax.zaxis.set_ticklabels([]);ax.zaxis.set_label('z')
# # pl.savefig(filename,bbox_inches='tight')
# # pl.close('all')
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment