E. Métral
17 June 2005

TRANSVERSE RESISTIVE-WALL IMPEDANCE FOR A CIRCULAR BEAM PIPE :
1) General multi-Layer formula from Burov-Lebedev (BL) implemented by Alex Koschik
2) 1- and 2-layer formulae derived by me (from Zotter's formalism) and implemented by me

The numbers in red have to be entered

Remarks: 
1) BL's formula is valid for (β = 1) only, whereas Zotter's formulae are valid for any β 
2) BL's formula is not valid at too high frequencies, whereas Zotter's formulae are
3) The 2-layer formula derived from Zotter's formalism also allows the possibility to have a wave velocity different from the beam velocity (multi-bunch case)

Off[General::spell];
Off[General::spell1];

1) General Constants

c = 2.997925*10^8;
e = 1.6*10^-19;
E0 = 0.938*10^9;
Z0 = 120*Pi;
μ0 = 4*Pi*10^-7;
ϵ0 = 8.84*10^-12; 
μr = 1;
θH = 0;
ϵr = 1;

2) Beam Data

Et = 270*10^9;

γb = Et/E0

βb = Sqrt[1 - 1/γb^2]

WaveVelocityFactor = 1;

3) Impedance data

L = 1;   
b = 0.0015;    
d1 = 0.025;   
ρ1 = 10*10^-6;  
μr1 = μr;
θH1 = θH;
ϵr1 = ϵr;
ρ2 = Infinity;(* Resistivity of layer 2 *)
μr2 = μr;
θH2 = θH;
ϵr2 = ϵr;

σ1 = 1/ρ1 // N

σ2 = 1/ρ2 // N

4) General (multi-layer) impedance formula from BL implemented by AK

ParametersFor2layersBL = 
        Table[{{σ1, σ2}, {μr1, μr2}, {ϵr1, ϵr2}, {b, b + d1}}];

ZtcBL[f_, myParam_] := 
      Module[{shb, chb, shbpr, chbpr, kappa, kappaTilde, kappaTrans, cc, cpr, 
          s, spr, n, t, kappaTrans10, 
          oneminusG}, 
shb[k_, r_, a_] := 
          k a (BesselI[1, k r] BesselK[1,  k a] -  
                BesselI[1,  k a] BesselK[1, 
                    k r]); 
chb[k_, r_, a_] := 
          k a (((1/2) (BesselI[0, k a] + 
                          BesselI[2, k a])) BesselK[1,  k r] -  
                BesselI[1,  
                    k r] ((1/2) (-BesselK[0, k a] - 
                          BesselK[2, k a]))); 
shbpr[
            k_, r_, a_] := 
          k a ((1/2) (BesselI[0, k r] + BesselI[2, k r]) BesselK[
                    1, k a] - 
                (1/2) BesselI[1, 
                    k a] (-BesselK[0, k r] - 
                      BesselK[2, k r])); 
chbpr[k_, 
            r_, a_] := 
          k a (-(1/4) (BesselI[0, k r] + 
                      BesselI[2, k r]) (-BesselK[0, k a] - 
                      BesselK[2, k a]) + 
                (1/4) (BesselI[0, k a] + 
                      BesselI[2, k a]) (-BesselK[0, k r] - 
                      BesselK[2, k r])); 
kappa[m_] := 
          Sqrt[I myParam[[1, m]] myParam[[2, 
                  m]] μ0*2*Pi*f]; 
kappaTilde[m_] := 
          If[myParam[[1, m]] == 0, 
            1/(myParam[[4, m]] myParam[[2, m]]), 
            kappa[m]/myParam[[2, m]]]; 
kappaTrans[m1_, 
            m2_] := kappaTilde[m1]/kappaTilde[m2]; 
cc[m_] := 
          chb[kappa[m], myParam[[4, m + 1]], 
            myParam[[4, m]]]; 
cpr[m_] := 
          chbpr[kappa[m], myParam[[4, m + 1]], 
            myParam[[4, m]]]; 
s[m_] := 
          shb[kappa[m], myParam[[4, m + 1]], 
            myParam[[4, m]]]; 
spr[m_] := 
          shbpr[kappa[m], myParam[[4, m + 1]], 
            myParam[[4, m]]]; 
n = 
          Length[myParam[[1]]]; 
t = 
          Table[0, {i, 1, n}]; 
t[[n - 1]] = 
          Simplify[-(cpr[n - 1] + 
                    kappaTrans[n, n - 1] cc[n - 1])/(spr[n - 1] + 
                    kappaTrans[n, n - 1] s[
                        n - 1])]; 
For[m = n - 2, 
          m ≥  1, 
          t[[m]] = 
            Simplify[-(cpr[m] - 
                      t[[m + 1]] kappaTrans[m + 1, m] cc[m])/(spr[
                        m] - t[[m + 1]] kappaTrans[m + 1, m] s[
                          m])]; m--]; 
kappaTrans10 = 
          kappaTilde[1] myParam[[4, 1]]; 
oneminusG = 
          (2/(1 - kappaTrans10 t[[1]])); 
I L (Z0/(2 π myParam[[4, 1]]^2)) 
          oneminusG
]

ZtcBLImplementedByAK[f_] := ZtcBL[f, ParametersFor2layersBL]

5) General (1-layer) impedance formula from Zotter implemented by EM

ParametersFor1layerZotter = {L, 
          b, ρ1, μr, θH, ϵr};

Clear[f];
ZtcZotter1layer[f_, Param_]:=Module[{σ, μprime, μ, ϵprime, ϵ, 
            k, δ, ν, x1, x2, s, Q1, P1, Q2, P2, 
            y, η2, α2, Qη, Qα, AA, BB, CC},
σ = 1/Param[[3]] // N;
μprime = 
              Param[[4]]*(1 + I*Tan[Param[[5]]]);
μ = μ0*μprime;
ϵprime = 
              Param[[6]] + σ/(I*2*Pi*f*ϵ0);
ϵ = ϵ0*ϵprime;
k = 2*Pi*f/(βb*c);
δ = Sqrt[1/(Pi*f*μ*σ)];
ν = 
              k*Sqrt[1 - ϵprime*μprime*(βb^2)];
x1 = k*Param[[2]]/γb;
x2 = ν*Param[[2]];
a = 10^-5;
s = k*a/γb;
Q1 = (D[BesselK[1, t], t] /. t -> x1)/
                BesselK[1, x1];
P1 = (D[BesselI[1, t], t] /. t -> x1)/
                BesselI[1, x1];
Q2 = (D[BesselK[1, t], t] /. t -> x2)/
                BesselK[1, x2];
P2 = (D[BesselI[1, t], t] /. t -> x2)/
                BesselI[1, x2];
y = ν*d;
η2 = 0;    
  
  
α2 = 0;   
  
  
Qη = (Q2 - η2*P2)/(1 - η2);
Qα = (Q2 - α2*P2)/(1 - α2);
AA = (k*μprime*Qη)/(ν*γb*
                    P1);
BB = (a^2)*βb*(γb^2)*
                BesselI[1, 
                    x1]/((Param[[2]]^2)*(BesselI[1, s]^2)*
                      BesselK[1, 
                        x1]*βb*((x1*x2*γb*ν)^2)*
                      P1*(P1 - Q1));
CC = ((ν*γb*x2 - 
                        k*x1)^2) - ((βb*x1*
                          x2)^2)*(ν*γb*P1 - 
                      k*μprime*Qη)*(ν*γb*P1 - 
                      k*ϵprime*Qα);
(I*Param[[1]]*
                  Z0/(Pi*Param[[2]]^2))*βb*(1 - AA)/(BB*
                    CC)
]

ZtcZotter1layerDerivedAndImplementedByEM[f_] := 
      ZtcZotter1layer[f, ParametersFor1layerZotter]

6) General (2-layer) impedance formula from Zotter implemented by EM

ParametersFor2layersZotter = 
        Table[{{σ1, σ2}, {μr1, μr2}, {θH1, θH2}, {ϵr1, ϵr2}, {b, b + d1}}];

Clear[f];
ZtcZotter2layers[f_, Param_]:=Module[{β, γ, σ2, σ3, μp2, μp3, μ2, μ3, ϵp2, ϵp3, ϵ2, ϵ3, 
            k, ν2, ν3, x1, x2, x3, x4, a, s, Q1, P1, Q2, P2, Q4, P4, 
            K32, I32, Q32, P32, η3, α3, Eq1, Eq2, Eq3, Eq4, 
            sol, α2sol, η2sol, E2sol, G2sol, α1, η1, 
            ZtcTotal, ZtcSC},
β = βb*WaveVelocityFactor;
γ = 1/Sqrt[1 - β^2];
σ2 = Param[[1, 1]];
σ3 = Param[[1, 2]];
μr2 = Param[[2, 1]];
μr3 = Param[[2, 2]];
θH2 = Param[[3, 1]];
θH3 = Param[[3, 2]];
ϵr2 = Param[[4, 1]];
ϵr3 = Param[[4, 2]];
μp2 = μr2*(1 + 
                    I*Tan[θH2]);
μp3 = μr3*(1 + I*Tan[θH3]);
μ2 = μ0*μp2;
μ3 = μ0*μp3;
Clear[f];
ϵp2 = ϵr2 + σ2/(I*2*Pi*
                      f*ϵ0);
ϵp3 = ϵr3 + σ3/(I*2*Pi*
                      f*ϵ0);
ϵ2 = ϵ0*ϵp2;
ϵ3 = ϵ0*ϵp3;
k = 2*Pi*f/(β*c);
ν2 = 
              k*Sqrt[1 - ϵp2*μp2*(β^2)];
ν3 = 
              k*Sqrt[1 - ϵp3*μp3*(β^2)];
x1 = k*Param[[5, 1]]/γ;
x2 = ν2*Param[[5, 1]];
x3 = ν2*Param[[5, 2]];
x4 = ν3*Param[[5, 2]];
a = 10^-5;
s = k*a/γ;
Q1 = (D[BesselK[1, t], t] /. t -> x1)/
                BesselK[1, x1];
P1 = (D[BesselI[1, t], t] /. t -> x1)/
                BesselI[1, x1];
Q2 = (D[BesselK[1, t], t] /. t -> x2)/
                BesselK[1, x2];
P2 = (D[BesselI[1, t], t] /. t -> x2)/
                BesselI[1, x2];
Q4 = (D[BesselK[1, t], t] /. t -> x4)/
                BesselK[1, x4];
P4 = (D[BesselI[1, t], t] /. t -> x4)/
                BesselI[1, x4];
K32 = BesselK[1, x3]/BesselK[1, x2];
I32 = BesselI[1, x3]/BesselI[1, x2];
Q32 = (D[BesselK[1, t], t] /. t -> x3)/
                BesselK[1, x2];
P32 = (D[BesselI[1, t], t] /. t -> x3)/
                BesselI[1, x2];
η3 = 0;    
α3 = 0;   
Eq1 = γ*ν2*x2*
                  E2*(1 - α2) + γ*ν2*x1*x2*β*
                  G2*(1 - η2)*P1 - k*x1*E2*(1 - α2) - 
                k*x1*x2*β*μp2*G2*(Q2 - η2*P2);
Eq2 = γ*ν2*x2*
                  G2*(1 - η2) + γ*ν2*x1*
                  x2*β*(Q1 - P1 + P1*E2*(1 - α2)) - 
                k*x1*G2*(1 - η2) - 
                k*x1*x2*β*ϵp2*E2*(Q2 - α2*P2);
Eq3 = ν3*x4*
                  E2*(K32 - α2*I32) + ν3*x3*x4*β*μp2*
                  G2*(Q32 - η2*P32) - ν2*x3*
                  E2*(K32 - α2*I32) - ν2*x3*x4*β*μp3*
                  G2*(K32 - η2*
                        I32)*(Q4 - η3*P4)/(1 - η3);
Eq4 = ν3*x4*G2*(K32 - η2*I32) + ν3*x3*
                  x4*β*ϵp2*
                  E2*(Q32 - α2*P32) - ν2*x3*
                  G2*(K32 - η2*I32) - ν2*x3*
                  x4*β*ϵp3*
                  E2*(K32 - α2*
                        I32)*(Q4 - α3*P4)/(1 - α3);
sol = 
              Flatten[Solve[{Eq1 == 0, Eq2 == 0, Eq3 == 0, 
                    Eq4 == 0}, {α2, η2, E2, G2}]];
α2sol = α2 /. sol;
η2sol = η2 /. sol;
E2sol = E2 /. sol;
G2sol = G2 /. sol;
α1 = 1 - E2sol*(1 - α2sol);
η1 = G2sol*(1 - η2sol);
ZtcTotal = -((I*L*Z0*k*BesselI[1, s])/(π*
                          a*β*γ))*((BesselK[1, 
                            s] - α1*BesselI[1, s]*
                            BesselK[1, x1]/
                              BesselI[1, x1])*(1 - β*βb)/
                        s + (D[BesselI[1, t], t] /. 
                          t -> s)*η1*(BesselK[1, x1]/
                          BesselI[1, x1])*(β - βb));
ZtcSC = -((I*L*Z0*k*BesselI[1, s])/(π*
                          a*β*γ))*((BesselK[1, s] - 
                        BesselI[1, s]*
                          BesselK[1, x1]/
                            BesselI[1, x1])*(1 - β*βb)/
                      s);
ZtcTotal - ZtcSC
]

ZtcZotter2layersDerivedAndImplementedbyEM[f_] := 
      ZtcZotter2layers[f, ParametersFor2layersZotter]

7) Comparison between the formulae at a certain frequency

fexample = 10^6;

ZtcBLImplementedByAK[fexample]

ZtcZotter1layerDerivedAndImplementedByEM[fexample]

ZtcZotter2layersDerivedAndImplementedbyEM[fexample]