(*Nastavit souradnice a zadat vektory E a H*)
Remove["Global`*"]
Needs["VectorAnalysis`"];
SetCoordinates[Cylindrical[r,fi,z]];
Evector = {0,Efi[r],0};
Hvector = {0,0,Hz[r]};
(*Zadam Maxwellovy rovnice a podminky*)
dospars= {μ-> 8000*4Pi*(10^(-7)), ω -> 2*π*50,γ ->5.6*(10^6),ϵ->8.884*(10^(-9))};
(*Polomer valce, minimum nemuze byt nula kvuli deleni nulou*)
rmin = (10^(-4));
rmax=0.1;
delka=1;
(*Maxwellovy rovnice*)
equation1 =Thread[Curl[Hvector]==γ*Evector ][[2]];
equation2 =Thread[Curl[Evector]==-μ*I*ω*Hvector ][[3]];

(*Analyticke reseni Maxwellovych rovnic bez podminek*)
solve1=DSolve[{equation1,equation2},{Efi[r],Hz[r]},r];
(*Z okrajove podminky lze urcit C*)
(*Efi(rmin) -> 0 => C[2]=0*)
solve2=solve1/.C[2]->0;
(*Polozim Hz[rmax] = Hz z reseni Maxwellek*)
equation3=H0==Hz[r]/.solve2[[1]]/.r->rmax;

(*Vypoctena konstanta*)
c1=Solve[equation3,C[1]][[1]]/.dospars;
(*Vlozim C[1] do vysledku a dostanu vysledne E a H*)
finalsolve=solve2/.c1/.dospars;

(*Intenzita magnetickeho pole na okraji valce*)
magnetickepole=25000;
(*Dosadim do vysledku magneticke pole na okraji valce*)
ecko=Efi[r]/.finalsolve[[1]]/.H0->magnetickepole;
hacko=Hz[r]/.finalsolve[[1]]/.H0->magnetickepole;
(*Vynesene grafy E a H uvnitr valce*)
Plot[Abs[ecko],{r,rmin,rmax}, AxesLabel->{"r [m]","E [(Vm^(-1))]"},PlotRange->All]
Plot[Abs[hacko],{r,rmin,rmax},PlotRange->Full, AxesLabel->{"r [m]","H [(Am^(-1))]"}]

(*Vypoctu teplo vyvinute v jednotce objemu*)
QV=Re[γ*ecko]/.dospars;
(*Zadane materialove parametry valce*)
λ=50;
ρ=7800;
cp=469;

(*Soucinitel prestupu tepla*)
alfa=5;
(*Pocatecni teploty*)
Tvzduch=20;
Tsteny=Tvzduch;
(*Emisivita povrchu valce a Stefan-Boltzmannova konstanta*)
ϵ=0.8;
σ=5.67*(10^(-8));

(*Fourier-Kirchhoffova rovnice*)
FK=ρ*cp*D[T[r,t],t]==Div[λ*Grad[T[r,t]]]+Abs[QV];
(*Doba ohrevu*)
tmax=7200;
(*Podminky na okrajich valce, pocatecni teplota*)
okraj11=D[T[r,t],r]==0/.r->rmin;
okraj21=-λ*D[T[r,t],r]==Tanh[1t]*((alfa*(T[r,t]-Tvzduch)+ϵ*σ*((((T[r,t]+273))^4)-(((Tsteny+273))^4))))/.r->rmax;
pocatecniT=T[r,0]==Tvzduch;
(*Reseni Fourier-Kirchhoffovy rovnice*)
tecko=NDSolve[{FK,okraj11,okraj21,pocatecniT},T,{r,rmin,rmax},{t,0,tmax},MaxStepSize->{Abs[rmin-rmax]/500,Automatic}][[1]]

(*Vyneseny vysledek ve 3D grafu*)
Plot3D[T[r,t]/.tecko,{r,rmin,rmax},{t,0,tmax},PlotRange->All,AxesLabel->{"r [m]","t [s]","T [°C]"}]
(*Graf rozlozeni teploty uvnitr valce*)
Plot[T[r,tmax]/.tecko,{r,rmin,rmax},AxesLabel->{"r [m]","T [°C]"}]

(*Vypise teplotu na povrchu valce*)
teplotapovrch=T[rmax,tmax]/.tecko
(*Vypocet tepelneho vykonu ohrevu*)
vykon=Pi*((rmax^2)-(rmin^2))*delka*ρ*cp*(teplotapovrch-Tvzduch)*(1/tmax)