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}]