(* \(README\ \(FIRST : \ 
                This\ code\ runs\ on\ Mathematica\ 5.0 . \ 
                    The\ input\ parameters\ are\ grouped\ into\ two\ \(sections : \ "MAIN INPUT LINE"\ and\ "INPUT LINE" . \ 
                        To\ change\ the\ values\ of\ these\ input\ parameters\ one\ may\ edit\ either\ the\ main\ input\ line\ or\ the\ input\ line . \ 
                        The\ main\ input\ line\ contains\ the\ values\ for\ μrel\ \((relative\ permeability)\)\)\)\), m (magnetic diple moment), σ (conductivity), R1 (interior radius), and R2(exterior radius) .  The five input line entries contain the different values of a (magnet radius) and L (magnet length) for the five cases (a) through (d) shown below .  All quatities are in SI units .  The execution of this program does take some time .  On a Pentium4 2.40 GHz computer with 1024MB of RAM this entire program took  minutes to complete .  The individual execution times are listed at the beginning of each section .  To return to the main text of " Electrodynamics of a Magent Moving through a Metallic Pipe," click here *)

Clear[k, κ, F, v, b0, b1, b2, b3, 
            b4, μrel, μ0, m, σ, R1, R2, a, L];

    μrel = 1.0; m = 1.00; σ = 5.0*10^7; R1 = 0.010; 
          R2 = 0.011;

μ0 = 
          4*π*10^-7; μ = μrel*μ0; α  = (Sqrt[1 - I*(v*μ*σ/k)]); β = α/μrel;
QI = 
            Im[((β^2*BesselK[1, R1 (Sqrt[(k^2)])]*
                        BesselK[1, 
                          R2 (Sqrt[(k^2)])] (BesselI[0, R1*α*(Sqrt[(k^2)])]*
                              BesselK[0, R2 *α*(Sqrt[(k^2)])] - 
                            BesselI[0, R2*α*(Sqrt[(k^2)])]*
                              BesselK[0, R1*α*(Sqrt[(k^2)])]) + 
                      BesselK[0, R1 (Sqrt[(k^2)])]*
                        BesselK[0, 
                          R2 (Sqrt[(k^2)])] (BesselI[1, R2*α*(Sqrt[(k^2)])]*
                              BesselK[1, R1 *α*(Sqrt[(k^2)])] - 
                            BesselI[1, R1*α*(Sqrt[(k^2)])]*
                              BesselK[1, R2*α*(Sqrt[(k^2)])]) + β*
                        BesselK[0, R1 (Sqrt[(k^2)])]*
                        BesselK[1, 
                          R2 (Sqrt[(k^2)])] (BesselI[1, R1*α*(Sqrt[(k^2)])]*
                              BesselK[0, R2 *α*(Sqrt[(k^2)])] + 
                            BesselI[0, R2*α*(Sqrt[(k^2)])]*
                              BesselK[1, R1*α*(Sqrt[(k^2)])]) - β*
                        BesselK[0, R2 (Sqrt[(k^2)])]*
                        BesselK[1, 
                          R1 (Sqrt[(k^2)])] (BesselI[1, R2*α*(Sqrt[(k^2)])]*
                              BesselK[0, R1 *α*(Sqrt[(k^2)])] + 
                            BesselI[0, R1*α*(Sqrt[(k^2)])]*
                              BesselK[1, 
                                R2*α*(Sqrt[(k^2)])]))/(β^2*
                        BesselI[1, R1 (Sqrt[(k^2)])]*
                        BesselK[1, 
                          R2 (Sqrt[(k^2)])] (BesselI[0, R2*α*(Sqrt[(k^2)])]*
                              BesselK[0, R1*α*(Sqrt[(k^2)])] - 
                            BesselI[0, R1*α*(Sqrt[(k^2)])]*
                              BesselK[0, R2*α*(Sqrt[(k^2)])]) + 
                      BesselI[0, R1 (Sqrt[(k^2)])]*
                        BesselK[0, 
                          R2 (Sqrt[(k^2)])] (BesselI[1, R2*α*(Sqrt[(k^2)])]*
                              BesselK[1, R1 *α*(Sqrt[(k^2)])] - 
                            BesselI[1, R1*α*(Sqrt[(k^2)])]*
                              BesselK[1, R2*α*(Sqrt[(k^2)])]) + β*
                        BesselI[0, R1 (Sqrt[(k^2)])]*
                        BesselK[1, 
                          R2 (Sqrt[(k^2)])] (BesselI[1, R1*α*(Sqrt[(k^2)])]*
                              BesselK[0, R2 *α*(Sqrt[(k^2)])] + 
                            BesselI[0, R2*α*(Sqrt[(k^2)])]*
                              BesselK[1, R1*α*(Sqrt[(k^2)])]) + β*
                        BesselI[1, R1 (Sqrt[(k^2)])]*
                        BesselK[0, 
                          R2 (Sqrt[(k^2)])] (BesselI[1, R2*α*(Sqrt[(k^2)])]*
                              BesselK[0, R1 *α*(Sqrt[(k^2)])] + 
                            BesselI[0, R1*α*(Sqrt[(k^2)])]*
                              BesselK[1, R2*α*(Sqrt[(k^2)])])))];
F[v_] := (μ0*m^2/(4*Pi^2))*
              NIntegrate[
                k^3*((Sin[k*L/2]/(k*L/2))^2)*((BesselI[
                            1, (Sqrt[k])^2*a]/((Sqrt[k])^2*a/2))^2)*2*QI, {k, 0, 
                  1000}];

Clear[a, L];

    a = 0.000001; L = 0.00000001;

P1 = 
          Plot[F[v],  {v, 0, 100}, 
            PlotStyle -> {RGBColor[0, 1, 0]}]

Clear[a, L];

   a = 0.65*R1;  L = 0.0000001;

P2 = 
          Plot[F[v],  {v, 0, 100}, 
            PlotStyle -> {RGBColor[1, 0, 0]}]

Clear[a, L];

   a = 0.96*R1;  L = 1.2*R1;

P3 = 
          Plot[F[v],  {v, 0, 100}, 
            PlotStyle -> {RGBColor[1, 0.5, 0]}]

Clear[a, L];

    a =  0.6*R1;  L = 2*a;

P4 = 
          Plot[F[v],  {v, 0, 100}, 
            PlotStyle -> {RGBColor[0, 0, 1]}]

Clear[a, L];

   a =  0.6*R1;  L = 4*a;

P5 = 
          Plot[F[v],  {v, 0, 100}, 
            PlotStyle -> {RGBColor[ .8, 0, 1]}]

G[v_] := 
      45*(R2 - R1)*m^2*
        v*μ0^2*σ/(1024*R1^4)
P6 = 
      Plot[G[v],  {v, 0, 48}, 
        PlotStyle -> {{GrayLevel[0.4], 
              Dashing[{0.05, 0.05}]}}]

Show[P1, P2, P3, P4, P5, P6, 
      Graphics[Text[
          StyleForm["a", FontSize -> 10, 
            FontFamily -> "Times"], {90, 2.7}]], 
      Graphics[Text[
          StyleForm["b", FontSize -> 10, 
            FontFamily -> "Times"], {90, 4.85}]], 
      Graphics[Text[
          StyleForm["c", FontSize -> 10, 
            FontFamily -> "Times"], {90, 5.5}]], 
      Graphics[Text[
          StyleForm["d", FontSize -> 10, 
            FontFamily -> "Times"], {90, 7.15}]], 
      Graphics[Text[
          StyleForm["e", FontSize -> 10, 
            FontFamily -> "Times"], {90, 8.5}]], Frame -> True, 
      FrameLabel -> {StyleForm["\nVelocity (m/s)", 
            FontSize -> 10, FontFamily -> "Times"], 
          StyleForm["Drag Force (N)", FontSize -> 10, 
            FontFamily -> "Times"]}, 
      FrameTicks -> {{0, 20, 40, 60, 80, 100}, {0, 2, 4, 6, 8, 10, 12, 
            14, 16}, {}, {}}, AxesOrigin -> {0, 0}, 
      PlotRange -> {{0, 100}, {0, 10.6}}, 
      AspectRatio -> 1]