diff --git a/ElectromagneticScattering/scattering.geo b/ElectromagneticScattering/scattering.geo
index 4d4851c0699c997effebf4bdb2e60063483b1526..ab6fdb629d90a4ddbb64f6d7f513c33e9adeb693 100644
--- a/ElectromagneticScattering/scattering.geo
+++ b/ElectromagneticScattering/scattering.geo
@@ -723,6 +723,7 @@ If(flag_cartpml==0)
Physical Volume(2) = {1073}; // Out - air
Physical Volume(3) = {1091}; // pml
Physical Point(2000) = {1}; // PrintPoint
+ Physical Surface(10) = {1055,1057,1059,1062,1064,1066,1068,1070};
EndIf
// Mesh.Algorithm = 1; // // 1=MeshAdapt, 5=Delaunay, 6=Frontal
@@ -752,3 +753,4 @@ Printf(Sprintf("n_max = %g;" ,n_max)) >> "scattering_tmp.py";
Printf(Sprintf("p_max = %g;" ,p_max)) >> "scattering_tmp.py";
Printf(Sprintf("siwt = %g;" ,siwt)) >> "scattering_tmp.py";
+Mesh.ElementOrder = 2;
diff --git a/ElectromagneticScattering/scattering.pro b/ElectromagneticScattering/scattering.pro
index 309c4e43e640a6bbe2b126127075aaea498faeec..8768a0c1d6ba0932935c3c89b4f333b9a693ae0d 100644
--- a/ElectromagneticScattering/scattering.pro
+++ b/ElectromagneticScattering/scattering.pro
@@ -21,6 +21,7 @@ Group {
Domain = Region[{Scat_In,Scat_Out}];
PMLs = Region[{PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
All_domains = Region[{Scat_In,Scat_Out,PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
+ SurfInt = Region[{}];
EndIf
If(flag_cartpml==0)
Scat_In = Region[1];
@@ -28,6 +29,7 @@ Group {
Domain = Region[{Scat_In,Scat_Out}];
PMLs = Region[3];
All_domains = Region[{Scat_In,Scat_Out,PMLs}];
+ SurfInt = Region[{10}];
EndIf
PrintPoint = Region[2000];
}
@@ -37,6 +39,7 @@ Function{
// For i In {0:#test()-1}
// Printf("%g",test(i));
// EndFor
+
I[] = Complex[0,1];
avoid_sing = 1.e-14;
mu0 = 4*Pi*100.0*nm;
@@ -192,6 +195,8 @@ Function{
me = ne*(ne+1) - Floor[pe];
Mnm_source~{pe}[] = Mnm[1,ne,me,XYZ[],k_Out];
Nnm_source~{pe}[] = Nnm[1,ne,me,XYZ[],k_Out];
+ Mnm_out~{pe}[] = Mnm[4,ne,me,XYZ[],k_Out];
+ Nnm_out~{pe}[] = Nnm[4,ne,me,XYZ[],k_Out];
source_M~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Mnm_source~{pe}[];
source_N~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm_source~{pe}[];
EndFor
@@ -235,10 +240,10 @@ Integration {
Case {
{ Type Gauss ;
Case {
- { GeoElement Point ; NumberOfPoints 1 ; }
- { GeoElement Line ; NumberOfPoints 4 ; }
- { GeoElement Triangle ; NumberOfPoints 6 ; }
- { GeoElement Tetrahedron ; NumberOfPoints 15 ; }
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line2 ; NumberOfPoints 4 ; }
+ { GeoElement Triangle2 ; NumberOfPoints 12 ; }
+ { GeoElement Tetrahedron2 ; NumberOfPoints 17 ; }
}
}
}
@@ -249,16 +254,16 @@ FunctionSpace {
{ Name Hcurl; Type Form1;
BasisFunction {
{ Name sn; NameOfCoef un; Function BF_Edge;
- Support Region[All_domains]; Entity EdgesOf[All]; }
+ Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
{ Name sn2; NameOfCoef un2; Function BF_Edge_2E;
- Support Region[All_domains]; Entity EdgesOf[All]; }
+ Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
If (is_FEM_o2==1)
{ Name sn3; NameOfCoef un3; Function BF_Edge_3F_b;
- Support Region[All_domains]; Entity FacetsOf[All]; }
+ Support Region[{All_domains,SurfInt}]; Entity FacetsOf[All]; }
{ Name sn4; NameOfCoef un4; Function BF_Edge_3F_c;
- Support Region[All_domains]; Entity FacetsOf[All]; }
+ Support Region[{All_domains,SurfInt}]; Entity FacetsOf[All]; }
{ Name sn5; NameOfCoef un5; Function BF_Edge_4E;
- Support Region[All_domains]; Entity EdgesOf[All]; }
+ Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
EndIf
}
Constraint {
@@ -316,6 +321,29 @@ Formulation {
}
}
EndFor
+ {Name VPWMN_helmholtz_vector; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [ (omega0/cel)^2*(epsilonr[]-epsilonr1[])*
+ ($isN ? Nnm[1,$NE,$ME,XYZ[],k_Out] : Mnm[1,$NE,$ME,XYZ[],k_Out])
+ , {u} ]; In Scat_In; Jacobian JVol; Integration Int_1;}
+ }
+ }
+ {Name VPWN_helmholtz_vector_test; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [ (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm[1,$NE,$ME,XYZ[],k_Out], {u} ]; In Scat_In; Jacobian JVol; Integration Int_1;}
+ }
+ }
+
EndIf
If (flag_study==RES_GREEN)
For ncomp In {0:2}
@@ -350,27 +378,35 @@ Resolution {
}
EndIf
If (flag_study==RES_TMAT)
- { Name res_VPWall_helmholtz_vector;
+ { Name res_VPWall_helmholtz_vector;
System {
- For pe In {1:p_max}
- { Name M~{pe}; NameOfFormulation VPWM_helmholtz_vector~{pe}; Type ComplexValue; Frequency Freq; }
- { Name N~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe}; Type ComplexValue; Frequency Freq; }
- EndFor
+ { Name T; NameOfFormulation VPWMN_helmholtz_vector; Type ComplexValue; }
}
Operation {
CreateDir[Str[myDir]];
Evaluate[Python[]{"scattering_init.py"}];
+
For pe In {1:p_max}
- Generate[M~{pe}];
- Solve[M~{pe}];
+ Evaluate[ $isN = 0 ];
+ Evaluate[ $PE = pe ];
+ Evaluate[ $NE = Floor[Sqrt[$PE]] ];
+ Evaluate[ $ME = $NE*($NE+1) - Floor[$PE] ];
+ If (pe==1)
+ Generate[T];
+ Solve[T];
+ EndIf
+ GenerateRHS[T];
+ SolveAgain[T];
PostOperation[VPWM_postop~{pe}];
- Generate[N~{pe}];
- Solve[N~{pe}];
+ Evaluate[$isN=1];
+ GenerateRHS[T];
+ SolveAgain[T];
PostOperation[VPWN_postop~{pe}];
EndFor
Evaluate[Python[]{"scattering_post.py"}];
}
}
+
EndIf
If (flag_study==RES_GREEN)
{ Name res_GreenAll_helmholtz_vector;
@@ -421,7 +457,7 @@ PostProcessing {
EndIf
If (flag_study==RES_PW)
- { Name PW_postpro ; NameOfFormulation PW_helmholtz_vector;NameOfSystem P;
+ { Name PW_postpro ; NameOfFormulation PW_helmholtz_vector; NameOfSystem P;
Quantity {
{ Name E_tot ; Value { Local { [{u}+Einc[]]; In All_domains; Jacobian JVol; } } }
{ Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
@@ -439,9 +475,10 @@ PostProcessing {
}
}
EndIf
+
If (flag_study==RES_TMAT)
For pe In {1:p_max}
- { Name VPWM_postpro~{pe}; NameOfFormulation VPWM_helmholtz_vector~{pe};NameOfSystem M~{pe};
+ { Name VPWM_postpro~{pe}; NameOfFormulation VPWMN_helmholtz_vector; NameOfSystem T;
Quantity {
{ Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
{ Name E_scat_sph ; Value { Local { [Vector[
@@ -452,9 +489,13 @@ PostProcessing {
In All_domains; Jacobian JVol; } } }
{ 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~{po}[]] ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+ { Name norma~{pe}~{po} ; Value { Integral { [ Mnm_out~{pe}[]*Conj[Mnm_out~{po}[]] ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+ EndFor
}
}
- { Name VPWN_postpro~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe};NameOfSystem N~{pe};
+ { Name VPWN_postpro~{pe}; NameOfFormulation VPWMN_helmholtz_vector; NameOfSystem T;
Quantity {
{ Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
{ Name E_scat_sph ; Value { Local { [Vector[
@@ -469,6 +510,7 @@ PostProcessing {
}
EndFor
EndIf
+
If (flag_study==RES_GREEN)
For ncomp In {0:2}
{ Name GreenAll_postpro~{ncomp}; NameOfFormulation GreenAll_helmholtz_vector~{ncomp};NameOfSystem M~{ncomp};
@@ -520,6 +562,7 @@ PostOperation {
}
}
EndIf
+
If (flag_study==RES_PW)
{Name PW_postop; NameOfPostProcessing PW_postpro ;
Operation {
@@ -540,6 +583,7 @@ PostOperation {
}
}
EndIf
+
If (flag_study==RES_TMAT)
For pe In {1:p_max}
{Name VPWM_postop~{pe}; NameOfPostProcessing VPWM_postpro~{pe} ;
@@ -557,6 +601,10 @@ PostOperation {
{sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} },
File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]],Sprintf["_r%g.pos",kr]]],
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
}
}
{Name VPWN_postop~{pe}; NameOfPostProcessing VPWN_postpro~{pe} ;
@@ -578,6 +626,7 @@ PostOperation {
}
EndFor
EndIf
+
If (flag_study==RES_GREEN)
For ncomp In {0:2}
{Name GreenAll_postop~{ncomp}; NameOfPostProcessing GreenAll_postpro~{ncomp} ;
@@ -656,4 +705,4 @@ If (flag_study==RES_GREEN)
R_ = {"res_GreenAll_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
-EndIf
\ No newline at end of file
+EndIf
diff --git a/ElectromagneticScattering/scattering_init.py b/ElectromagneticScattering/scattering_init.py
index 2f7621c6a6b329b909510f8b92f14a5f33529dcf..51761f3fdbebc097872ebb94cca385f5e6920df5 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
@@ -38,7 +38,7 @@ def field_VSH_expansion(post_filename):
FF_erCrossXnm_p = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
#####################################
##### sn, pn ,un
- ##### Brian's recurrence relations
+ ##### Brian Stout's recurrence relations (B_ arrays)
B_Pnpms = np.zeros((npts_theta,m_max+1,n_max+1))
B_unpms = np.zeros((npts_theta,m_max+1,n_max+1))
B_snpms = np.zeros((npts_theta,m_max+1,n_max+1))
@@ -88,7 +88,7 @@ def field_VSH_expansion(post_filename):
#####################################
##### sn, pn ,un
- ##### Brian's recurrence relations
+ ##### Brian Stout's recurrence relations (B_ arrays)
m_max = n_max
p_max = n_max*n_max+2*n_max
aM_nm = np.zeros(p_max,dtype=complex)
@@ -186,15 +186,15 @@ def field_VSH_expansion(post_filename):
B_Nnm_t = 1./(k_Out*r_sph) * dRicatti_dx_ofkr * B_Znm_t
B_Nnm_p = 1./(k_Out*r_sph) * dRicatti_dx_ofkr * B_Znm_p
- B_EdotconjYnm = E_scat_onsphere_sph_r*B_Ynm_r.conjugate()
- B_EdotconjZnm = E_scat_onsphere_sph_t*B_Znm_t.conjugate() + E_scat_onsphere_sph_p*B_Znm_p.conjugate()
- B_EdotconjXnm = E_scat_onsphere_sph_t*B_Xnm_t.conjugate() + E_scat_onsphere_sph_p*B_Xnm_p.conjugate()
- normalize_fhnm_X = 1./sph_bessel_n_ofkr
- normalize_fenm_Y = k_Out*r_sph/(sph_bessel_n_ofkr*np.sqrt(float(n)*(float(n)+1.)) )
- normalize_fenm_Z = k_Out*r_sph/dRicatti_dx_ofkr
- fenm_Y[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjYnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fenm_Y
- fenm_Z[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjZnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fenm_Z
- fhnm_X[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fhnm_X
+ B_EdotconjYnm = E_scat_onsphere_sph_r*B_Ynm_r.conjugate()
+ B_EdotconjZnm = E_scat_onsphere_sph_t*B_Znm_t.conjugate() + E_scat_onsphere_sph_p*B_Znm_p.conjugate()
+ B_EdotconjXnm = E_scat_onsphere_sph_t*B_Xnm_t.conjugate() + E_scat_onsphere_sph_p*B_Xnm_p.conjugate()
+ normalize_fhnm_X = 1./sph_bessel_n_ofkr
+ normalize_fenm_Y = k_Out*r_sph/(sph_bessel_n_ofkr*np.sqrt(float(n)*(float(n)+1.)) )
+ normalize_fenm_Z = k_Out*r_sph/dRicatti_dx_ofkr
+ fenm_Y[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjYnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fenm_Y
+ fenm_Z[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjZnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fenm_Z
+ fhnm_X[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fhnm_X
EdotconjMnm = E_scat_onsphere_sph_r*B_Mnm_r.conjugate() + E_scat_onsphere_sph_t*B_Mnm_t.conjugate() + E_scat_onsphere_sph_p*B_Mnm_p.conjugate()
EdotconjNnm = E_scat_onsphere_sph_r*B_Nnm_r.conjugate() + E_scat_onsphere_sph_t*B_Nnm_t.conjugate() + E_scat_onsphere_sph_p*B_Nnm_p.conjugate()
@@ -218,41 +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('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')