Plot electric potential and field

Question

I am trying to plot the electric potential of a dielectric cylinder along with its field. The potential is a piecewise function:

F[r_, f_] := 
 Piecewise[{{-2 e r Cos[f]/(1 + er), 
0 <= r < a}, {-e r Cos[f] + (er - 1)*a^2 *e*r^(-1)*Cos[f], 
r >= a}}](*electorstatic pottential*)
e = 500;(*outer electric field*)
a = 0.02;(*cylinder's radius*)
er = 2;(*relevant acceptance*)

What's the best way to plot this function (potential), with its derivative (electric field)?

I tried to plot it using ContourPlot, but it doesn't look nice (I am not a Mathematica expert apparently)

ContourPlot[-2 e r Cos[f]/(1 + er), {r, 0, a}, {f, 0, 2 Pi}]
ContourPlot[-e r Cos[f] + (er - 1)*a^2 *e*r^(-1)*Cos[f], {r, a, 
2 a}, {f, 0, 2 Pi}]
Show[%, %%]

What I am trying to achieve is a sophisticated plot, with the cylinder at the origin, where the equipotential lines and the electric field lines will be plotted. I found something like that on the net (specifically the page "Gradient field plots in Mathematica"), but I don't know how to modify it...

gradientFieldPlot[f_, rx_, ry_, opts : OptionsPattern[]] := 
 Module[{img, cont, densityOptions, contourOptions, frameOptions, 
   gradField, field, plotRangeRule, rangeCoords}, 
  densityOptions = 
   Join[FilterRules[{opts}, 
     FilterRules[Options[DensityPlot], 
      Except[{Prolog, Epilog, FrameTicks, PlotLabel, ImagePadding, 
        GridLines, Mesh, AspectRatio, PlotRangePadding, Frame, 
        Axes}]]], {PlotRangePadding -> None, Frame -> None, 
     Axes -> None, AspectRatio -> Automatic}];
  contourOptions = 
   Join[FilterRules[{opts}, 
     FilterRules[Options[ContourPlot], 
      Except[{Prolog, Epilog, FrameTicks, PlotLabel, Background, 
        ContourShading, PlotRangePadding, Frame, Axes, 
        ExclusionsStyle}]]], {PlotRangePadding -> None, Frame -> None,
      Axes -> None, ContourShading -> False}];
  gradField = ComplexExpand[{D[f, rx[[1]]], D[f, ry[[1]]]}];
  field = 
   DensityPlot[Norm[gradField], rx, ry, 
    Evaluate@Apply[Sequence, densityOptions]];
  img = Rasterize[field, "Image"];
  plotRangeRule = FilterRules[Quiet@AbsoluteOptions[field], PlotRange];
  cont = If[
    MemberQ[{0, 
      None}, (Contours /. FilterRules[{opts}, Contours])], {}, 
    ContourPlot[f, rx, ry, Evaluate@Apply[Sequence, contourOptions]]];
  frameOptions = 
   Join[FilterRules[{opts}, 
     FilterRules[Options[Graphics], 
      Except[{PlotRangeClipping, PlotRange}]]], {plotRangeRule, 
     Frame -> True, PlotRangeClipping -> True}];
  rangeCoords = Transpose[PlotRange /. plotRangeRule];
  Apply[Show[
     Graphics[{Inset[
        Show[SetAlphaChannel[img, 
          "ShadingOpacity" /. {opts} /. {"ShadingOpacity" -> 1}], 
         AspectRatio -> Full], rangeCoords[[1]], {0, 0}, 
        rangeCoords[[2]] - rangeCoords[[1]]]}], cont, 
     StreamPlot[gradField, rx, ry, 
      Evaluate@FilterRules[{opts}, StreamStyle], 
      Evaluate@FilterRules[{opts}, StreamColorFunction], 
      Evaluate@FilterRules[{opts}, StreamColorFunctionScaling], 
      Evaluate@FilterRules[{opts}, StreamPoints], 
      Evaluate@FilterRules[{opts}, StreamScale]], ##] &, 
   frameOptions]]

This can be run like that

gradientFieldPlot[(y^2 + (x - 2)^2)^(-1/
    2) - (y^2 + (x - 1/2)^2)^(-1/2)/2, {x, -1.5, 2.5}, {y, -1.5, 
  1.5}, PlotPoints -> 50, ColorFunction -> "BlueGreenYellow", 
 Contours -> 10, ContourStyle -> White, Frame -> True, 
 FrameLabel -> {"x", "y"}, ClippingStyle -> Automatic, Axes -> True, 
 StreamStyle -> Orange]

And the amazing output is

What I would like to do is

Answer 1

I'll just address the technical plotting aspect of the question.

The problem is that your cylinder potential is intended to be a function of polar coordinates r and f, but to plot it in Cartesian space you have to convert to x and y first. This is done as follows:

pot[r_, f_] := 
 Piecewise[{{-2 e r Cos[f]/(1 + er), 
    0 <= r < a}, {-e r Cos[f] + (er - 1)*a^2*e*r^(-1)*Cos[f], 
    r >= a}}](*electorstatic pottential*)
e = 500;(*outer electric \
field*)a = 0.02;(*cylinder's radius*)er = 2;

gradientFieldPlot[
 pot[Sqrt[x^2 + y^2], ArcTan[x, y]], {x, -.2, .2}, {y, -.2, .2}, 
 PlotPoints -> 50, ColorFunction -> "BlueGreenYellow", Contours -> 10,
  ContourStyle -> White, Frame -> True, FrameLabel -> {"x", "y"}, 
 ClippingStyle -> Automatic, Axes -> True, StreamStyle -> Orange]

The colors represent the field strength, and the orange lines are the field lines.

I called the potential pot instead of F.

Since your cylinder is so small, we should zoom in even more:

gradientFieldPlot[
 pot[Sqrt[x^2 + y^2], ArcTan[x, y]], {x, -.05, .05}, {y, -.05, .05}, 
 PlotPoints -> 50, ColorFunction -> "BlueGreenYellow", Contours -> 10,
  ContourStyle -> White, Frame -> True, FrameLabel -> {"x", "y"}, 
 ClippingStyle -> Automatic, Axes -> True, StreamStyle -> Orange]

So everything looks as expected for a dielectric cylinder in the homogeneous external field of a capacitor. The equipotential lines (white) are completely flat inside the cylinder, corresponding to a uniform field. The highest field strength is at the "poles" of the dielectric cylinder facing the (far-away) capacitor plates.