(* aperture.nb *)

Needs["Graphics`ComplexMap`"]

Clear[V,a,b,c,d,e,f,frb,fre,frf,chs,sgn9,
      nec,nef,neq,tcb,tcf,mult,tex1]
Off[General::spell1]

	                   (* Input data *)
	                        	
V = 1 (* Potential on the upper plate in Volt  
         potential on lower plate (= median plane) = 0 *);	
         		
a = 3 (* Shortest distance between the plates in mm :
         for 7 10^-5 < a < 10^5 *);

nec = 20 (* Number of lines used in 3D-image of 
            el. field`s absolute value: absef3D9 *);
            
nef = 11 (* Number of field lines between plates;
            variable needed in the pictures: 
            eflandeqpl9, efl9, absefateflu9, absefateflv9, 
            efuatefl9, efvatefl9 *);
            
neq = 7 (* Number of equipotential lines between plates;
            variable needed in the pictures: 
            eflandeqpl9, eqpll9, absefateqplu9, 
            absefataqplv9, efuateqpl9, efvateqpl9 *);
	
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]";   

                       (* Functions *)

cm9[z_, sgn9_] = 2 Exp[z/2] + sgn9 a^2/8 Exp[-z/2] ; 
                                 (* Conformal map *)
                                 
ef9[z_] = 1000 V/(Exp[z/2] + a^2/16 Exp[-z/2])/Pi ;
                        (* Function of el. field *)
                        
If[a < 1.7008433 10^2, chs = 1, chs = -1] ; 

             (* Variables needed in all pictures *)
             
c = cm9[b,1] ; d = cm9[b,-1] ;
frb = FindRoot[d - 3 a, {b,1,10}]; b = b/.frb ;
fre = FindRoot[cm9[e,-1] + 3 a,{e,-5,-15}]; e = e/.fre ; 
frf = FindRoot[cm9[-f,1] - c,{f,-10,-5}]; f = f/.frf ;

Print[SequenceForm[
"Distance between infinitely thin plate and median plane: a = ", a, " mm."]]

(* Section with black & white graphics *)

Print["Electric field on left side of upper plate "]
ppupleft9 = ParametricPlot[{cm9[-x,1],-2 chs mult/cm9[-x,-1]},
      {x,f,-2 Log[a/4] - chs .1}, 
      PlotRange -> {All,{0,-6 mult}},
      AxesOrigin -> {a,0}, 
      AxesLabel -> {"v[mm]",tex1}] ;
     
(* Calculations for a later picture *)    
ppupleft9inv = ParametricPlot[{-cm9[-x,1],-2 chs mult/cm9[-x,-1]},
         {x,f,-2 Log[a/4] - chs .01}, 
         PlotRange -> {{-a,-c},6 mult{-1,1}}, 
         AxesOrigin -> {-a,0}, 
         DisplayFunction -> Identity] ;
      
Print["Electric field on right side of upper plate "]
ppupright9 = ParametricPlot[{cm9[-x,1], -2 chs mult/cm9[-x,-1]},
      {x,-2 Log[a/4] + chs .01,-f - 4 Log[a/4]}, 
      PlotRange -> {{a,c},{0,6 mult}}, 
      AxesOrigin -> {a,0} , 
      AxesLabel -> {"v[mm]",tex1}] ;
     
Print["Electric field lines & equipotentials "]
eflandeqpl9 = CartesianMap[cm9[#,-1] &,
      {b,e,(e-b)/(2 nef-1)},{0,Pi,Pi/(neq-1)},
      PlotRange -> All, 
      AxesLabel -> {"u[mm]","v[mm]"}, 
      Epilog -> {Thickness[tcb],
      			{Line[{{0,1.05 a},{0,1.02 c}}],
                 Line[{{1.02 d,0},{-1.02 d,0}}]}}] ;   
             
Print["Electric field on median plate "]
ppdown9 = ParametricPlot[{cm9[x,-1], -2 mult/cm9[x,1]},{x,b,e},
      PlotRange -> {{d,-d},{0,-.7 mult}}, 
      AxesLabel -> {"u[mm]",tex1}];

Print["Electric field of upper plate without units"]
Show[GraphicsArray[{ppupleft9inv,ppupright9}, 
GraphicsSpacing -> 0] /.
     (Ticks -> Automatic) -> (Ticks -> None) /.
     (AxesLabel -> {"v[mm]",tex1}) -> (AxesLabel -> None) /.
     (PlotRange -> {{a,c},{0,6 mult}}) -> 
     (PlotRange -> {{a,c},6 mult{-1,1}}), 
     Epilog -> { Text["<- v left                right v ->",
     			      {1,.2}]},
     DisplayFunction -> $DisplayFunction ] ; 
     
Print["Electric fieldlines & equipotentials without units"]
Show[eflandeqpl9, Axes -> None] ; 

Print["Electric field of median plate without units"] 
Show[ppdown9,
     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 above upper plate"}]]
absef3D9 = ParametricPlot3D[Evaluate[
 		{Re[cm9[x + I y,-1]],Im[cm9[x + I y,-1]],
 		Abs[ef9[x + I y]]}],
  		{x,e,b},{y,0,Pi-.00001},
  		PlotPoints -> nec,
  		AxesLabel -> {"u[mm]","v[mm]","|E| [V/m]"}, 
  		BoxRatios -> {1,1,1}];

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

Print[ColumnForm[
	  {"Absolute electric fields on equipotential lines as",
      "functions of u (colours match the colours above)"}]]
absefateqplu9 = ParametricPlot[ Evaluate[Table[
		{Re[cm9[(k-1)/(neq-1) Pi I + x,-1]], 
		Abs[ef9[(k-1)/(neq-1) Pi I + x]]},
		{k,neq}]],{x,b,e}, 
		PlotRange -> {All,{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)"}]]
absefateqplv9 = ParametricPlot[ Evaluate[Table[
		{Im[cm9[(k-1)/(neq-1) Pi I + x,-1]], 
		Abs[ef9[(k-1)/(neq-1) Pi I + x]]},
		{k,neq}]],{x,b,e}, 
		PlotRange -> {0,3.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)"}]]
efvateqpl9 = ParametricPlot[ Evaluate[Table[
		{Re[cm9[(k-1)/(neq-1) Pi I + x,-1]], 
		 -Re[ef9[(k-1)/(neq-1) Pi I + x]]},
		{k,neq}]],{x,b,e}, 
		PlotRange -> {All,{0,-2 mult}},
		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)"}]]
efuateqpl9 = ParametricPlot[ Evaluate[Table[
		{Im[cm9[(k-1)/(neq-1) Pi I + x,-1]], 
		 -Im[ef9[(k-1)/(neq-1) Pi I + x]]},
		{k,neq}]],{x,b,e}, 
		PlotRange -> {-2 mult,2 mult},
		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]}, Reverse[hli], 
             {GrayLevel[0]},hli,{GrayLevel[0]}] ;
hlt = 	Table[Thickness[tcf],{k,2 nef-1}] ;
stylfe = Transpose[{hlc,hlt}] ;
dx  = 	e - b;

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

Print[ColumnForm[
	  {"Absolute electric fields on field lines as",
	  "functions of u (colours match the colours above)"}]]
absefateflu9 = ParametricPlot[ Evaluate[Table[
		{Re[cm9[b + (k-1)/(nef-1)/2 dx  + I y,-1]], 
		Abs[ef9[b + (k-1)/(nef-1)/2 dx  + I y]]},
		{k,2 nef-1}]],{y,0,Pi-.001}, 
		PlotRange -> {0,2.5 mult},
		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)"}]]
absefateflv = ParametricPlot[ Evaluate[Table[
		{Im[cm9[b + (k-1)/(nef-1)/2 dx  + I y,-1]], 
		Abs[ef9[b + (k-1)/(nef-1)/2 dx  + I y]]},
		{k,nef}]],{y,0,Pi-.001}, 
		PlotRange -> {0,3 mult},
		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)"}]]
efvatefl9 = ParametricPlot[ Evaluate[Table[
		{Re[cm9[b + (k-1)/(nef-1)/2 dx  + I y,-1]], 
		 -Re[ef9[b + (k-1)/(nef-1)/2 dx  + I y]]},
		{k,2 nef-1}]],{y,0,Pi-.001}, 
		PlotRange -> {0,-2 mult},
		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)"}]]
efuatefl9 = ParametricPlot[ Evaluate[Table[
		{Im[cm9[b + (k-1)/(nef-1)/2 dx  + I y,-1]], 
		 -Im[ef9[b + (k-1)/(nef-1)/2 dx  + I y]]},
		{k,2 nef-1}]],{y,0,Pi-.001}, 
		PlotRange -> {2 mult,-2 mult},
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"v [mm]","\!\(E\_u\) [V/m]"},
		PlotStyle -> stylfe] ;

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

up = .2; vp = .2;
(* 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. eflandeqpl9; but calculations below fail for points near to the corners. *) 

Which[vp < 0,     Print[vp " = vp < 0  is unphysical.",⚠],
      True,
      
frefp = FindRoot[cm9[efp,-1] - (up + I vp), 
         {efp, b - .1 + .1 I, e + .1 + (Pi - .1) I} ];

efp = efp /. frefp; imefp = N[Im[efp]];
error = Abs[N[cm9[efp,-1] - (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 = ef9[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 "]]; 


caleqpl9 = ParametricPlot[Evaluate[
	  {Re[cm9[I Im[efp] + s,-1]], 
	   Im[cm9[I Im[efp] + s,-1]]}], {s,b,e}, 
	  PlotStyle -> {{Thickness[tcf], Dashing[{.01}]}}, 
	  DisplayFunction -> Identity];
	  
calefl9 = ParametricPlot[Evaluate[
	  {Re[cm9[Re[efp] + I s,-1]], 
	   Im[cm9[Re[efp] + I s,-1]]}],
	  {s,0,Pi-.001}, DisplayFunction -> Identity,
	  PlotStyle -> Thickness[tcf]];
	  
Show[caleqpl9,calefl9, 
     DisplayFunction -> $DisplayFunction,
	 AxesLabel -> {"u","v"},
	 PlotRange -> All,
	 Epilog -> {Thickness[tcb],
               	  {Line[{{0,1.05 a},{0,1.02 c}}] ,
                   Line[{{1.02 d,0},{-1.02 d,0}}]}}]; ]]