(* sharpben.nb *)

Needs["Graphics`ComplexMap`"]
Off[General::spell1]

Clear[V,a,b,c,d,e,q,t,frc,frd,sgn8,mult,tex1]

		        (* Input data *)
		        
V = 1 (* Potential on upper plate in Volt, 
		 potential on lower plate = 0 *);
a = 3 (* Distance between vertical plate parts in mm :
         for 10^-5 < a < 10^5 *); 
b = a/2 (* Distance between horizontal plate parts in mm :
           for a/10 < b < a/1.0001 *);

nef = 11 (* Number of field lines between plates;
            variable needed in the pictures: 
            eflandeqpl8, efl8, absefateflu8,
            absefateflv8, efuatefl8, efvatefl8 *);
            
neq = 9 (* Number of equipotential lines between plates;
           variable needed in the pictures: 
           eflandeqpl8, eqpll8, absefateqplu8, 
           absefataqplv8, efuateqpl8, efvateqpl8 *);
            
nec = {30,15} (* Number of lines used in 3D-image of 
                 el. field`s absolute value: absef3D8 *);

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 *)
                  
mult = V 1000/a; tex1 = "E[V/m]"; q = a^2/b^2 ;

                     (* Functions *)

cm8[t_] = 2/Pi (b Log[(Sqrt[t+1]+Sqrt[t-q])/Sqrt[q+1]]-
                 a I Log[(Sqrt[t-q]+I a/b Sqrt[t+1])/
                          Sqrt[q+1]])+I (a/Pi Log[t]-b) ;
                                      (* Conformal map *)
                          
af80[t_, sgn8_] = Sqrt[t + sgn8]/Sqrt[t - sgn8 q] ;
ef8[t_] = -I V 1000/b Conjugate[af80[t,1]] ;
            (* Function for the el. field *)
             
t[z_] = Exp[z] ;            

            (* Variables needed in all pictures *)

e = Re[cm8[-t[c]]] ;
frc = FindRoot[e - 1.5 a, {c,0,10}]; c = c/.frc ;
frd = FindRoot[Im[cm8[-t[d]]] + e, {d,-5,-3}]; d = d/.frd ;

Print[ColumnForm[{SequenceForm[
"Distance between horicontal plate parts: a = ", a, " mm;"],
SequenceForm[
"Distance between vertical plate parts: b = ",N[b], " mm;"]}]]

(* Section with black & white graphics *)

Print["Electric field on horicontal part of upper plate "]
ppuphor8 = ParametricPlot[
	  {Re[cm8[-t[-x]]],-mult a/b Re[af80[t[-x],-1]]},
      {x,-c,-.000001}, 
      PlotRange -> {{0,1.02 e},{0,-1.01 mult a/b}},
      AxesLabel -> {"u[mm]",tex1}] ;
 
Print["Electric field on vertical part of upper plate "]
ppupvert8 = ParametricPlot[
      {Im[cm8[-t[-x]]], mult a/b Im[af80[t[-x],-1]]},
      {x,.000001,-d}, 
      PlotRange -> {{-1.02 e,0},{0,mult}},
      AxesLabel -> {"v[mm]",tex1}] ;
      
Print["Electric fieldlines & equipotentials"]
eflandeqpl8 = CartesianMap[cm8[t[#]] &, 
      {d,c,(c - d)/(2 nef-1)}, {0,Pi-.00001,(Pi-.00001)/(2 neq-1)}, 
      PlotRange -> All, 
      AxesLabel -> {"u[mm]","v[mm]"},
      Epilog -> { Thickness[tcb],
                { Line[{{0,-1.02 e},{0,0},{1.02 e,0}}], 
                  Line[{{a,-1.02 e},{a,-b},{1.02 e,-b}}] }} ] ;

Print["Electric field on vertical part of lower plate "]
ppdownvert8 = ParametricPlot[
	  {Im[cm8[t[x]]],-mult a/b Im[af80[t[x],1]]},
      {x,d,Log[q] - .000001}, 
      PlotRange -> {{-1.02 e,-b},{0,3 mult a/b}},
      AxesOrigin -> {-b,0}, 
      AxesLabel -> {"v[mm]",tex1}] ;
       
Print["Electric field on horicontal part of lower plate "]
ppdownhor8 = ParametricPlot[
      {Re[cm8[t[x]]],-mult a/b Re[af80[t[x],1]]},
      {x,Log[q] + .000001,c}, 
      AxesLabel -> {"u[mm]",tex1}, 
      PlotRange -> {{a,1.02 e},{0,-3 mult a/b}}] ;

Print["Electric field of upper plate without units"]
Show[GraphicsArray[{ppupvert8,ppuphor8}, 
GraphicsSpacing -> 0]/.
     (Ticks -> Automatic) -> (Ticks -> None)/.
     (AxesLabel -> {"u[mm]",tex1}) -> 
     (AxesLabel -> {None,None})/.
     (AxesLabel -> {"v[mm]",tex1}) -> 
     (AxesLabel -> {None,None})/.
     (PlotRange -> {{0,1.02 e},{0,-1.01 mult a/b}}) -> 
     (PlotRange -> {{0,1.02 e},{-1.01 mult a/b,mult}})/.
     (PlotRange -> {{-1.02 e,0},{0,mult}}) -> 
     (PlotRange -> {{-1.02 e,0},{-1.01 mult a/b,mult}}),
     Epilog -> {Text["v ->                       u ->",
                 {1,.115 a/b}] } ] ;
      
Print["Electric fieldlines & equipotentials without units"]
Show[eflandeqpl8, Axes -> None] ; 

Print["Electric field of median plate without units"] 
Show[GraphicsArray[{ppdownvert8,ppdownhor8}, 
GraphicsSpacing -> 0]/.
     (Ticks -> Automatic) -> (Ticks -> None)/.
     (AxesLabel -> {"u[mm]",tex1}) -> 
     (AxesLabel -> {None,None})/.
     (AxesLabel -> {"v[mm]",tex1}) ->
     (AxesLabel -> {None,None})/.
     (PlotRange -> {{a,1.02 e},{0,-3 mult a/b}}) -> 
     (PlotRange -> {{a,1.02 e},{-3 mult a/b,3 mult a/b}})/.
     (PlotRange -> {{-1.02 e,-b}, {0,3 mult a/b}}) -> 
     (PlotRange -> {{-1.02 e,-b}, {-3 mult a/b,3 mult a/b}}),
     Epilog -> {Text["v ->                       u ->",
                  {1,.13}] } ] ;

(* Section with colour graphics *)

Print[ColumnForm[{"Absolute value of electric field",
                  "in 3D seen from above upper corner"}]]
absef3D8 = ParametricPlot3D[Evaluate[
 		{Re[cm8[t[x + I y]]],Im[cm8[t[x + I y]]],
 		Abs[ef8[t[x + I y]]] }], 
  		{x,d,2c},{y,0,Pi-.001}, 
  		PlotPoints -> nec,
  		AxesLabel -> {"u[mm]","v[mm]","|E| [V/m]"}, 
  		BoxRatios -> {1,1,1},
  		ViewPoint->{-1.798, 3.333, 2}];

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 plates)"]
eqpl8 = ParametricPlot[ Evaluate[Table[
		{Re[cm8[t[(k-1)/(neq-1) Pi I + x]]], 
		 Im[cm8[t[(k-1)/(neq-1) Pi I + x]]]},
		{k,neq-1}]],{x,d,c},
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [mm]","v [mm]"},
		PlotStyle -> styleq,
		Epilog -> { Thickness[tcb],
                  { Line[{{0,-1.02 e},{0,0},{1.02 e,0}}], 
                    Line[{{a,-1.02 e},{a,-b},{1.02 e,-b}}]}}] ;

Print[ColumnForm[
	  {"Absolute electric fields on equipotential lines as",
      "functions of u (colours match the colours above)"}]]
absefateqplu8 = ParametricPlot[ Evaluate[Table[
		{Re[cm8[t[.99 (k-1)/(neq-1) Pi I + x]]], 
		Abs[ef8[t[.99 (k-1)/(neq-1) Pi I + x]]]},
		{k,neq}]],{x,d,c}, 
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [mm]","|E(u)| [V/m]"},
		PlotStyle -> styleq ];

Print[ColumnForm[
	  {"Absolute electric fields on equipotential lines as", 
	  "functions of v (colours match the colours above)"}]]
absefateqplv8 = ParametricPlot[ Evaluate[Table[
	   {Im[cm8[t[.99 (k-1)/(neq-1) Pi I + x]]], 
		Abs[ef8[t[.99 (k-1)/(neq-1) Pi I + x]]]},
		{k,neq}]],{x,d,c}, 
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"v [mm]","|E(v)| [V/m]"},
		PlotStyle -> styleq ];

Print[ColumnForm[
	  {"Electric fields (v-component) on equipotential lines as", 
	  "functions of u (colours match the colours above)"}]]
efvateqpl8 = ParametricPlot[ Evaluate[Table[
		{Re[cm8[t[.99 (k-1)/(neq-1) Pi I + x]]], 
		 Im[ef8[t[.99 (k-1)/(neq-1) Pi I + x]]]},
		{k,neq}]],{x,d,c}, 
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [mm]","\!\(E\_v\) [V/m]"},
		PlotStyle -> styleq ];

Print[ColumnForm[
	  {"Electric fields (u-component) on equipotential lines as",
	   "functions of v (colours match the colours above)"}]]
efuateqpl8 = ParametricPlot[ Evaluate[Table[
	   {Im[cm8[t[.99 (k-1)/(neq-1) Pi I + x]]], 
		Re[ef8[t[.99 (k-1)/(neq-1) Pi I + x]]]},
		{k,neq}]],{x, d,c}, 
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"v [mm]","\!\(E\_u\) [V/m]"},
		PlotStyle -> styleq ];

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

Print[nef," Electric field lines"]
efl8 = ParametricPlot[ Evaluate[Table[
		{Re[cm8[t[d + (k-1)/(nef-1) dx + I y]]], 
		 Im[cm8[t[d + (k-1)/(nef-1) dx + I y]]]},
		{k,nef}]],{y,0,Pi-.001}, 
		AspectRatio -> Automatic,
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [mm]","v [mm]"},
		PlotStyle -> stylfe,
		Epilog -> { Thickness[tcb],
                  { Line[{{0,-1.02 e},{0,0},{1.02 e,0}}], 
                    Line[{{a,-1.02 e},{a,-b},{1.02 e,-b}}]}}] ;

Print[ColumnForm[
	  {"Absolute electric fields on field lines as",
	  "functions of u (colours match the colours above)"}]]
absefateflu8 = ParametricPlot[ Evaluate[Table[
		{Re[cm8[t[d + (k-1)/(nef-1) dx + I y]]], 
		Abs[ef8[t[d + (k-1)/(nef-1) dx + I y]]]},
		{k,nef}]],{y,0,Pi-.001}, 
		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)"}]]
absefateflv8 = ParametricPlot[ Evaluate[Table[
		{Im[cm8[t[d + (k-1)/(nef-1) dx + I y]]], 
		Abs[ef8[t[d + (k-1)/(nef-1) dx + I y]]]},
		{k,nef}]],{y,0,Pi-.001}, 
		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)"}]]
efvatefl8 = ParametricPlot[ Evaluate[Table[
		{Re[cm8[t[d + (k-1)/(nef-1) dx + I y]]], 
		 Im[ef8[t[d + (k-1)/(nef-1) dx + I y]]]},
		{k,nef}]],{y,0,Pi-.001}, 
		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)"}]]
efuatefl8  = ParametricPlot[ Evaluate[Table[
	   {Im[cm8[t[d + (k-1)/(nef-1) dx + I y]]], 
		Re[ef8[t[d + (k-1)/(nef-1) dx + I y]]]},
		{k,nef}]],{y,0,Pi-.001}, 
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"v [mm]","\!\(E\_u\) [V/m]"},
		PlotStyle -> stylfe] ;

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

up = 1.5; vp = - .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. eflandeqpl8 *) 

Which[vp > 0,     Print[vp " = vp > 0  is unphysical.",⚠],
      up < 0 || (up > a && vp < -b),     
      Print[up " = up < 0 or", up " = up > ",a," and" 
      vp " = vp > ",-b,"is unphysical.",⚠],
      True,

frefp = FindRoot[cm8[t[efp]] - (up + I vp), 
        {efp, Sign[up - a] + I (Pi -.1)}];
efp = efp /. frefp; imefp = N[Im[efp]];
error = Abs[N[cm8[t[efp]] - (up + I vp)]];
				
Which[0 > imefp, 	 Print["The calculated point doesn`t fulfill the requirements!",\[FreakedSmiley]], 
     imefp > N[Pi], Print["The calculated point doesn`t fulfill the requirements!",\[FreakedSmiley]],
     error > 10^-6, Print["The error for the calculated point is too big !",☹], 
     True, Print["The calculated point fulfills the requirements.",☺];
     
valef = ef8[t[efp]]//N;
Print[
SequenceForm["At point (up,vp) = (",up ,", ",vp ,") mm : ",
"Potential = ",N[V Im[efp]/Pi], " V;"]];
Print[SequenceForm["Field (\!\(E\_u\),\!\(E\_v\)) = (", 
Re[valef],", ", Im[valef],") V/m "]];      


caleqpl8 = ParametricPlot[Evaluate[
	  {Re[cm8[t[I Im[efp] + s]]], Im[cm8[t[I Im[efp] + s]]] }], 
	  {s,d,c}, 
	  PlotStyle -> Dashing[{.01}], 
	  DisplayFunction -> Identity];
calefl8 = ParametricPlot[Evaluate[
	  {Re[cm8[t[Re[efp] + I s]]], Im[cm8[t[Re[efp] + I s]]] }],
	  {s,0,Pi-.001}, 
	  DisplayFunction -> Identity, 
	  PlotRange -> All];
	  
Show[caleqpl8,calefl8, 
      DisplayFunction -> $DisplayFunction, 
	  PlotRange  -> {{0,1.02e},{-1.02e,0}},
	  PlotRegion -> {{.01,.99},{.01,0.99}},
	  AxesLabel -> {"u [mm]","v[mm]"},
 	  Epilog -> { Thickness[tcb],
 	  			{ Line[{{0,-1.02 e},{0,0},{1.02 e,0}}], 
                  Line[{{a,-1.02 e},{a,-b},{1.02 e,-b}}] }} ]; ]]