From 951769f8e62e4991f49243c9aa67efad3c4b69a6 Mon Sep 17 00:00:00 2001 From: Guillaume Demesy <guillaume.demesy@fresnel.fr> Date: Wed, 4 Sep 2019 17:44:46 +0200 Subject: [PATCH] check - python and matplotlib problem --- ElectromagneticScattering/scattering.pro | 5 +- ElectromagneticScattering/scattering_init.py | 85 +++++++++++--------- 2 files changed, 48 insertions(+), 42 deletions(-) diff --git a/ElectromagneticScattering/scattering.pro b/ElectromagneticScattering/scattering.pro index ecc97ab..9c9107d 100644 --- a/ElectromagneticScattering/scattering.pro +++ b/ElectromagneticScattering/scattering.pro @@ -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 } } diff --git a/ElectromagneticScattering/scattering_init.py b/ElectromagneticScattering/scattering_init.py index 56609b5..b09a0ab 100644 --- a/ElectromagneticScattering/scattering_init.py +++ b/ElectromagneticScattering/scattering_init.py @@ -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') -- GitLab