/*-- 
Chauffage par induction d'une billette intérieure à un tube
1A ENSEM -- Apprentissage -- Base de l'électromagnétisme
G. Vinsard
*/

//-- Coefficients et dimensions physiques
real ConductiviteBillette=1/(3*1.0e-6);            // S/m
real PermeabiliteBillette=1.0;                     // relative
real RayonInterieurBillette=0.02;                  // m
real RayonExterieurBillette=0.03;                  // m
real HauteurBillette=0.05;                         // m
real ConductiviteTube=1/(3*1.0e-6);                // S/m
real PermeabiliteTube=100.0;                         // relative
real RayonInterieurTube=0.034;                     // m
real RayonExterieurTube=0.035;                     // m
real HauteurTube=0.08;                             // m
real RayonSpires=0.073/2;                          // m
real RayonSectionSpires=0.001/2;                   // m
real EspaceEntreSpires=0.001;                      // m 
real CourantInducteur=1000;                        // A
real frequence=5000;                               // Hz

//-- Parametre numerique qui peut etre ajuste
real DensiteMaillage=10000;                          // m^-1


//-- Parametres numeriques
real Lb=HauteurBillette/2.0;
real Lt=HauteurTube/2.0;
real Rti=RayonInterieurTube;
real Rtx=RayonExterieurTube;
real e=EspaceEntreSpires;
real a=RayonSectionSpires;
real Li=4*(2*a+e);
real LL=max(Li,Lb);
real Ri=RayonSpires;
real RayonExterieur=2*sqrt(Ri*Ri+LL*LL);            // m    

//-- Constantes
real mu0=4*3.141592653589793*1.0e-7;                // H/m permeabilite du vide
real ConductiviteCuivre=1/(2.0*1.0e-8);             // S/m

//-- Abreviations
real sigma=ConductiviteBillette;
real mur=PermeabiliteBillette;
real sigmat=ConductiviteTube;
real murt=PermeabiliteTube;
real R=RayonExterieurBillette;
real Rbi=RayonInterieurBillette;
real Rs=RayonSpires;
real f=frequence;
real I=CourantInducteur;
real Rinf=RayonExterieur;
real dens=DensiteMaillage;

//-- Coefficients intermédiaires
real omega=2*pi*f;
complex cf=1i*sigma*omega*mu0;
complex cft=1i*sigmat*omega*mu0;
complex czero=0;
complex cjs=I*mu0/(pi*a^2);

//-- Maillage
border b1(t=0,1){x=Rbi*t;y=0;label=1;};  dens=dens/10;int l1=max(dens*Rbi,2.0);dens=dens*10;
border b2(t=0,1){x=R+(Rti-R)*t;y=0;label=1;};  int l2=max(dens*(Rti-R),2.0);
border b7(t=0,1){x=Rti+(Rtx-Rti)*t;y=0;label=1;};  int l7=max(dens*(Rtx-Rti),2.0);
border b10(t=0,1){x=Rtx;y=Lt*t;label=4;}; int l10=max(dens*Lt,2.0);
border b13(t=0,1){x=Rti;y=Lt*t;label=4;}; int l13=max(dens*Lt,2.0);
border b12(t=0,1){x=Rti+(Rtx-Rti)*t;y=Lt;label=1;};  int l12=max(dens*(Rtx-Rti),2.0);
real Rinf0=Rs+(Rtx-Rbi);
border b3(t=0,1){x=Rtx+(Rinf0-Rtx)*t;y=0;label=1;};int l3=max(dens*(Rinf0-Rti),2.0);
border b17(t=0,1){x=Rinf0+(Rinf-Rinf0)*t;y=0;label=1;};dens=dens/10; int l17=max(dens*(Rinf-Rinf0),2.0);dens=dens*10;
border b5(t=0,1){x=0;y=t*Lb;label=2;}; dens=dens/10;int l5=max(dens*Lb,2.0);dens=dens*10;
border b6(t=0,1){x=0;y=Lb+(Rinf-Lb)*t;label=2;};dens=dens/10;int l6=max(dens*(Rinf-Lb),2.0);dens=dens*10;
border b8(t=0,1){x=R;y=Lb*t;label=3;}; int l8=max(dens*Lb,2.0);
//border b9(t=0,1){x=Rbi*t;y=Lb;label=3;}; int l9=l1;
border b14(t=0,1){x=Rbi+t*(R-Rbi);y=Lb;label=3;}; int l14=max(dens*(R-Rbi),2.0);
border b15(t=0,1){x=Rbi+t*(R-Rbi);y=0;label=3;}; int l15=l14;
border b16(t=0,1){x=Rbi;y=Lb*t;label=3;}; int l16=max(dens*Lb,2.0);
border b21(t=0,2*pi){x=Rs+a*cos(t);y=e/2+a+0*(e+2*a)+a*sin(t);label=6;};int l21=max(dens*2*pi*a,3.0);
border b22(t=0,2*pi){x=Rs+a*cos(t);y=e/2+a+1*(e+2*a)+a*sin(t);label=7;};int l22=max(dens*2*pi*a,3.0);
border b23(t=0,2*pi){x=Rs+a*cos(t);y=e/2+a+2*(e+2*a)+a*sin(t);label=8;};int l23=max(dens*2*pi*a,3.0);
border b24(t=0,2*pi){x=Rs+a*cos(t);y=e/2+a+3*(e+2*a)+a*sin(t);label=9;};int l24=max(dens*2*pi*a,3.0);
dens=dens/10.0;
border b11(t=0,pi/2){x=Rinf*cos(t);y=Rinf*sin(t);label=5;};int l11=max(dens*pi/2*Rinf,2.0);
dens=dens*10.0;
//plot(b1(l1)+b2(l2)+b3(l3)+b5(-l5)+b6(-l6)+b8(l8)+b9(-l9)+b11(l11)+b7(l7)+b10(l10)+b12(l12)+b13(l13));
mesh m=buildmesh(b1(l1)+b2(l2)+b3(l3)+b5(-l5)+b6(-l6)+b8(l8)+b11(l11)+b7(l7)//+b9(-l9)
+b14(l14)+b15(l15)+b16(l16)+b17(l17)
+b10(l10)+b12(l12)+b13(l13)
+b21(l21)+b22(l22)+b23(l23)+b24(l24));
plot(m);

//-- Indicateurs
fespace mP0(m,P0);mP0 reg=region; // Rituel ad hoc pour définir la fonction « reg » qui renvoie le label du sous-domaine dans lequel est le point dont les coordonnées sont ses arguments
int IndicateurBillette=reg((R+Rbi/2)/2.0,Lb/2.0);
int IndicateurInducteur=4;
int IndicateurAir=reg(R/2.0,(Lb+Rinf)/2.0);
int Ib=IndicateurBillette;
int IndicateurTube=reg((Rti+Rtx)/2.0,Lt/2.0);
int It=IndicateurTube;
int Ii1=reg(Rs,e/2+a+0*(e+2*a));
int Ii2=reg(Rs,e/2+a+1*(e+2*a));
int Ii3=reg(Rs,e/2+a+2*(e+2*a));
int Ii4=reg(Rs,e/2+a+3*(e+2*a));
int Ia=IndicateurAir;
//cout << "Indicateurs : " << "Billette = " << Ib << " ; inducteur = " << Ii1 << ", " << Ii2 << ", " << Ii3 << ", " << Ii4   << " ; air =      " << Ia << endl;

real delta,deltat;
cout << endl << "Dimension Triangle = " << 1000/dens << " mm ; Nb noeuds = " << m.nv << " ; Nb Triangles = " << m.nt << endl;
if (sigma > 0) {delta=sqrt(2.0/(omega*sigma*mu0*mur));
cout << "delta Billette = " << 1000*delta << " mm ; rayon Billette = " << 1000*R << " mm" << endl;};
if (sigmat > 0) {deltat=sqrt(2.0/(omega*sigmat*mu0*murt));
cout << "delta Tube = " << 1000*deltat << " mm ; Epaisseur Tube = " << 1000*(Rtx-Rti) << " mm" << endl;};

//-- Espaces fonctionnels
fespace mP1(m,P1);
mP1<complex> A,Ap;

//-- Calcul 
solve Calc(A,Ap)=
  int2d(m,Ia,Ii1,Ii2,Ii3,Ii4)((x*dy(A)*dy(Ap)+x*dx(A)*dx(Ap)+dx(A)*Ap+dx(Ap)*A+1/x*Ap*A))
 +int2d(m,Ib)((x*dy(A)*dy(Ap)+x*dx(A)*dx(Ap)+dx(A)*Ap+dx(Ap)*A+1/x*Ap*A)/mur)
 +int2d(m,Ib)(x*A*Ap*cf)
 +int2d(m,It)((x*dy(A)*dy(Ap)+x*dx(A)*dx(Ap)+dx(A)*Ap+dx(Ap)*A+1/x*Ap*A)/murt)
 +int2d(m,It)(x*A*Ap*cft)
 +int2d(m,Ii1)(x*Ap*cjs)
 +int2d(m,Ii2)(x*Ap*cjs)
 +int2d(m,Ii3)(x*Ap*cjs)
 +int2d(m,Ii4)(x*Ap*cjs)
 +on(2,5,A=czero);

//-- Exploitation
mP1 Ar=real(A), Ai=imag(A), rAr=x*Ar, rAi=x*Ai;
mP1<complex>  Ac=conj(A);
plot(rAr,cmm="r Re(a)");
plot(rAi,cmm="r Im(a)");

//complex ctest=int2d(m,Ii1)(A);cout << ctest << endl;

complex cS=-int2d(m,Ii1,Ii2,Ii3,Ii4)(x*A)*cjs*2*pi/mu0*1i*omega*2;
complex cPj=int2d(m,Ib)(x*A*Ac)*2*pi*sigma*omega*omega*2;
complex cPjt=int2d(m,It)(x*A*Ac)*2*pi*sigmat*omega*omega*2;
complex cQ=
  int2d(m,Ib)((x*dy(A)*dy(Ac)+x*dx(A)*dx(Ac)+dx(A)*Ac+dx(Ac)*A+1/x*Ac*A)/mur)
+  int2d(m,It)((x*dy(A)*dy(Ac)+x*dx(A)*dx(Ac)+dx(A)*Ac+dx(Ac)*A+1/x*Ac*A)/murt)
+ int2d(m,Ia,Ii1,Ii2,Ii3,Ii4)((x*dy(A)*dy(Ac)+x*dx(A)*dx(Ac)+dx(A)*Ac+dx(Ac)*A+1/x*Ac*A));
cQ=1i*omega*cQ/mu0*2*pi*2;

cout << "S = " << real(cS) << " + " << imag(cS) << " i "   << " ; Pj Billette = " << real(cPj) << " ; Pj Tube = " << real(cPjt)  << " ; Q = " << imag(cQ) << endl ;
complex Z=cS/I^2;
cout << "Z = " << real(Z) << " + " << imag(Z) << " i " << " (Ohm) " << endl;
cout << "1-Resistance due aux charges = " << 1000*real(Z) << " milli-Ohm "  << endl;
cout << "2-Inductance = " << 1.0e+6*imag(Z)/omega << " micro-Henry" << endl;
real Rind=1000*8*2*pi*Ri/(pi*a^2)/ConductiviteCuivre;
cout << "3-Estimation Resistance Inducteur = " << Rind << " milli-Ohm" << endl; 
cout << "4-Resistance 1 + 3 = " << Rind+1000*real(Z) << " milli-Ohm" << endl;

//cout << endl << "ok" << endl;