(* thicksem. nb *)

Needs["Graphics`ComplexMap`"]

Clear[V,a,b,c,d,e,f,g,p,ssp,dx,sp,frc,frd,
               nec,nef,neq,tcf,tcb,mult,tex1];
Off[General::spell1]

                   	(* Input data *)

V = 1  (* Potential on upper plate in Volt, 
		       potential on lower plate (= median plane) = 0 *);
a = 3  (* Distance between plates in mm : 
         for 10^-4 < a < 3 10^5 *); 
b = .5 (* Thickness of the upper plate in mm:
	         for ratios 6 < a/b < 700 
		     depending on the value of a *) ;     
          
nec = 30 (* Number of lines used in 3D-image of 
            el. field`s absolute value: absef3D4 *);
            
nef = 11 (* Number of field lines between plates;
            variable needed in the pictures: 
            eflandeqpl4, efl4, absefateflu4, absefateflv4, 
	           efuatefl4, efvatefl4 *);
            
neq = 7 (* Number of equipotential lines between plates;
           variable needed in the pictures: 
           eflandeqpl4, eqpll4, absefateqplu4, absefataqplv4, 
	          efuateqpl4, efvateqpl4 *);

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]";           
 
       (* Calculation of parameter depending on ratio a/b *)
                   
ssp=Solve[(sp-1)^2/(2 sp)-N[b]/a==0,sp];
sp=N[sp/.ssp];
If[First[sp]>Last[sp],sp=First[sp],sp=Last[sp]];p=sp^2;

                       (* Functions *) 
                               
cm4[t_] = (a (Sqrt[t + 1] Sqrt[t + p] + 
        (p + 1) Log[Sqrt[t + 1] + Sqrt[t + p]] - 
        2 sp Log[sp Sqrt[t + 1] + Sqrt[t + p]] + 
        sp Log[t]))/(Pi sp) - (b Log[p - 1])/Pi ;
	                          (* Conformal map *)
	                                         
ef4[t_]=sp/(Sqrt[t+1] Sqrt[t+p]); (* Function of el. field*)
t[z_]=Exp[z];

             (* Variables needed in all pictures *)

e=cm4[t[c]];f=cm4[t[d]];g=Im[ef4[-t[.5 Log[p]]]];
frc=FindRoot[e+1.5 (a+b),{c,-6,-5}];c=Re[c]/.frc;
frd=FindRoot[e+N[a+b]I-cm4[-t[d]],
     {d,Log[p]+.5,5},MaxIterations -> 300];d=Re[d]/.frd;

Print[ColumnForm[
   {SequenceForm["Distance between finitely thick plate and median plane: a = " ,a, " mm;"],
    SequenceForm["Thickness of thick plate: b = " ,b  ," mm."]}]]

(* Section with black & white graphics *)

Print["Electric field on top side of upper plate "]
ppuptop4=ParametricPlot[{Re[cm4[-t[-x]]],-mult Re[ef4[-t[-x]]]},
        {x,-d,-Log[p]-.00001},
		PlotRange -> {{e,0},{0,3 mult}},
		AxesLabel -> {"u[mm]",tex1}];

(* Calculations for a later picture *) ppuptop4inv=ParametricPlot[{-Re[cm4[-t[-x]]],-mult Re[ef4[-t[-x]]]},
       {x,-d,-Log[p]-.00001},
		PlotRange -> {{-e,0},3 mult {-1,1}},
		DisplayFunction->Identity];

Print["Electric field on front of upper plate "]
ppupfront4=ParametricPlot[{Im[cm4[-t[-x]]],-mult Im[ef4[-t[-x]]]},
        {x,-Log[p]+.00001,-.00001},
		PlotRange -> {{a,a+b},{0,7 mult}},
		AxesLabel -> {"v[mm]",tex1}];

Print["Electric field on bottom side of upper plate "]
ppupbot4=ParametricPlot[
       {Re[cm4[-t[-x]]],- mult Re[ef4[-t[-x]]]},{x,.000001,-c},
		PlotRange -> {{e,0},{0,-3 mult }},
		AxesLabel -> {"u[mm]",tex1}];

Print["Electric fieldlines & equipotentials "]
eflandeqpl4=CartesianMap[cm4[t[#1]]&,{c,d,(d-c)/(nef-1)},
        {0,Pi,Pi/(neq-1)},
		PlotRange -> All,
		AxesLabel -> {"u[mm]","v[mm]"},
		Epilog->{Thickness[.011],
				{Line[{{1.02 f,0},{1.02 e,0}}],
          Line[{{1.02 e,a},{0,a},{0,a+b},{1.02 e,a+b}}]}}];

Print["Electric field on median plate "]
ppdown4=ParametricPlot[{cm4[t[x]],-mult ef4[t[x]]},{x,c,d},
		PlotRange -> {{e,f},{-mult,0}},
		AxesLabel -> {"u[mm]",tex1}];

Print["Electric field of upper plate without units"]
Show[GraphicsArray[{ppupbot4,Show[ppupfront4,
     PlotRange -> {{a,a+b},3 mult {-1,1}},
     Ticks -> None,
	 AxesLabel -> None,
	 DisplayFunction -> Identity],
		            ppuptop4inv},
     GraphicsSpacing -> 0  ]/.
	 (Ticks->Automatic) ->(Ticks->None)/.                                          
	 (AxesLabel->{"u[mm]",tex1})->
	 (AxesLabel->{None,None})/.      
	 (AxesLabel->{"v[mm]",tex1})->
	 (AxesLabel->{None,None})/.                 
	 (PlotRange->{{e,0},{0,-3 mult}})-> 
	 (PlotRange->{{e,0},3mult {-1,1}}),
	 Epilog->{Text["Bottom u ->     Front v ->          Top <- u",{1.5,0}]}];

Print["Electric fieldlines & equipotentials without units"]
Show[eflandeqpl4,Axes -> None];
	
Print["Electric field of median plate without units"] 	
Show[ppdown4,Ticks -> None,AxesLabel -> {"u","E"}];

(* Section with colour graphics *)

Print[ColumnForm[
	  {"Absolute value of electric field",
	   "in 3D seen from above upper plate"}]]
absef3D4 = ParametricPlot3D[Evaluate[
 		{Re[cm4[t[x + I y]]],Im[cm4[t[x + I y]]],
 	     mult Abs[ef4[t[x + I y]]]}],{x,c,d},{y,0,Pi},
  		AxesLabel -> {"u[mm]","v[mm]","|E| [V/m]"}, 
  		BoxRatios -> {1,1.5,1}, 
  		PlotPoints -> nec,
		ViewPoint -> {-1.07, 2.3, 1.9}];

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

Print[neq," Equipotential lines (including plates)"]
eqpl4 = ParametricPlot[ Evaluate[Table[
     {Re[cm4[t[(k-1)/(neq-1)Pi I + x]]],
      Im[cm4[t[(k-1)/(neq-1)Pi I + x]]]},
		{k,neq}]],{x,c,d}, 
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [mm]","v [mm]"},
		PlotStyle -> styleq,
     	Epilog -> {Thickness[tcb],
                 { Line[{{1.02f,0},{1.02e,0}}], 
                   Line[{{1.02e,a},{0,a},{0,a+b},{1.02e,a+b}}] }} ];

Print[ColumnForm[
	  {"Absolute electric fields on equipotential lines as",
      "functions of u (colours match the colours above)"}]]
absefateqplu4 = ParametricPlot[ Evaluate[Table[
	{Re[cm4[t[(k-1)/(neq-1)Pi I + x]]],
    mult Abs[ef4[t[(k-1)/(neq-1)Pi I + x]]]},
		{k,neq}]],{x,c,d},
		PlotRange -> {0,2.5 mult}, 
		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)"}]]
absefateqplv4 = ParametricPlot[ Evaluate[Table[
	{Im[cm4[t[(k-1)/(neq-1)Pi I + x]]],
    mult Abs[ef4[t[(k-1)/(neq-1)Pi I + x]]]},
		{k,neq}]],{x,c,d},
		PlotRange -> {0,2.5 mult}, 
		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)"}]]
efvateqpl4 = ParametricPlot[ Evaluate[Table[
	 {Re[cm4[t[(k-1)/(neq-1)Pi I + x]]],
	  -mult Re[ef4[t[(k-1)/(neq-1)Pi I + x]]]},
		{k,neq}]],{x,c,d}, 
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [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)"}]]
efuateqpl4 = ParametricPlot[ Evaluate[Table[
	 {Im[cm4[t[(k-1)/(neq-1)Pi I + x]]],
	  -mult Im[ef4[t[(k-1)/(neq-1)Pi I + x]]]},
		{k,neq}]],{x,c,d}, 
		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"]
efl4 = ParametricPlot[Evaluate[Table[
          {Re[cm4[t[c + ((k - 1) dx)/(nef - 1) + I y]]], 
           Im[cm4[t[c + ((k - 1) dx)/(nef - 1) + I y]]]}, 
          {k, nef}]], {y, 0, Pi}, 
    PlotRegion -> {{0.01, 0.99}, {0.01, 0.99}}, 
    AxesLabel -> {"u [mm]", "v [mm]"}, 
    PlotStyle -> stylfe, 
    Epilog -> {Thickness[tcb], 
         {Line[{{1.02 f, 0}, {1.02 e, 0}}], 
          Line[{{1.02 e, a}, {0, a}, {0, a + b}, {1.02 e, a + b}}]}}];

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

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

up = 8; vp = 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. eflandeqpl4 *) 
				
Which[vp < 0, Print[vp " = vp < 0 is unphysical.",⚠],
      a < vp < a+b && up < 0, 
      Print[vp " = ",a," < vp < ",a+b," and 
            ",up," = up < 0 is unphysical.",⚠],
      True,
      
frefp = FindRoot[cm4[t[efp]] - (up + I vp), 
            {efp,  1.2  Sign[vp - a] + I Pi/2 }];
efp = efp /. frefp; imefp = N[Im[efp]];
error = Abs[N[cm4[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 = -mult ef4[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\)) = (", 
Im[valef],", ", Re[valef],") V/m "]]; 

caleqpl = ParametricPlot[Evaluate[
	  {Re[cm4[t[I Im[efp] + s]]], Im[cm4[t[I Im[efp] + s]]]}], 
	  {s,c,d}, 
	  PlotStyle -> Dashing[{.01}], 
	  DisplayFunction -> Identity];
	  
calefl = ParametricPlot[Evaluate[
	  {Re[cm4[t[Re[efp] + I s]]], Im[cm4[t[Re[efp] + I s]]]}],
	  {s,0,Pi}, DisplayFunction -> Identity];
	  
Show[caleqpl,calefl,
      DisplayFunction -> $DisplayFunction, 
      PlotRange -> All,
	  AspectRatio -> Automatic, AxesLabel -> {"u [mm]","v [mm]"},
 	  Epilog -> { Thickness[tcb],
                 {Line[{{1.02 f,0},{1.02 e,0}}],
                  Line[{{1.02 e,a},{0,a},{0,a+b},{1.02e,a+b}}]}}]; ]] 
                  
On[Genal::spell1]