(* plandcyl.nb *)

Needs["Graphics`ComplexMap`"]

Clear[V,aa,x0,r0,d,nec,nef,neq,tcb,tcf,g]
Off[General::spell1]

                      (* Input data *)

V = 1 (* Potential at circle (cylinder) in Volt,
 		 pontential on lower plate (= median plane) = 0 *) ; 
 		 
d = 5 (* Distance of cylinder centre and plane in mm :
         for 0 < d *) ; 
         
r0 = 3/5 d (* Radius of cylinder in mm :
          for 0 < r0 < d *) ;

nec = {10,20} (* Number of lines used in 3D-image of 
            el. field`s absolute value: absef3D13 *);
            
nef = 11 (* Number of field lines between plates;
            variable needed in the pictures: 
            eflandeqpl13, efl13, absefateflu13, 
            absefateflv13, efuatefl13, efvatefl13 *);
            
neq = 11 (* Number of equipotential lines between plates;
            variable needed in the pictures: 
            eflandeqpl13, eqpll13, absefateqplu13, 
            absefataqplv13, efuateqpl13, efvateqpl13 *);
            	
tcf = .005 (* Relative thickness of field or equipotential 
		      lines in color diagrams *);
		      
tcb = .011 (* Relative thickness of lines representing 
		   	  boundaries, electrodes, pole shoes, etc. *);

                 (* Global variables *)

x0 = ArcCosh[d/r0] //N ;
aa = r0 Sinh[x0] //N ;

                     (* Functions *) 
                     
cm13[z_] = aa I Tanh[z/2]; (* Conformal map *)
ef13[z_] = -V/x0 1000 Conjugate[2 Cosh[z/2]^2/(I aa)];
                           (* Function of el. field *)

            (* Variable needed in all pictures *)
            
g = 3.7 aa ;

Print[SequenceForm[
"Distance of cylinder center and median plane d = ", d, " mm;"]]
Print["Radius of cylinder = r0 = ",r0, " mm."]

(* Section with black & white graphics *)

Print[ColumnForm[{"Absolute electric field times signum",
                  "of potential on rim of cylinder"}]]
ppabscyl13 = ParametricPlot[
	  {(Arg[I (cm13[x0 + I y] - I d)])/Degree,
	   Sign[V] Abs[ef13[x0 + I y]]},
       {y,-Pi,Pi-.001}, 
       PlotRange -> {{-180,180},{0,Sign[V]Abs[ef13[x0]]}},
       AxesLabel -> {"ϕ [deg]","sgn(V).|E| [V/m]"}];

Print["Electric field lines & equipotentials "]     
eflandeqpl13 = CartesianMap[cm13[#] &, 
       {0,x0,x0/(neq-1)},{-Pi,Pi,2 Pi/(2 nef-1)},
       AxesLabel -> {"u [mm]","v [mm]"}, 
       PlotRange -> {{-g,g},{0,1.1 g}},
       Epilog -> {{GrayLevel[.9],
        		   Disk[{0,d},r0]},
        		  {Thickness[tcb], 
                   Line[{{-g,0},{g,0}}], 
                   Circle[{0,d},r0]}}];
        		 
Print["Electric field on median plate "]
ppdown13 = ParametricPlot[
	  {Re[cm13[I y]], -Im[ef13[I y]]}, 
      {y,-Pi+.001,Pi-.001}, 
       PlotRange -> {{-g,g},All}, 
       AxesLabel -> {"u [mm]","E [V/m]"}];

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

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

Print["Electric field of median plate without units"] 
Show[ppdown13, 
     PlotRegion -> {{.01,.99},{.01,.99}}, 
     Ticks -> None, 
     AxesLabel -> {"u","E"}];

(* Section with colour graphics *)

Print[ColumnForm[{"Absolute value of electric field",
                  "in 3D seen from right median plate"}]]
absef3D13 = ParametricPlot3D[Evaluate[
 		{Re[cm13[x + I y]],Im[cm13[x + I y]],
 		Abs[ef13[x + I y]]}],
  		{x,0,x0},{y,-Pi+.0001,Pi-.0001},
  		PlotRange -> {{-g,g},{-.1, g 1.1},All},
  		AxesLabel -> {"u[mm]","v[mm]","|E| [V/m]"}, 
  		BoxRatios -> {1,1,1}, 
  		PlotPoints -> nec,
  		ViewPoint -> {-1,.7,.5}];

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)"]
eqpl13 = ParametricPlot[ Evaluate[Table[
		{Re[cm13[k x0/(neq-1) + I y]], 
		 Im[cm13[k x0/(neq-1) + I y]]},
		{k,0,neq-1}]],{y,-Pi,Pi}, AspectRatio -> Automatic,
		AxesLabel -> {"u [mm]","v [mm]"},
		PlotRange -> {{-g,g}, {-.1, g 1.1}},
		PlotStyle -> styleq,
		Epilog -> {{GrayLevel[.9],Disk[{0,d},r0]},
		           {Thickness[tcb], 
                    Line[{{-g,0},{g,0}}], Circle[{0,d},r0]}}];

Print[ColumnForm[
	  {"Absolute electric fields on equipotential lines as",
      "functions of u (colours match the colours above)"}]]
absefateqplu13 = ParametricPlot[ Evaluate[Table[
		{Re[cm13[k x0/(neq-1) + I y]], 
		Abs[ef13[k x0/(neq-1) + I y]]},
		{k,0,neq-1}]],{y,-Pi,Pi},
		AxesLabel -> {"u [mm]","|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)"}]]
absefateqplv13 = ParametricPlot[ Evaluate[Table[
		{Im[cm13[k x0/(neq-1) + I y]], 
		Abs[ef13[k x0/(neq-1) + I y]]},
		{k,0,neq-1}]],{y,-Pi,Pi},
		AxesLabel -> {"v [mm]","|E(v)| [V/m]"},
		PlotRange -> {{-.1,g 1.1},All},
		PlotStyle -> styleq];

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

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

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"]
efl13 = ParametricPlot[ Evaluate[Table[
		{Re[cm13[s - .9999 Pi (k-1)/(nef-1) I ]], 
		 Im[cm13[s - .9999 Pi (k-1)/(nef-1) I ]]},
		{k,-nef+2,nef}]],{s,0,x0}, 
		PlotRange -> {0, 1.1 g}, 
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [mm]","v [mm]"},
		PlotStyle -> stylfe, 
		Epilog -> {{GrayLevel[.9],
        		    Disk[{0,d},r0]},
        		   {Thickness[tcb], 
                    Line[{{-g,0},{g,0}}], Circle[{0,d},r0]}}];

Print[ColumnForm[
	  {"Absolute electric fields on field lines as",
	  "functions of u (colours match the colours above)"}]]
absefateflu13 = ParametricPlot[ Evaluate[Table[
		{Re[cm13[s - .9999 Pi (k-1)/(nef-1) I ]], 
		Abs[ef13[s - .9999 Pi (k-1)/(nef-1) I ]]},
		{k,-nef+2,nef}]],{s,0,x0}, 
		PlotRange -> Automatic, 
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [mm]","|E(u)| [V/m]"},
		PlotStyle -> stylfe];

Print[ColumnForm[
	  {"Absolute electric fields on field lines as", 
	  "functions of v (colours match the colours above)"}]]
absefateflv13 = ParametricPlot[ Evaluate[Table[
		{Im[cm13[s - .9999 Pi (k-1)/(nef-1) I ]], 
		Abs[ef13[s - .9999 Pi (k-1)/(nef-1) I ]]},
		{k,-nef+2,nef}]],{s,0,x0}, 
		PlotRange -> Automatic, 
		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)"}]]
efvatefl13 = ParametricPlot[ Evaluate[Table[
		{Re[cm13[s - .9999 Pi (k-1)/(nef-1) I ]], 
		 Im[ef13[s - .9999 Pi (k-1)/(nef-1) I ]]},
		{k,-nef+2,nef}]],{s,0,x0}, 
		PlotRange -> Automatic, 
		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)"}]]
efuatefl13 = ParametricPlot[ Evaluate[Table[
		{Im[cm13[s - .9999 Pi (k-1)/(nef-1) I ]], 
		 Re[ef13[s - .9999 Pi (k-1)/(nef-1) I ]]},
		{k,-nef+2,nef}]],{s,0,x0},
		PlotRange -> {{-.01,1.1 g},Automatic},
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"v [mm]","\!\(E\_u\) [V/m]"},
		PlotStyle -> stylfe];

Clear[up,vp,frefp,efp,error,imefp,s,rim]

up = 5.1; vp = 7.5;
 (* 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. eflandeqpl13 *) 
				
rim = up^2 + (vp - d)^2 ;

Which[vp < 0,  
Print["vp = ",vp," < 0 is unphysical.",⚠],
      rim < r0^2, 
Print[SequenceForm["up^2 + (vp - d)^2 = ",rim,
                     " < r0^2  is unphysical.",⚠]],
      True, 

efp = 2 ArcTanh[(up + I vp)/(aa I)];
valef = ef13[efp]//N ;
Print[
SequenceForm["At point (up,vp) = (",up ,", ",vp ,") mm : ",
"Potential = ",N[V Re[efp]/x0], " V;"]];
Print[SequenceForm["Field (\!\(E\_u\),\!\(E\_v\)) = (", 
Re[valef],", ", Im[valef],") V/m "]];

caleqpl13 = ParametricPlot[Evaluate[
	  {Re[cm13[Re[efp] + I s]], Im[cm13[Re[efp] + I s]]}], 
	  {s,-Pi,Pi}, 
	  PlotStyle -> Dashing[{.01}], 
	  DisplayFunction -> Identity];
	  
calefl13 = ParametricPlot[Evaluate[
	  {Re[cm13[I Im[efp] + s]], Im[cm13[I Im[efp] + s]]}],
	  {s,0,x0}, DisplayFunction -> Identity]; 
	  
Show[caleqpl13,calefl13, 
     DisplayFunction -> $DisplayFunction,
	 AxesLabel -> {"u [mm]", "v [mm]"},
	 PlotRange -> {{-g,g}, {-.1, g .9}},
     Epilog -> {{GrayLevel[.9], Disk[{0,d},r0]},
        		{Thickness[tcb], 
                 Line[{{-g,0},{g,0}}], Circle[{0,d},r0]}}]; ]
                 
On[General::spell1]