(* This notebook is used to self consistently derive stationary equilibria  within the framework of Ideal Magnetohydrodynamics in cylindrical geometry .  
              The input quantities needed  as functions of the radial distance, rho, are the q - profile, the axial magnetic field, the axial and azimuthal  velocity components, and the mass density .  
              The calculated quantities are the azimuthal magnetic field, the pressure, the axial current density, the total axial current, the radial electric field, its radial variation (shear),  and the shear of the ExB velocity .  
In particular, reversed - magnetic - 
            shear equilibrium quantities  are evaluated here;  
          accordingly, the prescribed  q-profile  exhibits a maximum on the magnetic axis, rho = 0, and a minimum inside the plasma column .  Also, a peaked profile for the axial velocity and a Gaussian - 
            like profile for the azimuthal velocity are chosen. *)

Developer`ClearCache["Symbolic"]

q[rho_] := 
        qc (1 - a*rho^2  + b*rho^3) ;
Bz[rho_] := 
        Bz0 (1 + c (1 - rho^2))^(1/
                2);
uz[rho_] := uz0 (1 - rho^3)^3;
utheta[rho_] := 
        utheta0*rho (1 - rho) Exp[-((rho - rhomn)^2)/
              h];
m[rho_] := m0*((1 - rho^3)^3)

(* Btheta  = (Sqrt[e ])*rho*Bz[rho]/q[rho] *)

bth = 
        FullSimplify[(Sqrt[e ])*rho*
            Bz[rho]/q[rho]];
Btheta[rho_] := bth

Plot[
      bth /.  {a -> 8, b -> 10.67, c -> -0.0975, 
          e -> 0.11, Bz0 -> 1, qc -> 6}, {rho, 0, 
        1}]

(* 
        Calculation of the flow term in the pressure expression *)

Mthet = 
        FullSimplify[
          m[rho]*(utheta[rho]^2)/(Btheta[
                    rho]^2)];
FullSimplify[((1/mi0) - Mthet) (1/
              r0) ((Btheta[rho]^2)/
              rho)];
Int1 = Integrate[%, rho];
Int2 = Int1 /.  rho -> 0;

P = (-1/(2*
                    mi0)) (Bz0 (1 + 
                      c (1 - 
                            rho^2))  + (((Bz0 (Sqrt[e]) rho (Sqrt[1 + c - c (rho^2)]))/(qc - a qc (rho^2) + b qc (rho^3))))^2) + ((((Bz0^2) (1 + c))/(2 mi0))) - (Int1 - 
              Int2);
P1 = P /.  rho -> 1;
P0 = (P - P1) /.  rho -> 0;
Pn = (P - P1)/P0;
Plot[
      Pn /.  {a -> 8, b -> 10.67, c -> -0.0975, 
          m0 -> 8.35*10^(-8), e -> 0.11, h -> 0.001, 
          mi0 -> 4*Pi*10^(-7), utheta0 -> 4*^4, 
          r0 -> 1, Bz0 -> 1, qc -> 6, 
          rhomn -> 0.5}, {rho, 0, 1}]

(* Jz = (1/(mi0*R0*
                    rho)) (∂_rho) (Bz*(rho^2)/q) *)

j = 
        FullSimplify[(1/(mi0*R0*rho)) (∂_rho) (Bz[
                    rho]*(rho^2)/q[rho])];
j0 = j /.  rho -> 0;
jz = FullSimplify[j/j0]

Plot[
      jz /.  {a -> 8, b -> 10.67, c -> -0.0975}, {rho, 0, 
        1}]

(* I = 
          Integrate[
            Pi*(r0^2)*rho*Jzn[a, b, c, mi0, R0, Bz0, rho], {rho, 0, 1}] *)

Iz = 
      NIntegrate[
        Pi*(1^2)*rho*j /.  {a -> 8, b -> 10.67, 
            c -> -0.0975, Bz0 -> 1, qc -> 6, R0 -> 3, 
            mi0 -> 4*Pi*10^(-7)}, {rho, 0, 1}]

(*  Er = Btheta*uz - Bz*utheta *)

Er = FullSimplify[Btheta[rho]*uz[rho] - Bz[rho]*utheta[rho]]

Plot[
      Er /.  {a -> 8, b -> 10.67, c -> -0.0975, 
          e -> 0.11, h -> 0.001, Bz0 -> 1, qc -> 6, 
          utheta0 -> 4*^4, uz0 -> 1*^5, rhomn -> 0.5}, {rho, 
        0, 1}, PlotRange -> {{0, 1}, {6*^3, -6*^3}}]

SEr = (∂_rho) Er

Plot[
      SEr /.  {a -> 8, b -> 10.67, c -> -0.0975, 
          e -> 0.11, h -> 0.001, Bz0 -> 1, qc -> 6, 
          utheta0 -> 4*^4, uz0 -> 1*^5, rhomn -> 0.5}, {rho, 
        0, 1}, PlotRange -> {{0, 1}, {3*^5, -3*^5}}]

(* Shear of the E×B *)

k1 = (∂_rho) ((Er*Bz[rho])/(Bz[rho]^2  + 
                  Btheta[rho]^2));
k2 = (∂_rho) ((Er*
                  Btheta[rho])/((Bz[rho]^2  + 
                    Btheta[rho]^2)));
omega = (Sqrt[(k1^2  + k2^2)]);
Plot[
      omega /.  {a -> 8, b -> 10.67, c -> -0.00975, 
          e -> 0.11, h -> 0.001, Bz0 -> 1, qc -> 6, 
          utheta0 -> 4*^4, uz0 -> 1*^5, rhomn -> 0.5}, {rho, 
        0, 1}]

FindMinimum[-omega /.  {a -> 8, b -> 10.67, 
          c -> -0.0975, e -> 0.11, h -> 0.001, 
          Bz0 -> 1, qc -> 6, utheta0 -> 4*^4, 
          uz0 -> 1*^5, rhomn -> 0.5}, {rho, 0.45}]