Quiet@Remove["Global`*"];
curl[vect_]:=Curl[vect, {r,θ,z}, "Cylindrical"];
Quiet@<<PhysicalConstants`

ϵ0=(N@VacuumPermittivity)[[1]];
μ0=(N@VacuumPermeability)[[1]];

Jvect={0,0,Jz[r]};
Hvect={0,Hfi[r],0};
eqn1=Thread[curl[Jvect]==-I*ω*γ*μ*Hvect][[2]];
eqn2=Thread[curl[Hvect]==Jvect][[3]];
sol=DSolve[{eqn1,eqn2},{Jz[r],Hfi[r]},r][[1]];
sol2=Simplify[sol,{γ>0,μ>0,ω>0}]/.C[2]->0;
jecko=Jz[r]/.sol2;
dosC1=current==Integrate[jecko*2*Pi*r,{r,0,R}]//Solve[#,C[1]]&;
Jfinal=jecko/.dosC1[[1]];

ClearAll[Qv,Pz];
dos={γ->57.*(10^6),μ->μ0,ω->2.*Pi*f,R->0.01,length->1,current->1};
Qv[f_]=(1/γ)Abs[(Jfinal^2)]/.dos;
Pz[f_]:=NIntegrate[Evaluate[2*Pi*r*length*Qv[f]/.dos],Evaluate[{r,0,R/.dos}]];
R1[f_]:=((Pz[f])/(current^2))/.dos;
dataR={#,R1[#]}&/@Range[50,1000,10];
ListPlot[dataR,AxesOrigin->{0,0}]

ClearAll[a];
a[f_]=(Sqrt[(2/(ω*γ*μ))])/.dos;
dat1=dataR/.{f_,R1_}:>{f,(R1/(length/(γ*Pi*(R^2))))/.dos};
dat2=dataR/.{f_,R1_}:>{f,(R1/(length/(γ*Pi*((R^2)-(((R-a[f]))^2)))))/.dos};
ListPlot[{dat1,dat2}]
ListPlot[dat2,PlotRange->All]

a[50]

Plot[Abs[Jfinal/.dos/.f->1050],{r,0,0.01},AxesOrigin->{0,0}]

Quiet@Remove["Global`*"];
f=50;
a=(Sqrt[(2./(2*Pi*f*γ*μ))]);
μ=4.Pi*(10^(-7));
γ=57*(10^6);
cm=0.01;
D2=20cm;
D1=D2+2cm;
d2=D2-a;
d1=D1+a;
S1=(Pi/4)*(d1^2);
S2=(Pi/4)*(d2^2);
length=1;
Print["d1/length=",(d1/length)]
Print["d2/length=",(d2/length)]
nagaoka=0.94;
L1=S1*μ*((N1^2)/length)*nagaoka;
L2=S2*μ*((1^2)/length)nagaoka;
F=0.97;
M=S2*μ*((N1*1)/length)*F;
R2=Pi*(d2/(γ*a*length));
R1=Pi*((1.2*d1*(N1^2))/(γ*a*length))

T1=20;
T2=1000;
V=Pi*((D2^2)/4)*length;
c=385;
density=ro=8960;
m=V*density;
energy=m*c*(T2-T1);
time=15*60;
P=(energy/time);

U1=120;
ω=50*2.Pi;
eqnPrimary=U1==R1*i1+I*ω*L1*i1+I*ω*M*i2;
eqnSecondary=0==R2*i2+I*ω*L2*i2+I*ω*M*i1;
eqnPower=P==R2*((Abs[i2])^2);
eqnPower/.Solve[eqnSecondary/.Solve[eqnPrimary,i1][[1]],i2][[1]];
solN1=Solve[%][[2]]
Solve[{eqnPrimary,eqnSecondary}/.solN1][[1]];
Ploses=R1*Abs[(i1^2)]/.%/.solN1;
Print["loses in prim. win.= ",(Ploses/1000),"kW"]
Print["efficiency= ",(P/(P+Ploses))]
Tw1=45;
Tw2=85;
Print["mass flow= ",mdot/.Solve[Ploses==mdot*4186*(Tw2-Tw1)][[1]]," kg/s"]