(* finite.nb*)

Needs["Graphics`ComplexMap`"]

Clear[V,a,b,m,sn,elk,mprime,frmprime,
      nef,neq,nec,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 (* Twice the distance betweem the plates in mm :
         for 2 10^-5 < a < 10^5 *); 
b = 2 (* Length of upper plate in mm :
         for 10^-3 a < b < 8.75 a *);

nec = 20 (* Number of lines used in 3D-image of 
            el. field`s absolute value: absef3D11 *);

nef = 11 (* Number of field lines between plates;
            variable needed in the pictures: 
            eflandeqpl11, efl11, absefateflu11, absefateflv11, 
            efuatefl11, efvatefl11 *);

neq = 11  (* Number of equipotential lines between plates;
            variable needed in the pictures: 
            eflandeqpl11, eqpll11, absefateqplu11, 
            absefataqplv11, efuateqpl11, efvateqpl11 *);
            
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 the parameter for elliptic integrals *)

af110 = 2/Pi (EllipticK[mprime] EllipticE[ArcSin[1/Sqrt[mprime]*
              Sqrt[1 - EllipticE[mprime]/EllipticK[mprime]]],mprime] - 
              EllipticE[mprime] EllipticF[ArcSin[1/Sqrt[mprime]* 
              Sqrt[1 - EllipticE[mprime]/EllipticK[mprime]]],mprime]);

frmprime = FindRoot[af110 - b/a ,{mprime,.99,.999}]; 
mprime = mprime/.frmprime ; m = 1 - mprime ;

                            (* Global variables *)
                  
mult = V 1000/a; tex1 = "E[V/m]"; elk = Pi/EllipticK[m];

                               (* Functions *) 

sn = InverseJacobiSN[
      Sqrt[(1 - m EllipticK[mprime]/EllipticE[mprime]) / mprime], mprime] ; 

cm11[z_] = a/Pi ((EllipticE[mprime] - EllipticK[mprime]) z +
                   EllipticK[mprime] EllipticE[JacobiAmplitude[z,m],m]) ;
                                                      (* Conformal map *)      
                   
af111[z_] = 1/(EllipticE[mprime] - 
            m EllipticK[mprime] JacobiND[z,mprime]^2) ;

ef11[z_] = -N[elk mult /Conjugate[EllipticE[mprime]-
             m EllipticK[mprime] JacobiSN[z,m]^2] //Together//ExpandAll]; 
                                              (* Function of el. field *)
                                          
Print[ColumnForm[
{SequenceForm["Distance between finite plate and median plane: a/2 = ",.5 a, " mm."],
SequenceForm["Length of finite plate: b = ", b, " mm."]}]]

(* Section with black & white graphics *)            

Print["Electric field on top side of upper plate "]
ppuptop11 = ParametricPlot[{
           {Im[cm11[EllipticK[m] - 10^-5 + I y]],-elk mult af111[y]},
           {-Im[cm11[EllipticK[m] - 10^-5 + I y]],-elk mult af111[y]}},
           {y,sn + 10^-6, EllipticK[mprime] - 10^-6}, 
           PlotRange -> {{-1.02 b/2, 1.02 b/2},{0,4 elk mult}},
           AxesLabel -> {"u[mm]",tex1}] ;                             


Print["Electric field on bottom side of upper plate "]
ppupbottom11 = ParametricPlot[{
              {Im[cm11[EllipticK[m] - 10^-5 + I y]],-elk mult af111[y]},
              {-Im[cm11[EllipticK[m] - 10^-5 + I y]],-elk mult af111[y]}},
              {y,0,sn - 10^-6}, 
              PlotRange -> {{-1.02 b/2, 1.02 b/2},{0,-4 elk mult}},
              AxesLabel -> {"u[mm]",tex1}] ;          
                              
(* Calculations for a later picture *)  
ppupbottom11inv = ParametricPlot[{
           {Im[cm11[EllipticK[m] - 10^-5 + I y]],elk mult af111[y]},
           {-Im[cm11[EllipticK[m] - 10^-5 + I y]],elk mult af111[y]}},
           {y,0,sn - 10^-6}, 
           PlotRange -> {{-1.02 b/2, 1.02 b/2},{0,elk mult 4}},
           PlotStyle -> Dashing[{.02}], 
           DisplayFunction -> Identity] ;
               
Print["Electric field lines & equipotentials "]
eflandeqpl11 = CartesianMap[I cm11[#] &, 
        {0, EllipticK[m] - .0001, (EllipticK[m]-.0001)/(neq-1)},
        {-EllipticK[mprime] + .0001, EllipticK[mprime] - .0001,
        (2 EllipticK[mprime]-.0002)/(2 nef-1)},
       PlotPoints -> 25, 
       PlotRange -> {{-2 b,2 b},{-.2 a,2 a}},
       AxesLabel -> {"u[mm]","v[mm]"},
       Epilog -> { Thickness[tcb],
                  { Line[{{-b/2,a/2},{b/2,a/2}}], 
                    Line[{{-2b,0},{2b,0}}] }} ] ; 
                    
Print["Electric field on lower plate "]
ppdown11 = ParametricPlot[
           {Im[cm11[I y]],Re[ef11[I y]]},
           {y,-EllipticK[mprime],EllipticK[mprime]},
           PlotRange -> {{-2b,2b},{0,-elk mult}},
           PlotRegion -> {{.01,.99},{.01,.99}},
           AxesLabel -> {"u[mm]",tex1}];

Print[ColumnForm[
      {"Electric field of upper plate without units:",
       "Dashed line = bottom; Full line = top."}]]
Show[ppuptop11,ppupbottom11inv, 
         Ticks -> None,
         AxesLabel -> {"u","E"},
         PlotRegion -> {{.01,.99},{.01,.99}},
         DisplayFunction -> $DisplayFunction ] ; 

Print["Electric fieldlines & equipotentials without units"]
Show[eflandeqpl11, Axes -> None] ; 

Print["Electric field of lower plate without units"]
Show[ppdown11, Ticks -> None, AxesLabel -> {"u","E"}] ;

Print[ColumnForm[{"Absolute value of electric field",
                  "in 3D seen from above upper plate"}]]
absef3D11 = ParametricPlot3D[Evaluate[
 		{Re[cm11[x + I y]],Im[cm11[x + I y]],
 		 Abs[ef11[x + I y]]}],
  		{x,0,EllipticK[m]-.0001},
  		{y,-EllipticK[mprime]+.0001,EllipticK[mprime]-.0001},
  		PlotRange -> {0, 2 Abs[mult elk]},
  		AxesLabel -> {"v[mm]","u[mm]","|E| {V/m]"}, 
  		BoxRatios -> {1,1,1}, 
  		PlotPoints -> nec,
  		ViewPoint -> {2,.5,0}];

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

Print[neq," Equipotential lines (including plates)"]
eqpl11 = ParametricPlot[ Evaluate[Table[
		{Im[cm11[.9999(neq-k)/(neq-1) EllipticK[m] + I y]], 
		 Re[cm11[.9999(neq-k)/(neq-1) EllipticK[m] + I y]]},
		{k,neq}]],{y,-EllipticK[mprime]+.0001,EllipticK[mprime]-.0001},
		PlotRange -> {{-2.1b,2.1b},{-.2a,2a}},
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [mm]","v [mm]"},
		PlotStyle -> styleq,
		Epilog -> { Thickness[tcb],
                   Line[{{-b/2,a/2},{b/2,a/2}}],
                   Line[{{-2b,0},{2b,0}}]}];

Print[ColumnForm[
	  {"Absolute electric fields on equipotential lines as",
      "functions of u (colours match the colours above)"}]]
absefateqplu11 = ParametricPlot[ Evaluate[Table[
		{Im[cm11[.9999(k-1)/(neq-1) EllipticK[m] + I y]],
		 Abs[ef11[.9999(k-1)/(neq-1) EllipticK[m] + I y]]},
	    {k,neq,1,-1}]],
	    {y,-EllipticK[mprime]+.0001,EllipticK[mprime]-.0001}, 
		(*PlotRange -> {A,{}},*)
		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 u (colours match the colours above)"}]]
absefateqplv11 = ParametricPlot[ Evaluate[Table[
		{Re[cm11[.9999(k-1)/(neq-1) EllipticK[m] + I y]],
		 Abs[ef11[.9999(k-1)/(neq-1) EllipticK[m] + I y]]},
	    {k,neq,1,-1}]],
	    {y,-EllipticK[mprime]+.0001,EllipticK[mprime]-.0001}, 
		(*PlotRange -> {A,{}},*)
		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)"}]]
efvateqpl11 = ParametricPlot[Evaluate[Table[
     {Im[cm11[0.9999(k - 1)/(neq - 1)EllipticK[m] + I y]], 
      Re[ef11[0.9999(k - 1)/(neq - 1)EllipticK[m] + I y]]}, 
     {k, neq,1,-1}]], 
     {y, -EllipticK[mprime] + 0.0001, EllipticK[mprime] - 0.0001}, 
     PlotRegion -> {{0.01, 0.99}, {0.01, 0.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)"}]]
efuateqpl11 = ParametricPlot[ Evaluate[Table[
		{Re[cm11[.9999(k-1)/(neq-1) EllipticK[m] + I y]],
		 Im[ef11[.9999(k-1)/(neq-1) EllipticK[m] + I y]]},
		{k,neq,1,-1}]],
		{y,-EllipticK[mprime]+.0001,EllipticK[mprime]-.0001}, 
		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]}, 
        Reverse[hli], {GrayLevel[0]}]; 
hlt = 	Table[Thickness[tcf],{k,2nef-1}];
stylef = Transpose[{hlc,hlt}] ;

Print[2 nef-2," Electric field lines"]
efl11 = ParametricPlot[ Evaluate[Table[
		{Im[cm11[.9999(k-1)/(nef-1) EllipticK[mprime] I + x]], 
		 Re[cm11[.9999(k-1)/(nef-1) EllipticK[mprime] I + x]]},
		{k,-nef+2,nef}]],{x,0,EllipticK[m] - .0001},
		PlotRange -> {{-2b,2b},{-.2a,2a}}, 
		PlotRegion -> {{.01,.99},{.01,.99}},
		AxesLabel -> {"u [mm]","v [mm]"},
		PlotStyle -> stylef,
		Epilog -> { Thickness[tcb],
                    Line[{{-b/2,a/2},{b/2,a/2}}], 
                    Line[{{-2b,0},{2b,0}}]}] ;

\(Print[
      ColumnForm[

        	\ \ {"Absolute electric fields on field lines as", 
	\ \ \ \ 
          "functions of u (colours match the colours above)"}]]\)\(\(absefateflu11\  = \ 
      ParametricPlot[\ 
        Evaluate[
          Table[
		{
              Im[cm11[ .9999 \((k - 1)\)/\((nef - 1)\)\ EllipticK[mprime]\ I
                      \  + x]], 
		\ 
              Abs[ef11[ .9999 \((k - 1)\)/\((nef - 1)\)\ EllipticK[mprime] I
                      \  + x]]}, 
		{k, \(-nef\) + 2, nef}]], {x, 0, 
          EllipticK[m] -  .0001}, 
		
        PlotRange\  -> \ {{\(-2\) b, 2  b}, {0, 3\ \ elk\ mult}}, 
		
        PlotRegion\  -> \ {{ .01,  .99}, { .01,  .99}}, 	
		
        AxesLabel\  -> \ {"u [mm]", \ "|E(u)| [V/m]"}, \ \ 
		
        PlotStyle\  -> \ stylef]; \)\)

\(Print[
      ColumnForm[

        	\ \ {"Absolute electric fields on field lines as", 
	\ \ \ \ 
          "functions of u (colours match the colours above)"}]]\)\(\(absefateflv11\  = \ 
      ParametricPlot[\ 
        Evaluate[
          Table[
		{
              Re[cm11[ .9999 \((k - 1)\)/\((nef - 1)\)\ EllipticK[mprime]\ I
                      \  + x]], 
		\ 
              Abs[ef11[ .9999 \((k - 1)\)/\((nef - 1)\)\ EllipticK[mprime] I
                      \  + x]]}, 
		{k, nef, 1, \(-1\)}]], {x, 0, 
          EllipticK[m] -  .0001}, 
		
        PlotRegion\  -> \ {{ .01,  .99}, { .01,  .99}}, 	
		
        AxesLabel\  -> \ {"v [mm]", \ "|E(v)| [V/m]"}, 
		
        PlotStyle\  -> \ stylef]; \)\)

\(Print[
      ColumnForm[

        	\ \ {"Electric fields (v-component) on field lines as", \ 

          \ \ 	\ \ "functions of u (colours match the colours above)"}]]
      \)\(\(efvatefl11\  = \ 
      ParametricPlot[\ 
        Evaluate[
          Table[
		{
              Im[cm11[ .9999 \((k - 1)\)/\((nef - 1)\)\ EllipticK[mprime]\ I
                      \  + x]], 
		\ 
              Re[ef11[ .9999 \((k - 1)\)/\((nef - 1)\)\ EllipticK[mprime]\ I
                      \  + x]]}, 
		{k, \(-nef\) + 2, nef}]], {x, 0, 
          EllipticK[m] -  .0001}, \ 
		
        PlotRegion\  -> \ {{ .01,  .99}, { .01,  .99}}, 	
		
        AxesLabel\  -> \ {"u [mm]", \ \*""\!\(E\_v\)[V/m]""}, 
		
        PlotStyle\  -> \ stylef]; \)\)

\(Print[
      ColumnForm[

        	\ \ {"Electric fields (u-component) on field lines as", \ 

          \ \ 	\ \ "functions of v (colours match the colours above)"}]]
      \)\(\(efuatefl11\  = \ 
      ParametricPlot[\ 
        Evaluate[
          Table[
		{
              Re[cm11[ .9999 \((k - 1)\)/\((nef - 1)\)\ EllipticK[mprime]\ I
                      \  + x]], 
		\ 
              Im[ef11[ .9999 \((k - 1)\)/\((nef - 1)\)\ EllipticK[mprime]\ I
                      \  + x]]}, 
		{k, nef, 1, \(-1\)}]], {x, 0, 
          EllipticK[m] -  .0001}, \ 
		
        PlotRegion\  -> \ {{ .01,  .99}, { .01,  .99}}, 	
		
        AxesLabel\  -> \ {"v [mm]", \ \*""\!\(E\_u\)[V/m]""}, 
		
        PlotStyle\  -> \ stylef]; \)\)

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

up = 3; vp = 4; 
(* 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. eflandeqpl11 *) 

If[vp < 0, Print[vp" = vp < 0  is unphysical."], 
   
frefp = FindRoot[cm11[efp]-vp-I Abs[up],
        {efp,.7EllipticK[m]+.7I EllipticK[mprime],
             .3EllipticK[m]+.3I EllipticK[mprime]},
		    MaxIterations->50];
	
If[Sign[up]<0,efp = Conjugate[efp]/.frefp,
              efp = efp /. frefp ]; 
imefp = N[Im[efp]]; reefp = N[Re[efp]];
error = Abs[N[cm11[efp] - (vp + I up)]]; 

Which[!(0 < reefp < EllipticK[m]), 	 
 Print["The calculated point doesn`t fulfill the requirements!",\[FreakedSmiley]], 
      !(-EllipticK[mprime]< imefp <EllipticK[mprime]), 
 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 = ef11[efp]//N;
Print[
SequenceForm["At point (up,vp) = (",up ,", ",vp ,") mm : ",
"Potential = ",N[V Re[efp]/(EllipticK[m])], " V;"]];
Print[SequenceForm["Field (\!\(E\_u\),\!\(E\_v\)) = (", 
N[Im[valef]],", ", N[Re[valef]],") V/m "]]; 
 
caleqpl11 = ParametricPlot[Evaluate[
      {Im[cm11[Re[efp]+I s]],Re[cm11[Re[efp]+I s]]}],    
      {s,-.9999EllipticK[mprime],.9999EllipticK[mprime]}, 
      PlotStyle -> Dashing[{0.01}], 
      DisplayFunction -> Identity]; 
      
calefl11 = ParametricPlot[Evaluate[
      {Im[cm11[Im[efp]I+s]],Re[cm11[Im[efp]I+s]]}], 
      {s,0,.9999EllipticK[m]}, 
      DisplayFunction -> Identity]; 
      
Show[caleqpl11,calefl11, 
     DisplayFunction -> $DisplayFunction, 
     AxesLabel -> {"u [mm]", "v [mm]"},
     PlotRange -> All,
     Epilog -> {Thickness[tcb],
                Line[{{-b/2,a/2},{b/2,a/2}}],
                Line[{{-5b,0},{5b,0}}]}]; ]]
                
On[General::spell1]