(* twocyl.nb *)

Needs["Graphics`ComplexMap`"]

Clear[V,g,aa,x1,x2,r1,r2,d,d1,
      d2,nec,nef,neq,tcb,tcf]
Off[General::spell1]

V = 1. (* Potential at upper circle (cylinder) in Volt;
 		  potential on lower circle (cylinder) = 0 ! *) ; 
 		  
d = 6. (* Distance between upper cylinder centre and
          lower cylinder centre in mm :
          for 0 < d *) ; 
          
r1 = 1/6 d (* Radius of upper cylinder in mm :
           for 0 < r1+r2 < d *) ;
           
r2 = 1/3 d (* Radius of lower cylinder in mm :
           for 0 < r2+r1 < d *) ;

nef = 8 (* Number of field lines between plates;
           variable needed in the pictures: 
           eflandeqpl14, efl14, absefateflu14, absefateflv14, 
           efuatefl14, efvatefl14 *);
  
neq = 10 (* Number of equipotential lines between plates;
           variable needed in the pictures: 
           eflandeqpl14, eqpll14, absefateqplu14, 
           absefataqplv14, efuateqpl14, efvateqpl14 *);  

nec = {20,20}  (* Number of lines in perspectivic view, 
				  eflandeqpl14 *);

tcf = .005	(* Relative thickness of field or equipotential 
			   lines in color diagrams *) ;
			   
tcb = .011	(* Relative thickness of lines representing 
		   	   boundaries, electrodes, pole shoes, etc.	*) ;

           (* Calculation of global variables *)

frx1x2 = FindRoot[{Cosh[x1]+Coth[x2] Sinh[x1] == d/r1,
               Cosh[x2]+Coth[x1] Sinh[x2] == d/r2},
              {x1,1},{x2,1}];
x1 = x1/.frx1x2 ; x2 = -x2/.frx1x2 ;
aa = r1 Sinh[x1]//N ;
d1 = aa Coth[x1]; d2 = aa Coth[x2];
 
                   (* Functions *) 

cm14[z_] = aa I Tanh[z/2]; (* Conformal map *)
ef14[z_] = 2 V/(I aa(x1-x2)) 1000 Conjugate[Cosh[z/2]^2];
                                (* Function of el. field *)
                              
         (* Variable needed in all pictures *)

g = 3.7 aa;

Print[ColumnForm[
{SequenceForm["Radius of upper cylinder r1 = ",r1 ," mm with"],
 SequenceForm["centre at d1 = ",d1 ," mm above median plane;"]}]]
 
Print[ColumnForm[
{SequenceForm["Radius of upper cylinder r2 = ",r2 ," mm with"],
 SequenceForm["centre at d2 = ",d2 ," mm below median plane;"]}]]

(* Section with black & white graphics *)

Print[ColumnForm[{"Absolute electric field times signum",
                  "of voltage on rim of upper cylinder"}]]
ppupabscyl14= ParametricPlot[
	  {(Arg[I (cm14[x1 + I y] - I d1)])/Degree,
	   Sign[V] Abs[ef14[x1 + I y]]},
	  {y,-Pi,Pi-.001}, 
	  PlotRange -> {{-180,180},{0,Automatic}},
      AxesLabel -> {"ϕ [deg]","E [V/m]"}]; 

Print["Electric field lines & equipotentials "]    
eflandeqpl14 = CartesianMap[cm14[#] &,
       {x2,x1,(x1-x2)/(neq-1)}, {-Pi,Pi,2 Pi/(2 nef-1)},
       AxesLabel -> {"u[mm]","v[mm]"},
       PlotRange -> {{-g,g}, {-g,g}}, 
       Epilog -> {{GrayLevel[.9],
                   Disk[{0,d1},r1], Disk[{0,d2},r2]},
                  {Thickness[tcb], 
                   Circle[{0,d1},r1], Circle[{0,d2},r2]}}];   

Print[ColumnForm[{"Absolute electric field times signum",
                  "of voltage on rim of lower cylinder"}]]
ppdownabscyl14 = ParametricPlot[
	  {(Arg[-I (cm14[x2 + I y] - I d2)])/Degree,
	   Sign[V] Abs[ef14[x2 + I y]]},
	  {y,-Pi,Pi-.001}, 
	  PlotRange -> {{-180,180},{0,Automatic}},
      AxesLabel -> {"ϕ [deg]","E [V/m]"}];

Print[ColumnForm[{"Absolute electric field times signum of",
                  "voltage on rim of upper cylinder without units"}]]
Show[ppupabscyl14,
     PlotRegion -> {{.01,.99},{.01,.99}}, 
     AxesLabel -> {ϕ, E}, 
     Ticks -> None];

Print["Electric fieldlines & equipotentials without units"]
Show[eflandeqpl14, AxesLabel -> None, Ticks -> None];

Print[ColumnForm[{"Absolute electric field times signum of",
                  "voltage on rim of lower cylinder without units"}]]
Show[ppdownabscyl14, 
     PlotRegion -> {{.01,.99},{.01,.99}},
     AxesLabel -> {ϕ, E}, 
     Ticks -> None];

Print[ColumnForm[{"Absolute value of electric field",
                  "in 3D seen from right median plate"}]]
absef3D14 = ParametricPlot3D[Evaluate[
 		{Re[cm14[x + I y]],Im[cm14[x + I y]],
 		 Abs[ef14[x + I y]]}],
  		{x,x2,x1},{y,-Pi+.0001,Pi-.0001},
        PlotRange -> {{-g,g},{-g,g},All},
  		AxesLabel -> {"u[mm]","v[mm]","|E| [V/m]"}, 
  		BoxRatios -> {1,1,1}, 
  		PlotPoints -> nec,
  		ViewPoint -> {-2, 1, 2.3}];

hli =  	Table[Hue[N[1 - (k-1)/(neq - 2)]],{k,neq-2}] ;
hlc = 	Join[{GrayLevel[0]}, hli, {GrayLevel[0]}] ;
hlt = 	Table[Thickness[tcf],{k,neq}] ;
styleq = Transpose[{hlc,hlt}] ;

Print[neq," Equipotential lines (including plate and circle)"]
eqpl14 = ParametricPlot[Evaluate[Table[
		{Re[cm14[x2+k (x1-x2)/(neq-1) + I y]], 
		 Im[cm14[x2+k (x1-x2)/(neq-1) + I y]]},
		{k,0,neq-1}]],{y,-Pi,Pi}, 
		AspectRatio -> Automatic,
		AxesLabel -> {"u [mm]","v [mm]"},
		PlotRange -> {{-g,g}, {-g,g}},
		PlotStyle -> styleq,   
		Epilog -> 
		{{GrayLevel[.9],
         Disk[{0,d1},r1], Disk[{0,d2},r2]},
        {Thickness[tcb], 
         Circle[{0,d1},r1], Circle[{0,d2},r2]}}];

Print[ColumnForm[
	  {"Modified absolute electric fields on equipotential lines as", 
	  "functions of u (colours match the colours above)"}]]
absefateqplu14 = ParametricPlot[Evaluate[Table[
		{Re[cm14[x2+k (x1-x2)/(neq-1) + I y]], 
	     Sign[Im[cm14[x2+k (x1-x2)/(neq-1) + I y]]]*
		 Abs[ef14[x2+k (x1-x2)/(neq-1) + I y]]},
		{k,0,neq-1}]],{y,-Pi,Pi},
		AxesLabel -> {"u [mm]","Sign(v).|E(u)| [V/m]"},
		PlotRange -> {{-g,g},All},
		PlotStyle -> styleq];

Print[ColumnForm[
	  {"Absolute electric fields on equipotential lines as", 
	  "functions of v (colours match the colours above)"}]]
absefateqplv14 = ParametricPlot[Evaluate[Table[
		{Im[cm14[x2+k (x1-x2)/(neq-1) + I y]], 
		 Abs[ef14[x2+k (x1-x2)/(neq-1) + I y]]},
		{k,0,neq-1}]],{y,-Pi,Pi},
		AxesLabel -> {"v [mm]","|E(v)| [V/m]"},
		PlotRange -> {{-g,g},All},
		PlotStyle -> styleq];

Print[ColumnForm[
	  {"Electric fields (v-component) on equipotential lines as", 
	  "functions of u (colours match the colours above)"}]]
efvateqpl = ParametricPlot[Evaluate[Table[
		{Re[cm14[x2+k (x1-x2)/(neq-1) + I y]],
		 Im[ef14[x2+k (x1-x2)/(neq-1) + I y]]},
		{k,0,neq-1}]],{y,-Pi,Pi},
		AxesLabel -> {"u[mm]","\!\(E\_v\)V/m]"},
		PlotRange -> {{-g,g},All},
		PlotStyle -> styleq];

Print[ColumnForm[
	  {"Electric fields (u-component) on equipotential lines as",
	   "functions of v (colours match the colours above)"}]]
efuateqpl = ParametricPlot[Evaluate[Table[
		{Im[cm14[x2+k (x1-x2)/(neq-1) + I y]],
		 Re[ef14[x2+k (x1-x2)/(neq-1) + I y]]},
		{k,0,neq-1}]],{y,-Pi,Pi},
		AxesLabel -> {"v[mm]","\!\(E\_u\)V/m]"},
		PlotRange -> {{-g,g},All},
		PlotStyle -> styleq];

hli =	Table[Hue[N[1 - (k-1)/(nef-2)]],{k,nef-2}];
hlc = 	Join[{GrayLevel[0]}, hli, {GrayLevel[0]},
             Reverse[hli],{GrayLevel[0]}] ;
hlt = 	Table[Thickness[tcf],{k,2nef-1}] ;
stylfe = Transpose[{hlc,hlt}] ;

Print[2nef-1," Electric field lines"]
efl14 = ParametricPlot[Evaluate[Table[
		{Re[cm14[(k-1)/(nef-1) Pi I + x]], 
		  Im[cm14[(k-1)/(nef-1) Pi I + x]]},
		{k,-nef+2,nef}]],{x,x2,x1}, 
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [mm]","v [mm]"},
		PlotStyle -> stylfe,
		Epilog -> 
		{{GrayLevel[.9],
         Disk[{0,d1},r1], Disk[{0,d2},r2]},
        {Thickness[tcb], 
         Circle[{0,d1},r1], Circle[{0,d2},r2]}}];

Print[ColumnForm[
	  {"Absolute electric fields on field lines as",
	  "functions of u (colours match the colours above)"}]]
absefateflu14 = ParametricPlot[Evaluate[Table[
		{Re[cm14[(k-1)/(nef-1) Pi I + x]], 
		 Abs[ef14[(k-1)/(nef-1) Pi I + x]]},
		{k,-nef+2,nef}]],{x,x2,x1},
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [mm]","|E(u)| [V/m]"},(* units ?????? *)
		PlotStyle -> stylfe];

Print[ColumnForm[
	  {"Absolute electric fields on field lines as", 
	  "functions of v (colours match the colours above)"}]]
absefateflv14 = ParametricPlot[Evaluate[Table[
		{ Im[cm14[(k-1)/(nef-1) Pi I + x]], 
		 Abs[ef14[(k-1)/(nef-1) Pi I + x]]},
		{k,nef,1,-1}]],{x,x2,x1},
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"v [mm]","|E(v)| [V/m]"},
		PlotStyle -> stylfe];

Print[ColumnForm[
	  {"Electric fields (v-component) on field lines as", 
	  "functions of u (colours match the colours above)"}]]
efvatefl14 = ParametricPlot[Evaluate[Table[
		{Re[cm14[(k-1)/(nef-1) Pi I + x]], 
		 Im[ef14[(k-1)/(nef-1) Pi I + x]]},
		{k,-nef+2,nef}]],{x,x2,x1},
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [mm]","\!\(E\_v\) [V/m]"},
		PlotStyle -> stylfe];

Print[ColumnForm[
	  {"Electric fields (u-component) on field lines as", 
	  "functions of v (colours match the colours above)"}]]
efuatefl14 = ParametricPlot[Evaluate[Table[
		{Im[cm14[(k-1)/(nef-1) Pi I + x]], 
		 Re[ef14[(k-1)/(nef-1) Pi I + x]]},
		{k,-nef+2,nef}]],{x,x2,x1},
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"v [mm]","\!\(E\_u\) [V/m]"},
		PlotStyle -> stylfe];

Clear[up,vp,frefp,efp,s,rim1,rim2]

up = 0; vp = 1.75;
(* above you may insert values of the Cartesian coordinates 
u and v (in mm) of a point of your choice within the area 
shown in Fig. eflandeqpl14 *)
 
rim1 = up^2 + (vp - d1)^2 ; rim2 = up^2 + (vp - d2)^2 ;

Which[rim1 < r1^2,  
	  Print[rim1 "= up^2 + (vp - d1)^2 < r1^2 is unphysical.",⚠],
      rim2 < r2^2,
	  Print[rim2 "= up^2 + (vp - d2)^2 < r2^2 is unphysical.",⚠],  
	  True, 
	  
efp = 2 ArcTanh[(up + I vp)/(aa I)];
valef = ef14[efp]//N ;
Print[
SequenceForm["At point (up,vp) = (",up ,", ",vp ,") mm : ",
"Potential = ",N[V(Re[efp]-x2)/(x1-x2)], " V;"]];
Print[SequenceForm["Field (\!\(E\_u\),\!\(E\_v\)) = (", 
Re[valef],", ", Im[valef],") V/m "]];
	  
caleqpl14 = ParametricPlot[Evaluate[
	  {Re[cm14[Re[efp] + I y]], Im[cm14[Re[efp] + + I y]]}], 
	  {y,-Pi,Pi}, 
	  PlotRange -> All,
	  PlotStyle -> Dashing[{.01}], 
	  DisplayFunction -> Identity];
calefl14 = ParametricPlot[Evaluate[
	  {Re[cm14[I Im[efp] + s]], Im[cm14[I Im[efp] + s]]}],
	  {s,x1,x2}, DisplayFunction -> Identity ];
Show[caleqpl14, calefl14, 
     DisplayFunction -> $DisplayFunction,
	 AxesLabel -> {"u [mm]","v [mm]"},
	 PlotRange -> {{-g,g}, {-g,g}}, 
	 Epilog -> {{GrayLevel[.9],
                 Disk[{0,d1},r1], Disk[{0,d2},r2]},
                {Thickness[tcb], 
                 Circle[{0,d1},r1], Circle[{0,d2},r2]}}];   ]

On[General::spell1]